Nacim Belkhir

Johann Dréo

Pierre Savéant
email: pierre.saveant@thalesgroup.commarc.schoenauer@inria.fr

Marc Schoenauer

Feature Based Algorithm Configuration: a Case Study with Differential Evolution

Keywords: Empirical Study, Black-box Continuous Optimization, Problem Features, Algorithm Configuration, Empirical Performance Model, Differential Evolution

Algorithm Configuration is still an intricate problem especially in the continuous black box optimization domain. This paper empirically investigates the relationship between continuous problem features (measuring different problem characteristics) and the best parameter configuration of a given stochastic algorithm over a bench of test functions -namely here, the original version of Differential Evolution over the BBOB test bench. This is achieved by learning an empirical performance model from the problem features and the algorithm parameters. This performance model can then be used to compute an empirical optimal parameter configuration from features values. The results show that reasonable performance models can indeed be learned, resulting in a better parameter configuration than a static parameter setting optimized for robustness over the test bench.

Introduction

Today, it is widely acknowledged that the quest of a universal black box optimization algorithm is vain, even if the No Free Lunch Theorem [START_REF] Wolpert | No free lunch theorems for optimization[END_REF] has been questioned in the continuous framework [START_REF] Auger | Continuous Lunches are Free plus the Design of Optimal Optimization Algorithms[END_REF]. However, many algorithms exist, more or less specific to different classes of optimization problems, and the new grail of optimizers has now turned toward Algorithm Selection, as formulated by Rice [START_REF] Rice | The algorithm selection problem[END_REF], or Algorithm Configuration, that can be considered as yet another (meta-)optimization problem [START_REF] Hoos | Programming by optimization[END_REF]. In both cases, the choice (of an algorithm, or of the parameters of a given algorithm) is made w.r.t. the user's preference, aka a performance criterion (e.g., quality of the solution obtained in a given CPU cost, the smallest CPU cost to reach a given solution quality, the probability to reach a given quality, with given thresholds, etc.).

A first approach to Algorithm Configuration is to optimize this performance criterion once and for all using a specific algorithm, e.g., SMAC [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF]. But this results in a single configuration, and even if several problems are used to compute the performance criterion, the generalization of the results to other problems might be problematic.

More recent approaches are based on a description of the objective function in some feature space, and try to learn a mapping from this feature space onto the space of parameter configurations of the algorithm at hand, based on examples of the behavior of several configurations on a training set of objective functions. And the most successful approach for learning such a mapping is to first learn an empirical model of the algorithm performance (that predicts the performance criterion for a given set of features and an algorithm configuration). When a new problem arises (i.e. a new set of features), finding the algorithm configuration that is predicted to have the best performance is then straightforward. This approach, initially proposed in [START_REF] Leyton-Brown | Learning the empirical hardness of optimization problems: The case of combinatorial auctions[END_REF], has demonstrated successful results in different combinatorial optimization domains [START_REF] Hutter | Performance prediction and automated tuning of randomized and parametric algorithms[END_REF][START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF][START_REF] Xu | Satzilla: portfolio-based algorithm selection for sat[END_REF].

In continuous domains, however, though several feature sets have been proposed [START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF], and successfully demonstrated to accurately classify problem instances [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF], only Algorithm Selection problems have actually been tackled [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF][START_REF] Muñoz | A meta-learning prediction model of algorithm performance for continuous optimization problems[END_REF].

The present work addresses the Algorithm Configuration problem for continuous domains, building an Empirical Performance Model (EPM) based on the problem features in continuous search spaces cited above. The approach is experimentally validated with the original version of Differential Evolution [START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF], that has few hyper-parameters, but is known to be highly sensitive to their setting. The set of objective functions used for this validation is the well-known BBOB test bench [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF].

The paper is organized as follows. Section 2 presents the general idea of featurebased Empirical Performance Model and subsequent algorithm configuration. Section 3 surveys the problem features in the case of continuous optimization. Section 4 then details the DE case study and the BBOB testbench used in this work. Section 5 describes the different experiments, and the corresponding results are detailed and discussed in Section 6, before the concluding Section 7.

Algorithm Configuration with an Empirical Performance Model

Context and notations The general context is the Black Box Optimization of objective functions f ∶ Ω ↦ IR. An algorithm A is given together with its control parameters θ ∈ Θ. We assume that the objective functions can be described by some features ψ ∈ Ψ. The goal of Instance-based Algorithm Configuration is to find automatically, for a given objective function f described by its features ψ f , the best possible configuration of A, i.e., values θ * f ∈ Θ such that running A with parameters θ * f on f leads to optimal performances w.r.t. a given performance measure ϕ.

Empirical Performance Model

The first step is to build an Empirical Performance Model (EPM) φ that approximate ϕ on Ψ × Θ. A is run to optimize different functions f i (described by their features ψ fi ) using different parameter configurations θ j . This allows to compute the exact values ϕ(ψ fi , θ j ) for different pairs (i, j)1 . The set of all ((ψ fi , θ j ), ϕ(ψ fi , θ j )) is a training set that can be used as input to any standard regression method to learn a model φ for ϕ. Note that building such a model is done once, and hence its computer cost is not a critical issue.

Empirical Optimal Configuration When a new objective function g is to be optimized with A, its features ψ g are computed, and the optimization of φ(ψ g , θ g ) on parameter space Θ leads to the empirical optimal parameters of A for g. Here, the cost of computing the features ψ g , in terms of number of calls to g, is here of utter importance. In particular, it should be compared to the cost of running a full 'ad hoc' meta-optimization of A parameters for g (using e.g., SMAC [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF]).

The remaining of the paper is concerned with objective functions defined on some continuous domain D ⊂ IR d for a given dimension d ∈ IN. Different features, taken from the literature, will first be discussed, before the case study is detailed.

Features for Continuous Optimization

In this section, a single objective function f ∶ IR d ↦ IR is considered, and the feature space Ψ is a vector space of p real-valued features. The black-box context implies that features should be computed from samples of f2 , i.e. n pairs (x i , f (x i ))), specifically gathered for that purpose. The set of values {f (x i ) i = 1, . . . , n} is denoted Y.

A first set of 55 features is taken from [START_REF] Mersmann | Exploratory landscape analysis[END_REF]. These features are grouped into six classes: the 3 y-Distribution features are related to the distribution of the values in Y, the 18 Levelset features to the relative position of Y w.r.t a given threshold, the 9 Meta-Model features rely on meta-modeling of the sample set w.r.t linear and quadratic regression models, the 14 Curvature features on some numerical estimation of the Hessian and gradient of the problem, the 4 Convexity features on the empirical probability of convexity, and the 7 Local Search features on the ratio of local optima and global optima, estimated using some iterated local search procedure. The y-Distribution, the Levelset and the Meta-Model features can all be evaluated on the same sample dataset, hence their cost altogether is n, the number of samples. However, some additional evaluations are required for the other feature classes, that depend on the previous samples. The orders of magnitudes are about 10 3 × d for the Convexity features, around 10 4 × d for the Curvature and the Local Search features.

