Maxime Gobert
email: maxime.gobert@umons.ac.be

Nouredine Melab
email: nouredine.melab@univ-lille.fr

Jan Gmys
email: jan.gmys@inria.fr

Daniel Tuyttens
email: daniel.tuyttens@umons.ac.be

Adaptive Space Partitioning for Parallel Bayesian Optimization

Keywords: Parallel Bayesian Optimization, Design Space partitioning, Acquisition Process, Black-Box Optimization

This paper presents a new approach to provide large and well-chosen batches of points to evaluate in parallel in the context of batch-parallel Global Optimization. The method combines Bayesian Optimization and Design Space Partitioning in order to independently select different candidates according to the expected information they could provide if included in the design of experiment. The information in question is given by a figure of merit computed thanks to the chosen bayesian model. The presented algorithm is designed to be flexible and adaptable to various objective functions since it automatically adapts the space partition during the optimization process. Three benchmark functions are investigated in this study and reveal mixed results, but also interesting features for future works.

INTRODUCTION

This paper deals with Global Optimization and more precisely Black-Box Optimization, min D f (x), where D is a box constrained domain. The term Black-Box refers to the fact that the objective function is not analytically known and therefore neither its gradient, no information is known a priori. Furthermore, it is assumed that the evaluation of the objective function comes at a high simulation cost, as it is typically the case when f (x) is the outcome of a numerical simulation for an input vector x ∈ D. A classical approach for dealing with computationally expensive objective functions is to build a cheaper-to-evaluate approximation of the objective function called surrogate (model) or metamodel. The basic principle of surrogate-assisted optimization consists in first building an initial surrogate model using few points obtained with some sampling technique (e.g. Latin Hypercube Sampling). The surrogate is used to determine promising candidate solutions to evaluate with the (real) objective function which, in turn, are used to iteratively refine the surrogate.

Our focus is set on Bayesian Optimization (BO) which has proved to be very efficient in this context. The original concept is introduced by Kushner et al. [START_REF] Kushner | A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise[END_REF] and generalized by Mockus et al. [START_REF] Močkus | On bayesian methods for seeking the extremum[END_REF]. In BO, a prior distribution is set over the sample points, which is commonly a gaussian prior. After the evaluation of the sample, the prior is updated as a posterior distribution used to compute a figure of merit that is designed to provide a valuable candidate through its optimization. This figure of merit is computed thanks to the Acquisition Function (AF), also named Infill Criterion (IC) and helps finding the areas to sample in the landscape to be optimized.

Facing time-consuming objective functions, parallel computing can be used as a complementary way, in addition to surrogate models, to further speed up the optimization process. However, the joint usage of parallel processing and surrogatebased optimization is challenging, as it requires to be able to provide a well-chosen batch of candidates for parallel evaluation. We call Acquisition Process (AP) the way to provide one or multiple candidates. It is usually constituted of one or several AFs and their respective optimization. Finally, we define a cycle as the following sequence of instructions: update the surrogate model, perform AP and evaluate the selected candidate(s) simultaneously.

The objective of this paper is to present a new AP based on space partitioning optimization in order to efficiently select a large and valuable batch of candidates from the surrogate model that will be exactly evaluated in parallel. After presenting related works and precisely defining the objectives of the new approach in Section I, the latter is precisely described in Section II. Then, a benchmark analysis is presented in Section III before concluding on benefits of the approach and future works.

I. RELATED WORKS AND OBJECTIVES

Efficient Global Optimization (EGO) is developed by Jones et al. in [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF]. It uses Gaussian Processes (GP) surrogate models, which has the ability to provide a measure of uncertainty on the prediction. The Expected Improvement (EI) is used as AF and is designed to find a trade-off between exploring the regions where the surrogate accuracy is low, and exploiting in area where predicted values are close to the best known solution.

A. Related Works

Dealing with parallel BO, one major difficulty is to provide useful candidates to improve either the surrogate model, either the objective value. BO and in particular EGO provide a way by using an AF. However this process is inherently sequential and must be slightly modified to exploit parallel computing. In [START_REF] Carraro | A multi-points criterion for deterministic parallel global optimization based on gaussian processes[END_REF] Ginsbourger et al. several approaches for the parallelization of EGO are presented. To select q candidates, the analytical q-points EI (q-EI) is defined and optimized through time-consuming Monte Carlo simulations. An alternative is proposed to emulate q-EI: each time a candidate is provided by the optimization of the AF, the model is updated with an arbitrary chosen objective value attributed to the candidate (i.e. without any access to the real simulator). In this multi-point AP, q points are sequentially selected by repeating the singlepoint AP. Due to the relevance of the q-EGO algorithm for the present work, it will be further detailed in Section III-A1.

Motivated by a wider access to parallel machines, several other multi-point APs have been designed. Wang et al. [START_REF] Wang | Parallel bayesian global optimization of expensive functions[END_REF] proposed another way to select q points using infinitesimal perturbation analysis to construct a gradient estimator of the multi-point EI surface coupled with a multi-start procedure to find the set of points to evaluate exactly. The method appears to be fast while providing the bayes-optimal set of points. Another method from Marmin et al. [START_REF] Marmin | Differentiating the multipoint expected improvement for optimal batch design[END_REF] introduces the analytical form of the multi-point EI gradient to be able to optimize the function with gradient information. This approach reduces the computational cost of q-points selection compared to sequential heuristics or Monte Carlo sampling of [START_REF] Carraro | A multi-points criterion for deterministic parallel global optimization based on gaussian processes[END_REF].

Nevertheless, it remains costly in high dimensions because it is still a d * q dimensional problem, d being the number of decision variables. Based on the same idea of providing a batch of points at each cycle, Kandasamy et al. [START_REF] Kandasamy | Asynchronous parallel bayesian optimisation via thompson sampling[END_REF] use Thompson Sampling to maximize the probability of reward (i.e. improving the target). Shah et al. [START_REF] Shah | Parallel predictive entropy search for batch global optimization of expensive objective functions[END_REF] extend the Entropy Search from Hennig et al. [START_REF] Hennig | Entropy search for informationefficient global optimization[END_REF] for parallel computation. In [START_REF] Wang | Multi-point efficient global optimization using niching evolution strategy[END_REF] and [START_REF] Dawei Zhan | Balancing global and local search in parallel efficient global optimization algorithms[END_REF] niching Genetic Algorithms are used to select simultaneously ensemble of points. Another approach consists in using different complementary AFs or ensembles of surrogate models. Wang et al. [START_REF] Wang | Efficient global optimization using multiple infill sampling criteria and surrogate models[END_REF] uses n infill criteria coupled to multipoint proposal inspired by the Kriging Believer heuristic from [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF] to get n * q candidates per cycle. Lyu et al. [START_REF] Lyu | Batch Bayesian optimization via multi-objective acquisition ensemble for automated analog circuit design[END_REF] and Feng et al. [START_REF] Feng | A multiobjective optimization based framework to balance the global exploration and local exploitation in expensive optimization[END_REF] use multi-objective optimization on different AFs to get a Pareto set of optimal candidates.

B. Objectives

Many methods have been developed trying to either associate different models or AFs to propose several points, either optimizing the multi-point version or a single criterion. However, APs based on partitioning the decision space have only been sparsely addressed in the literature. We can mention Villanueva et al. [START_REF] Villanueva | Dynamic design space partitioning for optimization of an integrated thermal protection system[END_REF] in which a multi-agent strategy is used to explore different subdomains independently, following the machine learning formalism of reinforcement learning. Agents that performs well are rewarded and possibly split while other may merge. Wang et al. uses clustering techniques in [START_REF] Wang | Fuzzy clustering based hierarchical metamodeling for design space reduction and optimization[END_REF] to locate distinct local optima.

Parallel versions of EGO-like algorithms often suffer from a time-consuming AP or costly metamodel updates, and a major difficulty stays in balancing the optimization process [START_REF] Briffoteaux | Parallel surrogate-assisted optimization: Batched bayesian neural networkassisted ga versus q-ego[END_REF]. Trying to push these limits back, we propose a new algorithm which is able to select large batches of points to evaluate in parallel in a moderate amount of time, while keeping a balance between exploration and exploitation. We use a space partition managed by a self-organizing binary tree in order to perform simultaneously different local APs in each subdomain. The Binary Space Partitioning tree (BSP tree) allows to decompose the global AP in several ones. It also balances the optimisation process using a decision heuristic dealing with where to intensify the decomposition process, and where to sample less frequently. The strength of this method lies in its adaptability. It is able to provide as many candidates as needed, remains fast to execute and parallelizable (sub-APs are independent from each other). Finally, the partition tree adapts automatically according to the rule explained in Section II-B2.

II. PRESENTATION OF THE METHOD

As for many BO algorithm, the method rests on the Expected Improvement (EI) criterion. To be computed, EI requires that the metamodel provides a value of the uncertainty of the prediction. Thanks to its statistical nature, the Kriging equations provide both prediction of the output at the desired location and its variance.