A set of 16 Dispersion features was originally proposed in [START_REF] Lunacek | The dispersion metric and the cma evolution strategy[END_REF]. They are based on comparisons of the distances between best samples from different percentiles of the overall sample (in terms of solution quality) to the mean or median distance between all samples. Finally, 5 Information Content features were proposed in [START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF], giving information about the global structure of the landscape.

Recent works [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF] successfully demonstrated that these features could be used in order to classify the optimization problems w.r.t their classes in BBOB (that will be introduced in Section 5) for the Algorithm Selection Problem. More recently, in [START_REF] Bossek | Learning feature-parameter mappings for parameter tuning via the profile expected improvement[END_REF] a subset of these features were used in order to improve the process of a parameter tuning algorithm, relying on the SMBO method [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF].

A Case Study in Continuous Domain: DE on BBOB

Differential Evolution and its Parameters Differential Evolution [START_REF] Storn | Differential evolution-a simple and efficient heuristic for global optimization over continuous spaces[END_REF] is a popular continuous optimization algorithm that encountered many successes. It is also known for its simplicity, at least in its original version, that comes at the price of a large sensitivity to its parameter setting: this is the reason why it has been chosen here, making it easier to see big differences of results for different parameter settings. Several advanced versions of DE exist, that clearly outperform the original version, but comparing our results with theirs is left for further work.

DE generates new individuals from the current population by adding to each individual in turn a difference vector between two other individuals, and recombining the result with another individual from the population. The original version of DE has only four static parameters:

• the population size NP ∈ IN;

• the strategy S ∈ {best1bin, randtobest1bin, best2bin, rand2bin, rand1bin} controls how to choose the endpoints of the difference vector; • F ∈ [0, 2] controls the intensity of the difference vector;

• the crossover rate CR ∈ [0, 1].
In this work the population size NP is kept to the default value 15×d recommended by the authors 3 . Note that the recommendation for the other parameters is S = best1bin, F = 0.8, and CR = 0.9, and will be used as one of the baseline (Section 6).

Test Bench

The following experiments consider the noiseless test functions from the Black Box Optimization Benchmark (BBOB4 ) [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF]. The BBOB test bench is made of 24 analytically defined functions defined on [-5, 5] d , with known global optima and known difficulties (e.g., non-separability, nulti-modality, etc). They have been manually classified in five classes of problems. In this work, only three of the BBOB classes are considered: 5 separable functions, 4 uni-modal functions with low or moderate conditioning, and 5 uni-modal functions with high conditioning (functions F1 to F14). Dimensions 2, 3, 5, 10 are considered for all functions. As advocated in the original framework, any independent run on a function is actually done on a variant, in order to get over a possible algorithm bias. Variants are obtained from the original function by a translation of the position of the optimum and -for the non-separable functions -a rotation of the coordinate system.

Performance Measure Following the COCO/BBOB framework, the performance measure used here is the Expected Run Time 5 (ERT) needed to reach the optimal objective value with a given precision. Let RT s be the average running time of successful runs, and p s the empirical probability of success (out of the 15 independent runs). The ERT is defined as ERT = ERT (f, θ) = RT s p s if the results were obtained with DE configuration θ optimizing test function f . Features From the features briefly introduced in Section 3, different set of problem features are considered. All features have been computed using the R package kindly publicly provided by Pascal Kerschke 6 . Features are distinguished by their costs:

ψ ⋆ includes all features from Section 3, with an initial sample of size k = 2000 × d. However, as discussed in Section 3, the actual cost is much larger because of features requiring additional evaluations. Regression Model Preliminary experiments, not discussed here, lead to consider a Random Forests regression model (in accordance with [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF]). A grid search with tenfold cross-validation has been run on the meta-parameters of the Random Forest. The implementation of the scikit-learn python library has been used throughout this work 7 , using 10 trees of maximal depth 200.

-ψ • k : k × d

Experiments

Dataset A 40-steps discretization is used for F ∈ [0, 2[ and CR ∈ [0, 1], resulting in 5 × 40 × 40 different configurations. For each of the 14 functions of the test bench and for each dimension d ∈ {2, 3, 5, 10}, each one of its 15 variants is optimized with these 8000 DE configurations, and the ERT is computed. The initial dataset is hence made of 14 × 8000 entries per dimension, or 440 000 entries in total, considering all dimensions.

Dimensions As discussed, the dimension d can be considered as a particular feature, available "for free" (without any function evaluation), or as part of the problem definition -and there are as many problems as dimensions. These two points of view will be compared here: the EPM will be learned either using only the entries of the dataset of the same dimension, or all entries, and the dimension will then be used as an additional feature in the feature vector.

Cross-validation

All experiments are based on a leave-one-out procedure: one of the 14 functions in the test bench is completely removed from the dataset (all dimensions and, of course, all variants). An EPM is then learned, and the left-out function considered a "unknown". The only exception is the robust baseline described below.

Baselines Different DE configurations are computed for each function, and used as a baseline for comparison with the results of the proposed approach. The default configuration recommended by DE authors (Section 4) is the first obvious baseline. However, it is likely to perform poorly across the whole test bench.

At the other extreme, the specific configuration found by some meta-optimizers for a given problem is likely to give the best overall results: two such meta-optimization were performed for each function: on the one hand, the best configuration encountered while computing the full dataset using the grid described above is saved; on the other hand, SMAC [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF] is applied to each function, using the ERT performance measure as fitness. The best of both configurations is reported as adhoc configuration.

Finally, one single SMAC optimization is performed using the average ERT over all 14 functions as performance: the idea behind this is to try to find some robust configuration that would give good results on all functions simultaneously. The resulting configuration is termed robust and considered as the reference (Section 6).

Experiment Costs All DE runs were allowed a maximum budget of 10 4 × d function evaluations -though of course some runs did stop earlier, having reached the target precision. And all results have been averaged over 15 independent variants for each function and dimension.

Furthermore, both adhoc baseline configurations have the same cost, because the budget given to SMAC was purposely chosen to match that of the grid search, i.e., 8000 × 15 runs of at most ×10 4 × dim evaluations each.

On the other hand, the robust configuration did cost 14 times more, as each of its iterations required to evaluate all the 14 functions -but has to be done only once.

Compared to that, the cost of finding the best empirical configuration using an EPM is the cost of the features: 500 or 2000 in the case of ψ • 500 and ψ • 2000 features, or around 2.10 4 × d for the full set of features ψ * .

Results

The series of experiments described above are presented in this section from two points of view: first, the different EPM are analyzed and compared to the ground truthand the Empirical Optimal Configuration is compared to the true optimal configuration in parameter space. Then, the actual results of DE optimization using the Empirical Optimal Configuration are compared to those of the different baselines, keeping in mind the actual costs of the different approaches (see last paragraph above).

EPM Analyses

Due to space constraints, only few typical figures are displayed8 (Figure 1) and will be discussed here. There are indeed some strong similarities between top and bottom colormaps for Figures 1a, 1c, and1d, even though they correspond to different functions, dimensions, feature sets, and dimension handling modes. On the other hand, the two plots for F4 (Figure 1b) are very different, and here the EPM fails to capture even an approximate shape of the true ERT landscape.

But beyond such comparisons, the optimal configurations of both plots are shown (on both plots too, to ease the comparison), displaying very different situations: in Figure 1a and Figure 1c, both optima are rather close (1c), or at least are both in the same color area of the true ERT; on the opposite, in Figure 1b are far from one another, and the Empirical Optimal Configuration lies in a region of very poor true ERT: the performance of these configurations used within DE for the corresponding function will be poor too (see forthcoming Section 6.2).

Empirical Optimal Configurations at Work

For each computed EPM (described in Section 5), an Empirical Optimal Configuration is obtained by optimization on the parameter space, and the ERT of this configuration is obtained by running DE on each of the 15 variants of the target function. Tables 1a and1b shows, for each function, the ratios of the ERTs of these different configurations against the robust configuration defined in Section 5 (the smaller the better). The first two columns are the other baselines, θ d are the parameters values recommended by DE authors, and θ L the adhoc configuration with best results (see Section 5). The two series of 3 columns correspond to the 3 different feature sets ψ • 500 , ψ • 2000 , and ψ ⋆ (as defined in Section 4) and the two types of dimension handling (learning performed only on the same dimension as testing, or on all dimensions at once -see Section 5).

As it could be expected, the adhoc configuration is a clear winner, and the default values recommended by DE authors a clear loser. The results of dimension 5 (as well as those in dimensions 2 and 3, not shown due to space constraints) bring several good news: most proposed approaches perform better than the robust configuration, and at g θ d θ L train on dim 5 train on all dims least one does, for all functions but F3. For some functions (F11, and also F8 and F10), some feature-based approaches even get close to the best adhoc configurations, never being worse than twice that best performance, except for F4 -as could be foreseen on Figure 1b. When it comes to compare the different EPM settings, learning only from the single target dimension gives better results than learning for all dimensions together -and in the former case, using all available features does improve over only using the cheap features. The situation is not so clear in dimension 10: in several cases, the Empirical Optimal Configuration cannot even reach the target in the allocated budget -a situation for which an example was given in Figure 1d. However, when an optimum can be found, similar conclusions to the dimension 5 case can be drawn, though not as contrasted.

ψ • 500 ψ • 2k ψ ⋆ ψ • 500 ψ • 2k ψ ⋆ F1
F3⋆ 100 0.2 ∞ ∞ 1 ∞ ∞ ∞ F4⋆ 10 0.3 ∞ ∞ 2 ∞ ∞ ∞ F5⋆ 10 0.2 ∞ ∞ ∞ ∞ ∞ ∞ F6

Discussion

Our results suggest that the learned EPM can be similar to the actual performance map, at least in a large part of the parameter space. Nonetheless, when a particular feature is not included in the learning set, it can be very hard for the EPM to achieve good accuracy and performance prediction, as witnessed with function F4 (Figure 1b): only F3 and F4 are multi-modal functions, and they have very different structures (apart from being separable, and hence belong to the first BBOB class). This was clear, too, with the F5 function (Figure 1d): F5 is the only linear function of the test bench, hence the EPM was learned without any linear function in the training set. However, the global accuracy of the EPM for F5 is rather good -not worse than F1 for instance (Figure 1a). But unfortunately, the small region of the parameter space where the EPM differs from the true ERT is the region that contains the optimal configuration.

This clearly demonstrates that a good accuracy over the parameter space of the EPM w.r.t. the true performances, such as the one being optimized by the learning algorithm (Random Forests here), is not required to reach the ultimate goal of the Algorithm Configuration process -find a quasi-optimal configuration for unknown instances. The only important property of the EPM is to be able to robustly identify good-performing regions of the parameter space. This opens several new possible research paths. At the level of the learning algorithm, the best regions of the parameter space could be weighted more than other parts of the space; at the extreme, rank-based learning could be used rather than regression of ERT values. At the level of the sampling, only good configurations could be used -e.g., the configurations encountered while running SMAC to find the true optimal configuration.

No clear differences can be seen between EPM learned using ψ • 500 or ψ • 2000 , except for some functions, in dimension 10, where EPMs learned with ψ • 2000 solve the function while those learned from ψ • 500 don't. On the other hand, using the full set of features ψ * does help, both in dimension 5 for the dimension-specific models, and in dimension 10 where it succeeds in reaching the target precision where the other models fail (e.g., F3 and F4, the only multi-modal functions of the test bench). While not surprising, this demonstrates that both the training set and the set of features should cover all the foreseeable difficulties of the unknown forthcoming instances. Any limited test bench (including BBOB) might hence be insufficient to learn a general-purpose configurator.

Finally, the fact that learning the EPM for a specific dimension leads to better results was to be expected. While this makes difficult to build an universal EPM, it does not prevent from any practical use of this approach, as the dimension is usually known (and constant) in most real world applications.

Conclusion

This paper has investigated the computation and use of an Empirical Performance Model (EPM) in the context of continuous black box optimization and has demonstrated that it is possible to learn a reasonable approximation of the real performance. More importantly, it was demonstrated that an efficient parameter configuration can be extracted from the learned EPM by optimizing the predicted performance, given a set of features on a new unknown function. In particular, it was possible to obtain empirical configurations that outperform a static parameter setting optimized for an average performance, over the whole test bench at the same overall cost. However, some open issues remain related to the robustness of the results, and deeper analyses are necessary to better understand (and avoid) some rare cases where the approach fails.

Several paths for further research are suggested by this work, both at the level of the learning algorithm and of the sampling of the parameter search space, as discussed in Section 6.3. A promising direction is to embed the EPM as a parameter control mechanism within the optimization process itself, assuming that the features can be efficiently approximated using a rather small number of samples, (e.g., w.r.t. an approximation of the objective function, as proposed in [START_REF] Belkhir | Surrogate Assisted Feature Computation for Continuous Problems[END_REF]). This would open a new perspective on the on-line parameter tuning grail.

  is the size of the initial sample, and only features not requiring any additional function evaluation beyond those of the initial sample are considered (i.e., Meta-Model, Information Content, Levelset, y-Distribution and Dispersion). Results with k = 500 or k = 2000 are presented in the following.

10 Fig. 1 :

 101 Fig. 1: Examples of comparisons between the true ERT (top) and the EPM (bottom) for 4 different functions and dimensions. Both EPMs of F1 (a) and F4 (b) have been learned only on samples from dimension 5, with features ψ • k=2000 while those of F11 (c) and F5 (d) have been learned on samples of all dimensions, and with ψ • k=500 . Each subplot shows performances colormaps (without interpolation) of log 10 (ERT d), for one DE strategy, with the 2 other DE parameters F and CR on the axes. The true optimal configurations are plotted as white stars and the Empirical Optimal Configurations as white small circles ○.

Table 1 :

 1 Percentages of the ERTs of different Empirical Optimal Configurations w.r.t. that of the robust configuration defined in Section 5, for dimensions 5 and 10. See text for details. Best results are printed in blue bold face when smaller than 100, in red italic when larger than 100; Worst results are printed in light gray; the ⋆ symbol indicates that the robust configuration never reached the target, and was artificially attributed an ERT of 15 times the maximum budget of one run, and ∞ indicates that the corresponding Empirical Optimal Configuration never reached the target.

		157 12.9 33 29 46 56 82	56
	F7	99 0.22 24 35 34 35 35	∞
	F8 114 17 ∞ ∞ ∞ ∞ ∞	29
	F9 107 15.4 ∞ 28 29 27 29	17
	F10 107 12.9 17 17 20 17 28 226
	F11 114 9.7 15 16 15 15 16	16
	F12 429 7.8 ∞ 21 10 ∞ 10	63
	F13⋆ 10 0.1 ∞ 2 ∞ ∞ ∞	1
	F14 100 14.4 ∞ ∞ 16 ∞ ∞	28
		(b) Dimension 10

The same θj need not have been tried for all fi.

d, the dimension of the search space, can be considered as the only external feature -or the Algorithm Configuration can be conducted anew for each dimension (more in Section 5).

http://www1.icsi.berkeley.edu/ ˜storn/code.html

http://coco.gforge.inria.fr

Measured as the number of function evaluations.

http://github.com/flacco

http://scikit-learn.org/

Additional plots are available at https://drive.google.com/open?id= 0B9GuQcCjvwtFdkotR1h1N3dlOG8


Nacim Belkhir

Johann Dréo

Pierre Savéant
email: pierre.saveant@thalesgroup.commarc.schoenauer@inria.fr

Marc Schoenauer

Surrogate Assisted Feature Computation for Continuous Problems

Keywords: Empirical Study, Black-box Continuous Optimization, Surrogate Modelling, Problem Features

A possible approach to Algorithm Selection and Configuration for continuous black box optimization problems relies on problem features, computed from a set of evaluated sample points. However, the computation of the features proposed in the literature require a rather large number of such sample points, unlikely to be practical for expensive real-world problems. On the other hand, surrogate models have been proposed to tackle the optimization of expensive objective function. It is proposed in this paper to use surrogate models to approximate the values of the features at reasonable computational cost. Two experimental studies are conducted, using the well-known BBOB framework as testbench. First, the effect of sub-sampling is analyzed. Then, a methodology to compute approximate values for the features using a surrogate model is proposed, and validated from the point of view of retrieving BBOB classes. It is shown that when only small computational budgets are available, using surrogate models as proxies to compute the features can be beneficial.

Introduction

Different optimization algorithms, or, equivalently, different parameterizations of the same algorithm, will in general perform non-uniformly on different classes of optimization problems. Indeed, it is today widely acknowledged in the domain of optimization at large that the quest for the general optimization algorithm, that would solve all problems at best, is vain, as proved by different works around the No Free Lunch theorem [START_REF] Auger | Continuous Lunches are Free plus the Design of Optimal Optimization Algorithms[END_REF][START_REF] Wolpert | No free lunch theorems for optimization[END_REF]. Hence, tackling an unknown optimization problem amounts to choosing the best algorithm among a given set of possible algorithms (algorithm selection), and/or the best parameter set for a given algorithm (algorithm configuration).

Such a choice can be considered as an optimization problem in itself, thus pertaining to the Programming by Optimization (PbO) paradigm proposed by [START_REF] Hoos | Programming by optimization[END_REF]:

given a new optimization problem instance (an objective function to minimize or maximize), a set of algorithms with domains for their parameters, and a performance measure, find the best1 algorithm or parameter setting to solve this problem. However, such a meta-optimization problem is in general difficult to solve (hierarchical search space, multi-modal landscape, . . . ) and thus requires running different algorithms with different parameterizations, where each of these runs will in turn call the objective function a large number of times: in total, the number of calls to the objective function will be huge, making such an approach intractable in most real-world situations with expensive objective functions.

The PbO approach can however be applied to classes of objective functions: the best algorithm/parameter setting can be learned off-line, once and for all for a given class of functions, and the optimal setting (algorithm + parameters) applied to all members of that class. This is the case for instance in operational contexts, where the same type of problem, but with slightly different settings, has to be solved again and again. However, in the general case of black-box optimization problems, very little domain knowledge is known (the type of search space for sure, maybe some dimensionning parameters) and such characteristics are not sufficient to reliably choose an algorithm and its parameters.

Another approach is then to compute some characteristics of objective functions, aka features, without much domain knowledge, and, thanks to a large example base of algorithms performances on known objective functions, to learn a performance model of several algorithms and their parameters. A well-known success in that direction has been obtained in the SAT domain [START_REF] Xu | Satzilla: portfolio-based algorithm selection for sat[END_REF], but dozens of features had been proposed in the literature to describe SAT problems and try to understand what makes a SAT problem hard for this or that algorithm. This situation is quite unique, and is an appeal for research regarding the design and study of features in other domains.

In particular, in the domain of continuous optimization (the search space is a subspace of IR d , for some d), several recent works [START_REF] Lunacek | The dispersion metric and the cma evolution strategy[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF] have proposed many different features to try to understand the landscape of continuous optimization and, ultimately, solve the Algorithm Selection and/or Configuration problem. Note that a large body of mathematical programming algorithms exist, and are proved to be optimal for specific classes of objective functions: Linear Programming should be applied if the function (and the constraints) are linear; gradient-based algorithms should be applied if the function is convex, twice differentiable and well conditioned, etc. But these are rare exceptions in the real world, where the general problem remains open, and the feature-based approach seems worth investigating, in particular considering the promising initial results obtained by [START_REF] Mersmann | Exploratory landscape analysis[END_REF] (more in Section 2).

Unfortunately, the computation of all proposed features is based on many sample points, i.e., values of the objective function at given (generally random) points of the search space. In real-world situations, where the objective function is expensive and the computational budget limited, the computation of the features as proposed in [START_REF] Mersmann | Exploratory landscape analysis[END_REF] might simply be impossible. A first research question is hence to study how badly the features behave when the size of the sample used to compute them decreases.

Yet, a prominent approach has already been proposed to handle expensive objective functions in optimization. It relies on surrogate models, i.e., models of regression or interpolation of the objective function built upon sample points gathered during the run of the algorithm, and used from time to time in the optimization algorithm in lieu of the actual objective function.

Building on this idea, the present work investigates the use of surrogate models to compute problem features: based on a small sample set, a surrogate model is developed. The features of the surrogate model can then be easily computed, since the cost of evaluating a surrogate model is negligible compared to the original objective function. The second issue investigated in this paper relates to the accuracy of the features of the surrogate model as approximations of the features of the original objective function, and will be studied here experimentally, using the well-known BBOB test set of functions [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF]. Note that this paper will not directly tackle the algorithm selection or algorithm configuration problem, left for further work. The accuracy of feature sets will be assessed first by a direct comparison with the features computed using a very large sample set, as well as by their abilities to correctly recover the BBOB (hand-designed) classes, as in [START_REF] Mersmann | Exploratory landscape analysis[END_REF] (see Section 4 for more details on the BBOB testbench).

The paper is organized as follow, Section 2 surveys the features proposed in the literature. Section 3 will introduce the methodology proposed in this paper. Section 4 describes the experimental context, while Section 5 describes the experimental protocol in detail, including the performance measures used to validate the results presented in Section 6. Finally these results are summarized and hints for further works are given in Section 7.

Problem Features for Continuous Optimization

Let F be the objective function at hand, defined on a domain D ⊂ IR d for some given dimensions d: d is the only high-level feature given a priori as domain knowledge. All (other) features considered here will be computed from a sample set X = (x 1 , . . . x p ) of points of D, and the corresponding values of the objective function Y = (y 1 , . . . , y p ) (with y i = F(x i ) for all i).

A first series of features used here is taken from [START_REF] Mersmann | Exploratory landscape analysis[END_REF], where low-level features are computed directly from (X , Y). High-level features,or fitness landscape properties(e.g. multi-modality,separability, convexity,etc), are then built from different statistics on these low-level features. A total of 62 low-level features is computed, giving in turn 8 high-level features. The low-level features are grouped in the following classes:

-Distribution: these features consider the distribution of objective values in Y, computing the skewness and the kurtosis of the distribution, but also estimating the number of peaks of the distribution. -Level Set: from a given threshold in the objective function values (e.g., some quantiles), a classification technique such as LDA, QDA or MDA is trained to predict the position w.r.t. the threshold. The distribution of the misclassification errors for the different classifiers are used as proxies to the multi-modality of the objective function. -Meta-Model : First, linear and quadratic regression models for (X , Y) are computed. The resulting R 2 coefficients for the accuracy of the models, as well as some statistics on the relative sizes of the coefficients of the models give some indication on the shape of the landscape. -Convexity: estimate the probability of convexity and linearity of the objective function by selecting two random points from the sample X and generating a new sample point in the segment and comparing the objective value of the new sampled point and the same convex combination of the two initial objective values. The probability of convexity is computed by averaging the number of trials where the computed difference is lower than a predefined negative threshold, while the probability of linearity is computed by considering all absolute differences that are smaller than the absolute value of the predefined threshold. -Curvature: considers the numerical approximation of the gradient in each point of a sub-sample of X by the Richardson's extrapolation method, and the resulting features consider the basic statistics of the respective derivatives, the condition number of the similar numerical approximation of the Hessian. -Local Search: a local search algorithm, e.g., Nelder-Mead, is run from starting points in a sub-sample of X ; the final solutions of all runs are clustered in order to identify the local optima of the objective function, the basin sizes are approximated by the number of local searches which terminate there, and other statistics gathered during the different runs give other indicators of the landscape properties.

Furthermore, a series of features termed Dispersion was originally proposed in [START_REF] Lunacek | The dispersion metric and the cma evolution strategy[END_REF]. Their computation analyzes the distance between candidate solutions for a percentile of the best solutions of the optimization problem by comparing them to the mean or median distance between solutions in the whole initial sample. Different percentiles are considered, giving in total 16 features.

Finally, features related to Information Content have been proposed in [START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF]. They are related to the number of the binary decisions needed to find the information, such that each candidate solution is binarized with respect to the fitness value of their nearest neighbors fitness. These metric giving information about the smoothness or a global structure of the objective function. A total of 5 dispersion features are considered in [START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF].

The recent works that defined those features [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF] successfully demonstrated that these features could be used in order to classify the optimization problems with respect to their classes in BBOB (see definition in Section 4). However, some of these features (local search, curvature, convexity) require additional samples. On the other hand, in [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF], only features that can be estimated on a fixed initial sample are considered, and the results demonstrate that such limited number of features (33 features) can nevertheless be used to correctly identify the BBOB classes. It should be emphasized, however, that all the abovementioned works consider large samples of candidate solution (between 500 × d and 5000 × d), which is not practical if the objective function is expensive.

3 Surrogate Assisted Feature Computation

Sub-sampling

Real world applications where 3D numerical simulations are involved, result in very expensive objective functions, for which the computation of a single value might take a few hours. For such functions, whereas the features described in previous Section 2 might help choosing the right algorithm, and/or the right parameters of a given algorithm, thus saving computation time for its optimization, the computation of these features should not cost more than its optimizationwhich might be the case if around 1000 × d evaluations are needed, as in [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF].

A first solution is of course to simply use a small sample set. The price to pay will be a poor approximation of the actual features, which might result in a wrong choice for the optimization algorithm or its parameters: as can be expected (see Section 6.1), a too small sample set results in a poor approximation with a high variance w.r.t. the choice of the sample set.

Surrogate Modeling

Coming from another field, numerical engineers have tackled this problem using Response Surface Methods for many decades now: after few iterations of any optimization algorithm, many points of the search space have already been evaluated, and this sample set can be used to build a surrogate model of the objective function, that can in turn be used in the optimization algorithm as a proxy for the actual objective, being costless to compute (see e.g., [START_REF] Jin | Comparative studies of metamodeling techniques under multiple modeling criteria[END_REF][START_REF] Jin | A Comprehensive Survey of Fitness Approximation in Evolutionary Computation[END_REF] for surveys in the engineering domain and in in Evolutionary Computation domain, respectively). The most critical issue to be addressed when using surrogate modelling techniques is the choice of the function space where to look for a model. Several approaches have been proposed in the literature to solve such regression problem, from Neural Networks to Gaussian Processes (aka Kriging) to Support Vector Machines and Regression Random Forests, and it is beyond the scope of this paper to describe them in detail.

However, whereas a surrogate model to be used within an optimization process should be accurate in the region that is currently searched (generally around the current best solution found), because it will be updated as the search progresses, or a surrogate model whose purpose is to compute features of the function at hand should be globally approximating the objective function: both goals are different, and this should be taken into account.

Accuracy vs Efficiency

The baseline of the approach proposed in this paper is to use small sample sets to build a surrogate model of the objective function at hand, to compute the features of that surrogate model using as many samples as needed, and to use these features in lieu of the unreachable actual features of the true objective function.

However, this approach must be validated -experimentally -against the simple sub-sampling approach (compute the features on small sample sets). Furthermore, some parameters of the approach (the model for the surrogate and its hyper-parameters, the number of samples to be used for evaluating the features of the surrogate model) must also be tuned. This raises the question of what measure should be used to assess the quality of the approximated features.

An obvious measure is simply the error made on the feature values: for each feature, the "exact" value can be computed using as large sample sets as needed (say, same sizes than in [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF]), the L2 norm of the error vector, difference between the approximated end "exact" values for all features gives a first idea of how good the approximation is. However, the ultimate reason for computing the features is to solve the algorithm selection and/or configuration problem. And it might be the case that the error in solving the latter problem varies differently from the L 2 error on the feature values (e.g., some features that are poorly approximated are also not so important for the algorithm selection and/or configuration). Because there is not yet any standard algorithm selection or classification problem that would allow us to do such a validation, some obvious proxy classification problems that are used in [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF], retrieving the known classes of problems manually defined on BBOB testbench (described in next Section). Next Section will give the technical details of these experiments, and their results will be presented in the following Section.

Experimental Settings

BBOB testbench

All the experiments presented here use test functions from the Black Box Optimization Benchmark (BBOB) 2 [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF]. The BBOB benchmark contains 24 analytically defined continuous objective functions with known global optimum. All these functions are defined on the d-dimensional domain [-5, 5] d , with d ∈ {2, 3, 5, 10, 20, 40}. In order to avoid possible biases, 15 variants of each function are considered, obtained from the original function by a translation of the position of the optimum, plus a rotation of the coordinate system for the nonseparable functions. Noisy variants are also available in BBOB, but only the noiseless functions will be considered here.

These 24 functions have been manually classified in five classes of problem, with 5 separable functions, 4 uni-modal functions with low or moderate conditioning, 5 uni-modal functions with high conditioning, 5 multi-modal functions with regular global structure, 5 multi-modal functions with weak global structure.

Sample Sets

All sample sets are drawn uniformly on [-5, 5] d for a BBOB function in dimension d (see above). In the remaining of the paper, all sample set sizes are normalized w.r.t. the dimension d of the definition domain of the objective function. For the sake of brevity, we will only mention the ratio between the sample set size and the dimension: "a sample of size k" will actually mean "a sample of size k × d". In all experiments, k ∈ {30, 50, 100, 500, 1000, 2000} (the largest value 2000 was decided after the first experiments, see Section 6.1).

Features

As discussed in Section 2, only part of the features described in [START_REF] Lunacek | The dispersion metric and the cma evolution strategy[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF] will be used here, namely: 3 Distribution features, 9 Meta-Model features, 16 Dispersion features, and 5 Information Contents features. All considered features are implemented as an R package publicly available at http://github.com/flacco, thanks to Pascal Kerschke.

Surrogate Models

As discussed in Section 3.2, several approaches to surrogate modelling will be used and compared: Gaussian Processes, Random Forests, Support Vector Machines with polynomial or RBF kernel, denoted respectively GP, RF, SVM P and SVM RBF . But learning a surrogate model requires some hyper-parameters to be tuned: a grid search is performed in the hyper-parameter space (4 parameters for GP, 5 for the other models), using a 5-folds cross-validation procedurerandomly re-sample (Bootstrap) the sample set, using 80% for training the surrogate and the remaining 20% for testing-for 300 iterations, optimizing the approximation accuracy on the test set. All surrogate modelling procedures are implemented w.r.t the python scikit-learn library3 .

Experimental Protocol and Validation

Experimental protocol and notations For a given objective function F (one trial of one instance of one d-dimensional function from BBOB), the basic experiment goes as follows: one sample sets of given size is drawn from the definition domain of F; the features are computed on this sample set, and a surrogate model using one of the chosen modelling techniques. The sample set is then completed with more samples, using the surrogate model in lieu of the original function. Approximate features are then computed using this extended sample set.

An immediate validation of such approximated feature values can be made by comparing them to the "exact" values: a proxy for these values will be the features computed with the largest initial sample set, of size 2000 (see Section 6.1). However, the global validation (see below) requires to compute such approximated features for all BBOB functions, and for several different sample set sizes.

Let X s , s = 1, . . . , S be some sample sets from the domain of F4 . Features Φ(F s ) (vector of IR F if F is the number of features) are computed for F from sample set (X s , Y s ) (with y i = F(x i ) for all (x i , y i ) ∈ (X s , Y s ), see Section 2). Let us denote by Φ(F * ) the features computed from the largest sample set (of size 2000).

Each sample set X s is also used to learn some surrogate models F s t using different surrogate modelling techniques t = T 1 , . . . , T T . For each (s, t), the set of features Φ( F ) can then be compared to the "exact" features Φ(F * ), and their accuracy assessed using the L 2 norm in IR F of the errors on the feature values (Err(

F s ) = ||Φ(F s )-Φ(F * )|| 2 , Err( F s t ) = ||Φ( F s,s t ) -Φ(F * )|| 2 )
. However, as discussed in Section 3.3, another comparison is needed between the approximate features and the "exact" values, that relates to the ability of the approximate features to be sufficient to correctly classify the BBOB classes. Such validation requires the computation of all approximate features with same sizes of sample sets for all instances of functions of BBOB testbench. Classification Efficiency However, measuring the efficiency of a set of approximated features as a whole goes through using them as input for learning a classifier in order to discriminate the five BBOB classes. This is done using a 5-folds cross-validation, repeated 100 times, procedure as follows. Let us denote Cl(F) the class (in 1..5) a given function F belongs to. For a given sample set size s, the example set for the classification task consists of (Φ(F s ), Cl(F)) pairs when dealing with features computed on F and (Φ( F s,s t ), Cl(F)) when dealing with surrogate model t built on F (t ∈ {GP, RF, SV M P , SV M RBF }). Such example set is made of 5 trials × 5 instances × 24 functions × 5 dimensions. Out of these 3000 examples, 80% are randomly chosen without replacement, equally distributed in the 5 classes, to build the training set, on which a Random Forest classifier is trained (with default hyper-parameters from scikit-learn).

The accuracy of the resulting classifier should then be assessed on the remaining 20% of the global example set. However, different scenarii are possible in real-world situations. A first scenario is when the training phase is done on easy functions, for which it is possible to compute the features with large enough sample sets, and the unknown functions on which to perform algorithm selec-tion/configuration (the test phase) are all expensive. An "orthogonal" scenario is when the functions available for training are also expensive. In the latter case, only approximate features will be available for training, either computed on small samples, or computed using a surrogate of the functions used for training.

Two situations similar to the ones described above will be experimented with here, involving different sample sizes for the training of the classifier and its test.

When no surrogate model is involved (Section 6.1), an experiment studying the efficiency of the approximate features as a basis for classification (Section 6.1) is defined with only two parameters: the size s train of the sample set used to learn the features for the training set, and the size s test of the sample set used to learn the features for the test set. The classification accuracy of the resulting classifier will be denoted Ef f (s train , s test ).

But when analyzing the efficiency of surrogate assisted feature computation (Section 6.2), an experiment is defined with 5 parameters: the type T of surrogate model (in {GP, RF, SV M P , SV M RBF }), and, for both the training features and the test features, the sizes of the original sample sets used to learn the surrogate models (respectively s org train and s org test ), and the additional number of points added to these original sample sets using the surrogate models (respectively s surr train and s surr test ). The classification accuracy of the resulting classifier will be denoted Ef f (T, s org train , s

Results

Two series of experiments will be presented here. The first one (Section 6.1) doesn't involve any surrogate model, and aims at studying how the features diverge from their "exact" baseline values when the size of the sample decreases. The goal of the second series (Section 6.2) is to check whether using a surrogate model built on the same available small sample set to complement it can help to cope with such divergence. In both series, the divergence with the baseline values will be assessed by the accuracy of the approximated values (individual comparison for a given feature and a given function, and their aggregation in the L 2 error -Sections 6.1 and 6.2 respectively), and by the efficiency of the whole set of approximate features, using them to discriminate the 5 BBOB classes, as explained in previous Section -Sections 6.1 and 6.2 respectively).

Effects of Sub-sampling

Accuracy of Sub-Sampled Features Five test functions (F1, F8, F13, F17, F23, one per BBOB class) are used here to assess the effect of sub-sampling on the feature values. Figure 1 shows some typical feature behaviors on those 5 functions: whereas plot (a) display a smooth behavior, where feature values stabilize for s ≥ 100, both other plots show that even s = 2000 might still be somehow too small for the multi-modal functions F17 and F23. However, most features on most functions exhibited a smooth behavior, and were stable for s ≥ 500, justifying the decision to take Φ(F 2000 ) as Φ(F * ). It is also clear from Figure 1 that small sample sizes (e.g., 30) will provide a poor approximation of the feature values, and might not allow to discriminate among different functions.

Efficiency of Sub-Sampled Features Figure 2 displays the efficiency of the approximated features to discriminate among BBOB classes (see Section 5). Each line corresponds to a sample size s train used to train the classifier, and each point corresponds to a different sample size s test used to compute the features of the test instance to be classified. Two conclusions can be drawn from this figure. First, there is no reason to use a larger test sample size s test than the size s train that was used to train the classifier, as the efficiency does not increase after s test has reached s train . Second, if you know that only a limited budget will be available at test time (i.e., all new instances will be very expensive), then you should train the classifier with a small budget too: for a given s test , the best efficiency is obtained by the classifier trained with s train = s test .

A possible explanation, to be investigated deeper in further work, is that sub-sampling does not only increase the variance of the feature values, but it also induces some bias that might be also tracked by using the same sample size for training than for testing.

Surrogate Assisted Feature Approximation

In this Section, we try to improve the accuracy and efficiency of the approximate features by adding samples computed using some surrogate model, as described in Section 3. function alone as displayed in Figure 1. However, it should be noted that the different surrogate models can give different behaviors on the same feature. An alternative point of view and a comparison with the approximated features directly computed using the initial small sample with exact objective values is given in Figure 5. Here, the Random Forest surrogate model is used to add 2000 points to the initial sample set. It is clear (and results on other feature confirm this trend) that Random Forests give a much smaller error than Support Vector Machines (with both polynomial and RBF kernel), and even more so with Gaussian Processes (not shown here for space reason). More interestingly, using the Random Forest surrogate model often result in more accurate approximate features than computing their values only on the few available exact values (see the s = 30 histograms on Figure 5).

Accuracy of Surrogate-Assisted Features

But whereas adding as many points as possible (2000 here) using the surrogate model seems a natural way to go, one can wonder if this is always the best thing to do. And the answer is no: Figure 6 For each (s, s ), the error Err( F s,s RF ) is plotted (the darker the higher). obtained when adding s (y-axis) sample points using a surrogate model (here, a Random Forest), to a sample set of size s (x-axis). In many cases, the error increases with s , or displays an unstructured behavior, in particular for small values of s. Only for large-enough values of s (500 and more, not shown here) does the error decrease.

Efficiency of Surrogate-Assisted Features Let us now look at the efficiency of the approximated features to discriminate among BBOB classes (as explained in Section 5). Figure 3 displays Ef f (s, s) (the upper hull of the plots on adjacent Figure 2) as well as the different Ef f (T, s, * ; s, * ), for T ∈ [GP, SV M P andSV M RBF , RF ]. It is again obvious here that the Random Forest model outperforms all others, which is consistent with previous results (as well as with all other not presented results). From now on, only RF surrogate model will be considered. But another observation that can be made here is that there is no gain to be expected by using the surrogate model during the learning phase too, as both plots for Ef f (s, s) and Ef f (RF, s, * ; s, * ) are almost identical.

A final experiment will try to answer the main question that motivated this work: can the use of a surrogate model improve the efficiency of the approximated features in the context of expensive objective functions? Figure 7 And indeed, there is some advantage in using the surrogate model during the test phase, the more so for small training budgets. Furthermore, this advantage of using the surrogate is statistically significant, as witnessed in Table 1 where the same data are given together with the standard deviations. Figures in bold are statistically better than the corresponding non-bold figures according to a Wilcoxon signed rank test with 95% confidence.

Discussion and Conclusion

In the context of continuous black-box optimization, this paper has proposed a methodology to compute features describing the characteristics of the problem at hand, as proposed in [START_REF] Bischl | Algorithm selection based on exploratory landscape analysis and cost-sensitive learning[END_REF][START_REF] Mersmann | Exploratory landscape analysis[END_REF][START_REF] Munoz | Exploratory landscape analysis of continuous space optimization problems using information content[END_REF], using surrogate models to cope with expensive objective functions: to-date methods to compute such features rely on large sample sets of evaluated points, that are not practically available when dealing with expensive real-world problems. The performance of approximated features have been measured both with their accuracy, related to the error on their values when compared to values computed on very large sample sets, and with their efficiency, ability to train a classifier that can correctly discriminate the five classes of the test functions in the BBOB testbench [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF].

The paper has first experimentally studied the loss of accuracy of the features due to sub-sampling, identifying a reasonable sample size beyond which the computed features exhibit stable values with small variance -2000× the problem dimension. The study of the efficiency of the sub-sampled features led to a first conclusion: if the budget at test time is going to be small (expensive objective function), then the budget allocated to training should be small too.

Experiment involving surrogate models first surprisingly demonstrated that, in the context of the present work, only Random Forest surrogate models gave satisfactory results when used to augment the sample set used to compute the features, whereas for instance Gaussian Processes, very popular today when it comes to expensive optimization per se, gave the worst results of all.

But the most interesting observation is that when only few samples are available, using a surrogate model built on these small sample sets to augment the sample set can lead to better classification results (on BBOB classes) than the sub-sampled features directly computed on the small available sample sets. However, if some large training budgets are available, i.e., if there exists cheap functions that are representative enough of the expensive real-world objective functions that will be encountered later, then using a surrogate model on these expensive unknown function does not seem to be beneficial.

There are several directions for further work. First, several features that have been proposed in the literature have been discarded because they require additional function evaluations (e.g., the local search, curvature, and convexity features [START_REF] Mersmann | Exploratory landscape analysis[END_REF]). Surrogate models might also be useful in order to compute those features at low cost. Such progress would then make even more important to use some feature selection method rather than using all features for the classification task at hand.

Different paths can also be explored regarding the way the surrogate models are computed. First, as of today, only global surrogate models are considered. But it might be interesting to restrict the scope of the surrogate learning to better compute local search features for instance. Similarly, the performance measure used to construct the surrogate models is the standard approximation error on the known samples. But in the context of feature computation, some other measures could be considered -ultimately, surrogate models should be optimized for the quality of the approximate features they allow to compute.

Only global results on all features have been presented here. In particular, it could be the case that different surrogate models are the best choice to compute different features (e.g., SVM for convexity features, random forests for multimodality features, etc). A deeper insight on feature accuracy is needed here. Along the same line, some results presented here suggest that there might be a systematic bias induced by sub-sampling -something that requires deeper analysis too.

Last but not least, the link between those features and Algorithm selection and configuration remains to be established. Indeed, retrieving BBOB classes is a much easier problem than learning an empirical performance model, as in [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF] for instance. Or those features might be used together with some latent features identification method, as in [START_REF] Mısır | Algorithm selection as a collaborative filtering problem[END_REF].

Fig. 1 :

 1 Fig. 1: Examples of effect of sub-sampling (x-axis) on feature values for the different test functions in dimension d = 5.

Figure 4 Fig. 2 :Fig. 3 :

 423 Fig. 2: Ef f (s train , s test ) vs s test for different values of k train . Fig. 3: Ef f (s, s) (black circles) and Ef f (T, s, * ; s, * ) for different T .

Fig. 4 :

 4 Fig. 4: Values of 3 features vs initial sample size s, for the 5 test functions (d = 20), using a surrogate model (GP for (a,b,c), RF for (d,e,f) ) to add s = 2000 points to the sample set.

Fig. 5 :

 5 Fig. 5: Accuracy error for 3 features on F13 (d = 5), for different values of s (x-axis). For each s, Err(s, s), Err( F s,2000 SV M P ) and Err( F s,2000 RF ) are plotted.

Fig. 6 :

 6 Fig. 6: Accuracy error of 4 features on F1, d=2 (left), and on F17, d=20 (right).

Table 1 :

 1 displays 4 plots, Mean and Standard deviation of the efficiency of sub-sampled features Ef f (s train , s test (columns Obj. Fn.) and surrogate assisted features Ef f (RF, s train , 0; s test , * ) (columns RF). Statistically sig nificantly better results (Wilcoxon signed test with 95% confidence) are in bold.corresponding to different values of s train . On each plot, the black continuous line is Ef f (s train , s test ), i.e., the corresponding line of Figure2, and the dotted grey line shows Ef f (RF, s train , 0; s test , * ), i.e., the efficiency obtained when using the RF surrogate model to augment the sample set with 2000 new samples.

	100			s train = 30				100		s train = 100
		90									90
	Efficiency	70 80								Efficiency	70 80
		60						Obj. fn. RF		60						Obj. fn. RF
		50	30	50	100	s test	500	1000	2000		50	30	50	100	s test	500	1000	2000
		100		s train = 500				100		s train = 2000
		90									90
	Efficiency	70 80								Efficiency	70 80
		60						Obj. fn. RF		60						Obj. fn. RF
		50	30	50	100	s test	500	1000	2000		50	30	50	100	s test	500	1000	2000
			strain = 30				strain = 100				strain = 500				strain = 2000
	stest Obj. fn.	RF		stest Obj. fn.	RF	stest Obj. fn.	RF		stest Obj. fn.	RF
	30 85.1±1.6 87.4±1.5	30 66.9±3.7 82.2±1.5	30 59.5±3.7 65.2±3.1	30 54.8±3.5 60.9±3.1
	50 87.3±1.4 90.0±1.3	50 82.6±3.0 89.3±1.2	50 73.8±2.1 69.8±3.1	50 66.6±3.0 67.1±3.3
	100 86.5±1.5 90.2±1.2	100 94.6±0.9 95.1±1.0	100 88.3±2.0 76.3±2.8	100 78.6±3.6 73.1±3.0
	500 87.0±1.6 85.6±2.0	500 95.2±1.1 94.3±1.0	500 98.6±0.5 87.6±1.8	500 95.3±1.1 84.0±2.4

Fig. 7: Comparison, for different values of s train , of the efficiency of sub-sampled features Err(s train , s test (continuous black line) with that of surrogate assisted features Ef f (RF, s train , 0; s test , * ) (grey dotted line).

The performance measure generally involves time-to-solution (CPU time, or number of function evaluations) and precision/accuracy of the solution returned by the algorithm (precision of the solution for continuous optimization problems, number of constraints violated for Constraint Programming problems, etc).

http://coco.gforge.inria.fr

http://scikit-learn.org/

By abuse of notation, s will denote both the sample set and the (normalized) size of the sample set.