A. Metamodel

Kriging intends to predict the value of a certain quantity according to its location. It uses spatial correlation between the locations and their measure of interest (the output) to build the Best Linear Unbiased Predictor. The theoretical basis of the Kriging model will not be presented here, but can be found for example in [START_REF] Noel | Statistics for spatial data[END_REF]. Nevertheless some parameters need to be defined for better understanding. Let us call x = (x 1 , . . . , x d ) ∈ D ⊂ R d previously mentioned as the location, and y(x) ∈ R the objective value at this specific location. The value is assumed to be a random variable that writes as the sum of a trend function µ(x) and a centered gaussian process z(x): y(x) = µ(x) + z(x). Denoting y(x (i) ) = y (i) , we write y = (y (1) , . . . , y (n) ) the known values of the outputs at locations X = (x (1) , . . . , x (n) ). For this study, the trend is chosen to be constant, this kind of kriging is also known as Ordinary Kriging. The spatial correlation is expressed between two points using a kernel function that needs to be defined as an hyper-parameter of the model.

The Kriging prediction is the Best Linear Unbiased Predictor, hence the prediction of the output ŷ(x * ) = ŷ * at

B. Acquisition Process

The acquisition process is dedicated to find a batch of points to simultaneously evaluate with the objective function (e.g. the simulator). The particularity of our AP is that it is decomposed in several sub-processes responsible for providing one candidate each.

1) Binary Space Partition: In the proposed AP, the search space D is recursively decomposed according to a binary separation of space. Let us suppose that the whole domain is D 1 . The partitions are managed by a binary tree where the root node (at tree level 1) contains D 1 . The next level (level 2) has two nodes (node 2 and node 3) that contain D 2 and D 3 such that D 1 = D 2 ∪ D 3 . For each node k, the property is respected such that D k = D 2k ∪ D 2k+1 and thus for each level of the tree, the union preserves the entire domain.

An example is presented in Fig. 1, where the domain is split in four subdomains. Let us call F n the family of subdomain indexes at cycle n, such that k∈Fn D k = D and ∀i, j ∈ F n D i ∩ D j = ∅. Regarding Fig. 1, F n = {2, 6, 7} and F n = {3, 4, 5} are acceptable sets. Thanks to this kind of decomposition, it is easy to perform one sub-AP in each subdomain, while keeping knowledge of the entire domain.

As soon as all the candidates are computed, they are collected and sorted according to the chosen figure of merit. A subset of them will be exactly evaluated while the rest is discarded. We call the number of candidates evaluated in parallel the batch size n batch .

This strategy intends to reinforce the global aspect of optimization by keep sampling in a priori less interesting areas of D (from the IC point of view). Nevertheless, in order to avoid sampling with clearly no gain, and thus waste computational power, the total number of candidates will be greater than the batch size. Furthermore, to balance the exploration and exploitation process, the tree has the possibility to evolve and split further some nodes resulting in intensifying search into the best subdomain -always in terms of the chosen figure of merit.

Even though the subdomains are distinct, it may happen that several candidates are very close to a shared boundary, and thus to each other. In that scenario, the candidates receive a small random perturbation so that the area is still sampled twice, but most of all the stability of the kriging model is preserved.

2) Evolution strategy: For this study, the number of candidates provided by the AP before the selection is a multiple of the batch size. Indeed, the batch size is fixed equal to the number of available cores (1 evaluation per core), thus each computing unit performs the same number of sub-APs and the parallel load is balanced. As stated in the previous paragraph, one candidate is chosen in each subdomain, thus we have as many candidates as leaves in the tree. Let us call n leaves that number, and n batch the batch size. Consequently n leaves = r * n batch where r ∈ N \{0} .

The tree is updated once a cycle to take into account the new information. The supposed best subdomain, according to the IC, is decomposed further to intensify the search in this area. Nevertheless, as we decide to keep n leaves constant, this splitting step is only performed if two domains are merged. In terms of the BSP-tree, the leaf with the highest figure of merit is split, and the parent node with the lowest one loses its leaves to become a leaf itself. This process is illustrated in Fig. 2.

This example directly follows the one from Fig. 1, one candidate is chosen in each leaf (i.e the ones indexed by F 0 = {4, 5, 6, 7}). Each node is attributed the best figure of merit of its children, this number is denoted in Fig. 2 by IC node . In Fig. 2a, D 4 possesses the best value among the nodes indexed by F 0 . Consequently, the node will be split if model ← learn model(X, y)

D 1 root D 2 IC 2 = best(IC 4 , IC 5 ) D 4 IC 4 D 5 IC 5 D 3 IC 3 = best(IC 6 , IC 7 ) D 6 IC 6 D 7 IC 7 (a)
6:

B candidates ← ∅ 7:
for leaf in F i do parallelizable loop 8:

c ← find best candidate(D leaf )

9:

B candidates ← B candidates ∪ c 10:
end for 11:

B candidates ← selection(B candidates )

12:

(T , F i+1 ) ← update tree(T , B candidates , F i )

13:

X ← X ∪ B candidates 14:

y ← y ∪ f (B candidates ) parallel evaluation 15: end for 16: return x max the merge operation can be performed. Regarding the parents of the leaves, D 3 possesses the worst value, meaning that area does not need as many attention and thus it will be merged. Eventually, F 1 = {3, 5, 8, 9}.

In case of a non-allowed operation, the tree is kept identical for the next cycle. For instance, regarding Fig. 2b, it may happen that for the next cycle D 3 is still the worst subdomain and can't be merged with D 2 . However, this kind of exception is relatively rare when dealing with large trees.

3) Splitting the subdomain: Without loss of generality, let us set D 1 = [0, 1] d . We must decide how to separate the hypercube. For practical reasons, the choice is made for now to keep hyper-rectangular domains, so that a splitting operation is characterized by the axis/dimension to be split, and the range at which the section is done. For example, if D 1 = [0, 1] 2 is split according to the first axis in the middle of the segment, it comes to

D 2 = [0, 0.5] × [0, 1] and D 3 = [0.5, 1] × [0, 1].
The heuristic for this study is arbitrary splitting: based on the idea that it is preferable to have dimensions of the same order of magnitude, axes are split one after another in a cyclic way. The chosen axis is determined by the level of depth of the tree node, i.e. the axis along which a subdomain is split is given by: d split = depth(node) mod(d). The initial tree is formed using this heuristic until it reaches the desired depth.

C. BSP-EGO

Named after EGO algorithm, this methods is called BSP-EGO. It uses a global kriging metamodel and the presented BSP-based AP with one sub-AP per sub-region. BSP-EGO is outlined in Algorithm 1. First, the initial sample is created, namely the couple (X, y). Then, the tree T is initialized at a predefined depth, deduced from the user-defined hyperparameter n leaves . At the beginning of each cycle, starting at line 4 of Algorithm 1, the model is updated and then, the AP begins in line 7. In each leaf of the tree, marked in F i , a candidate is proposed using classical optimization on the chosen AF. These AF maximisations are independent and can be performed in parallel. Candidates are gathered and the n batch most promising ones, according to the EI criterion, are selected (line 11). Then the tree and the leaves (i.e. their indexes) are updated. A cycle ends with the parallel evaluation of the selected candidates and their insertion into the database (lines 13 and 14). BSP-EGO is implemented in C + +, compiled with GCC 8.3.0, using OpenMPI 3.1 for the parallel evaluation of candidate solutions and the BayesOpt C++ library [START_REF] Martinez-Cantin | Bayesopt: A bayesian optimization library for nonlinear optimization, experimental design and bandits[END_REF] for the Kriging metamodel.

III. BENCHMARK ANALYSIS

A. Experimental protocol

In this section, experimental results obtained with BSP-EGO are presented and compared to the q-EGO algorithm from Ginsbourger et al. [START_REF] Carraro | A multi-points criterion for deterministic parallel global optimization based on gaussian processes[END_REF]. One implementation is available in the R package DiceOptim presented in Roustant et al. [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]. All experimentation are performed using the Chetemi cluster from Grid5000 [START_REF] Balouek | Adding virtualization capabilities to the Grid'5000 testbed[END_REF]. The cluster is composed of 15 nodes with 2 Intel Xeon E5-2630 of 10 cores each. Among the 300 available cores, only n batch are reserved for experiments.

1) Reference approach: The used q-EGO algorithm is outlined in Algorithm 2 and differs from BSP-EGO in its AP. model ← learn model(X, y)

Algorithm 2 q-EGO pseudo-code

5:

B candidates ← ∅ 6:
X tmp ← X 7:

y tmp ← y 8:
for k in 1 : q do sequential loop 9:

x (k) ← argmax D (IC(x))

10:

B candidates ← x (k)
11:

y (k) ← best value (y) 12: X tmp ← X tmp ∪ x (k)
13: x 2 i -exp

y tmp ← y ∪ y (k)
1 n n i=1 cos(2πx i ) + 20 + exp(1) 0 Rosenbrock [-32, 32] 6 d i=1 100(x i+1 -x 2 i ) 2 + (1 -x i ) 2 0 Alpine [0, 10] 6 d i=1 √ x i sin(x i ) -490.21
To be able to provide a batch of distinct candidates without accessing the exact evaluation, q-EGO iteratively performs q maximisations of the AF, each followed by a metamodel update with heuristically assigned objective values. Indeed, assigning a value to a candidate allows to update the model and transform the landscape of the infill criterion so that it knows that the area is sampled (or will be), resulting in a low figure of merit for this location. The used heuristic is named Constant Liar (CL). It assigns to all candidates the best current objective value. As shown in lines 8 to 15 of Algorithm 2 this process is repeated q times until the batch B candidates contains the q = n batch candidates. As a consequence, this method involves q sequential updates of the model (line 14), and a temporary learning set denoted as (X tmp , y tmp ) in Algorithm 2.

Several alternatives to CL are presented in Ginsbourger et al. [START_REF] Ginsbourger | Kriging is well-suited to parallelize optimization[END_REF] such as Kriging Believer where the model is trusted and the lie of CL is replaced with the prediction of the model. According to preliminary experiments, CL seems to perform better on our benchmark test-bed. The creation of the model as well as the optimization process are left to the DiceKriging package [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF] concerning the q-EGO part.

2) Hyper-parameters settings: For both approaches, the metamodel is an Ordinary Kriging model -i.e with constant trend -with Matern 5 2 covariance kernel. The study is performed with n batch = 8, 16, 32.

Concerning BSP-EGO specific parameters, the initial tree is build with 2 * n batch leaves and this number is preserved all through the 48 cycles attributed to each experimentation. Only one tree update is allowed per cycle.

For each experimentation, we compare the average over 10 repetitions of the objective value evolution over the 48 cycles. The comparisons are done using the same initial sample containing 64 points generated randomly so that each approach starts with the same best known value at cycle 0.

3) Benchmark functions: The performance analysis is performed on three classical benchmark functions. They are chosen to be representative of current challenges in global optimization. The first one is the Ackley function: it is optimized in the domain D = [-32, 32] 6 and its only global minimum is 0. The Ackley function is also known as the "well-function" since it has (in its 2D landscape) a clear minimum in a well in the middle of a noisy plane. The noise can be difficult to manage when dealing with surrogate models and IC. Second choice is the Rosenbrock function, optimized as well in D = [-32, 32] 6 with 0 as unique global minimum. Also named Valley or Banana function because of its 2D representation, the function is unimodal, and the global minimum lies in the valley. Even though the valley is easy to spot, algorithms may have convergence problems. The last one for this study is the Alpine function, named after its multimodal 2D landscape. It is optimized in D = [0, 10] 6 and has a single global minimum of -2.808 d ≈ -490.21. Its difficulty lies in the number of local minima of same amplitude in which the process may be stuck. These functions are summarized in Table I.

B. Experimental Results

Fig. 3, 4 and5 show the evolution of the best target (y min ) in unction of the current cycle for the three test-functions Ackley, Rosenbrock and Alpine respectively. The black-dotted curves represent the results of the q-EGO reference approach. The results from BSP-EGO approach are labeled as redcrossed curves. Additionally, Fig. 6, 7, 8 show boxplots of the best objective found after 48 cycles.

Fig. 3 shows that the space partitioning approach presents a good behavior regarding the results from the Ackley benchmark function. For each batch size (n batch ∈ {8, 16, 32}) the BSP-EGO curve is beneath the q-EGO one, meaning that at any cycle of the optimization, BSP-EGO provides on average a better solution. This is confirmed by the boxplot in Fig. 6 where we can see that almost any execution of BSP-EGO performs better than every q-EGO one. Furthermore, adding parallelism allows to provide a better final value for the same number of cycles. As shown in Fig. 6, while q-EGO stagnates, the BSP approach continues to improve the objective value as n batch increases.

For sake of readability, the figures associated with the Rosenbrock function do not display the entire evolution of the objective value. Indeed the display window is focused on the final values to be able to tell the difference between the curves. Thus, the high output values are not visible for the first cycles.

Fig. 4 displays the results for the Rosenbrock function and shows that the q-EGO curves are always on top of BSP-EGO ones. Even though q-EGO performs better on this benchmark function, it has to be noted that the gap shrinks when increasing the batch size. As for the Ackley function, increasing the batch size and thus the degree of parallelism favors BSP-EGO and enables the latter to find better solutions. The evolution of the best found solution with n batch is clearly Regarding our last benchmark function, Alpine, we can observe in Fig. 5 that BSP-EGO performs better than q-EGO for n batch = 8, 16, while for a batch size of 32 q-EGO outperforms its concurrent at the end of the optimization. Regarding the average curves, BSP-EGO appears to be faster than q-EGO to find a local minimum but has difficulties to escape from it. One hypothesis that could explain this behavior is that when increasing the batch size, the number of leaves is also increased, allowing bigger sampling in a relatively small region and therefore over-intensifying the search. Nevertheless, the boxplots in Fig. 8 indicate that the final results of the two approaches are very similar and they both managed to find the global minimum.

The recursive splitting allows to sample more areas of interest which seems beneficial for the Ackley function where classical BO might face difficulties to exploit the optimal area. Furthermore, increasing the degree of parallelism allows to improve the quality of the best found solution. This is especially visible regarding Fig. 4, which shows that BSP-EGO becomes better with wider parallelism. The two other benchmark functions also seem to benefit from the larger The study suggests that the proposed AP can perform as good as q-EGO while being much faster as it requires only one model update instead of q. However, it has to be noted that the cost of fitting the kriging model still becomes prohibitive with the size of the design of experiment which grows especially fast rising the batch size.

IV. CONCLUSION

BO and parallel computing are common ways to deal with the optimization of highly time-consuming function/simulators. However making them act in concert still raises difficulties. We propose an approach to provide n batch candidates based on a recursive design space decomposition managed by Fig. 8: Box-plots of the final value for BSP-EGO and q-EGO -Alpine function a binary decomposition tree. A performance analyse is led on three benchmark functions chosen to represent known difficulties in global optimization. The proposed method challenges the q-EGO approach from [START_REF] Roustant | DiceOptim: Two R packages for the analysis of computer experiments by kriging-based metamodeling and optimization[END_REF]. According to the presented results in Fig 3,[START_REF] Carraro | A multi-points criterion for deterministic parallel global optimization based on gaussian processes[END_REF][START_REF] Wang | Parallel bayesian global optimization of expensive functions[END_REF], the proposed partition-based AP seems to perform fairly compared to the state-of-the-art method q-EGO.

In addition, the developed AP presents interesting qualities such as its much smaller execution time, which can still be reduced using parallelism for the AP. In the current version, BSP-EGO only performs the exact evaluations in parallel. This AP allows to propose arbitrarily big batches of valuable points with only a small extra cost. BSP-EGO is presented in this paper in its first version, without hyper-parameter calibration. Especially parameters from the tree management need to be further investigated (e.g. number of updates or leaves). Heuristics are being developed based on correlation coefficients to determine the most significant axis to chose as split direction. Regarding the figures of merit coming from the n batch candidates we observe that the value of the last candidates are low. This observation is an argument to explore a multi-criteria approach based on complementary AF. Last, we are planning to investigate the presented BSP-approach with different surrogate models, or with a Kriging model fitted on a subset of points to tackle the prohibitive cost of fitting the Kriging model on large data samples and be able to increase the parallel potential of BSP-EGO. All these ideas remain to be exploited and should hopefully improve the newly presented AP.

7 Fig. 1 :

 71 Fig. 1: Partitioning of D through BSP tree

1 D 2 D 4 D 8 D 9 D 5 D 3 (b) Step 2 :Fig. 2 :

 124895322 Fig. 2: Illustration of one tree update

  Input f : function to optimize n init : initial sample size n cycle : number of cycles D: search space 1: X ← initial sampling(D, n init ) 2: y ← f (X) 3: for i in 1 : n cycle do 4:

-n batch = 32 Fig. 3 :-n batch = 32 Fig. 4 :Fig. 5 :

 3233245 Fig. 3: Evolution of the best found target in function of the cycle -Ackley function

Fig. 6 :Fig. 7 :

 67 Fig.6: Box-plots of the final value for BSP-EGO and q-EGO -Ackley function

  Alpine -n batch = 32

TABLE I :

 I Characteristics of benchmark functions

	f	D		f (x) =	y min
	Ackley	[-32, 32] 6 -20.exp -0.2	1 n	d i=1

ACKNOWLEDGMENT

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several Universities as well as other organizations (see https://www.grid5000.fr).