Borja Calvo
email: borja.calvo@ehu.es

Ofer M Shir
email: ofersh@telhai.ac.il

Josu Ceberio
email: josu.ceberio@ehu.es

Carola Doerr
email: carola.doerr@mpi-inf.mpg.de

Hao Wang
email: h.wang@liacs.leidenuniv.nl

Thomas Bäck

Jose A Lozano
email: ja.lozano@ehu.eus

Bayesian Performance Analysis for Black-Box Optimization Benchmarking

Keywords: Computing methodologies → Randomized search, • Mathematics of computing → Statistical paradigms, performance measures, Bayesian inference, Plackett-Luce model, black-box optimization, benchmarking, evolutionary algorithms

The most commonly used statistics in Evolutionary Computation (EC) are of the Wilcoxon-Mann-Whitney-test type, in its either paired or non-paired version. However, using such statistics for drawing performance comparisons has several known drawbacks. At the same time, Bayesian inference for performance analysis is an emerging statistical tool, which has the potential to become a promising complement to the statistical perspectives offered by the aforementioned p-value type test. This work exhibits the practical use of Bayesian inference in a typical EC setting, where several algorithms are to be compared with respect to various performance indicators. Explicitly, we examine performance data of 11 evolutionary algorithms (EAs) over a set of 23 discrete optimization problems in several dimensions. Using this data, and following a brief introduction to the relevant Bayesian inference practice, we demonstrate how to draw the algorithms' probabilities of winning. Apart from fixed-target and fixed-budget results for the individual problems, we also provide an illustrative example per groups of problems. We elaborate on the computational steps, explain the associated uncertainties, and articulate considerations such as the prior distribution and the sample sizing. We also present as a reference the classical p-value tests.

INTRODUCTION

Benchmarking environments and competitions are designed to compare algorithms' performance over a set of predefined representative problems. In Evolutionary Computation (EC), benchmarking black-box optimization problems evolved differently in continuous versus discrete domains. While the Black-box Optimization Benchmarking suite (BBOB, [START_REF] Hansen | Comparing Results of 31 Algorithms from the Black-box Optimization Benchmarking BBOB-2009[END_REF]) constitutes an established testing framework for evaluating performance of continuous optimizers, the discrete domain, on the other hand, has not had the benefit of an equivalent suite. Attempts to establish such an environment have lately become prominent [START_REF] Doerr | IOHprofiler: A Benchmarking and Profiling Tool for Iterative Optimization Heuristics[END_REF][START_REF] Rapin | Nevergrad -A gradient-free optimization platform[END_REF][START_REF] Wasik | Optil.io: Cloud Based Platform For Solving Optimization Problems Using Crowdsourcing Approach[END_REF][START_REF] Weise | Optimization Benchmarking[END_REF]. Wherever performance comparisons are sought, based on empirical data, they call for statistical assessments, to evaluate whether the observed performance gaps can be supported by an appropriate estimator for the true, underlying performance distribution, i.e., a distribution which assigns a probability to each possible result of the algorithms. Within the typical EC context, arrays of recorded optimization runs merely constitute samples of the corresponding performance distribution.

Statistical tests fall into the category of statistical inference methods, which in this case seek to deduce measures about the unknown population 1 (the set of all possible results), when given only a sample drawn from that population. The most commonly used statistics in EC are of the Wilcoxon-Mann-Whitney-test type [START_REF] Wineberg | Introductory Statistics for EC: A Visual Approach[END_REF], which facilitates non-parametric pairwise performance comparisons, in either paired or non-paired scenarios [START_REF] Gopal | 100 Statistical Tests[END_REF]. In its non-paired version, the test concentrates on a null hypothesis, H 0 , stating that the compared two populations share the same mean. Furthermore, this test assumes that both populations share the same shape and spread of their underlying continuous frequency distributions [START_REF] Gopal | 100 Statistical Tests[END_REF]. To test this hypothesis, a sample from the populations, denoted as s and referring to two arrays of recorded runs of the algorithms, is used to get a certain statistic R(s), whose probability distribution under H 0 (Pr H 0 (R)) is known, given the set of assumptions. Then, the so-called p-value could be calculated, representing the probability of observing such an extreme value for the statistic, p = Pr H 0 (R ≥ R(s)), subject to all the assumptions (including H 0 ).

Small p-values, e.g., p ≤ 0.05, are taken as evidence of H 0 being false, that is, providing confidence that the averages of the populations significantly differ. In other words, such p-value are used as a statistical support to claims concerning, for example, outperformance. However, often times this probability measure is misinterpreted as the probability of H 0 being true (or, alternatively, 1 -p misinterpreted as the probability of the counterpart hypothesis to hold, i.e., the averages being different) [START_REF] Greenland | Statistical tests, P values, confidence intervals, and power: a guide to misinterpretations[END_REF].

Recently, Bayesian inference [START_REF] Dennis | Introduction to Probability and Statistics from a Bayesian Viewpoint: Inference[END_REF] has been proposed as an alternative to analyze algorithms' performance within the Machine Learning community [START_REF] Benavoli | A Bayesian Wilcoxon signed-rank test based on the Dirichlet process[END_REF][START_REF] Benavoli | Imprecise Dirichlet process with application to the hypothesis test on the probability that X ≤ Y[END_REF], and has since been gradually adopted therein [START_REF] Benavoli | Time for a Change: a Tutorial for Comparing Multiple Classifiers Through Bayesian Analysis[END_REF]. In this alternative, rather than testing a certain hypothesis, the methods focus on estimating relevant information about the underlying performance distribution, represented by a set of parameters θ (e.g., subscribing to a performance difference between algorithms, or the probability of this difference being in a certain interval). They directly assess the distribution of θ conditioned on a sample s drawn from the performance distribution, which formally reads Pr(θ |s).

The current study aims to demonstrate the effectiveness of Bayesian inference analysis for algorithms' performance ranking. To this end, following a short introduction to this analysis, we report on its application to a specific benchmarking campaign, which comprises a set of 23 discrete optimization problems. Technically, the actual benchmarking execution relies on the recently announced IOHprofiler [START_REF] Doerr | IOHprofiler: A Benchmarking and Profiling Tool for Iterative Optimization Heuristics[END_REF], and more concretely on the data of the experiments described in [START_REF]Anonymous for double-blind reasons[END_REF].

BAYESIAN INFERENCE FOR ALGORITHM RANKING ANALYSIS

Comparing heuristic optimization algorithms implies handling the natural uncertainty associated with empirical results. Usually, this uncertainty stems from two main sources: the stochastic nature of the algorithms, and the extremely large number of possible instances. Therefore, in order to draw sound conclusions, it is not only essential to craft the design of the experimentation very carefully. It is equally important to apply statistical tools to the analysis of the results. Statistical tools are applied to do inference, i.e., to extract information about a certain underlying probability distribution using a sample drawn from that distribution. If this distribution represents, for example, the probability associated with the possible results obtained by two stochastic algorithms when run on a certain problem instance, then a sample of such a distribution would be a certain number of random runs of the algorithms on that instance.

Based on that sample of results, statistical inference methods are used to obtain information about the performance of the algorithms on that particular instance. The most commonly used statistical tools for this purpose are statistical tests. These methods require the formulation of a null hypothesis and provide a measure on whether or not this hypothesis can be refuted with a certain probability. A typical example for a null hypothesis is the statement that the average performance of two or more algorithms is the same for a given set of problems. The statistical test uses the empirical evidence, i.e., the results obtained by the algorithms on one or more instances, to decide whether or not the null hypothesis can be rejected with an appropriate degree of confidence. Therefore, this type of statistical tool provides us with two possible answers. If the null hypothesis is rejected, then we state that "there are statistically significant differences between the algorithms". Otherwise, we cannot make a statement, since there is not enough statistical evidence to ensure that the average performance is different -which is not the same as stating that there are no differences.

Therefore, statistical tests can either confirm that there are differences or simply tell us that we have not enough data to confirm the differences. At this point it is important to note that, in almost all realistic use cases, we already know that there are differences between the algorithms (no matter how small), as the algorithms are of different design. In this sense, some authors argue that using statistical tests when comparing algorithms is not informative enough to draw sound conclusions from the experimental data [START_REF] Benavoli | Time for a Change: a Tutorial for Comparing Multiple Classifiers Through Bayesian Analysis[END_REF]. One of the weak points of statistical tests is that, when the null hypothesis is not rejected, it is not a trivial task to explain ite.g., either by the small sample size or by the similar performance of the algorithms. As an alternative to statistical tests, Bayesian methods provide us with interpretable information, distinguishing between the estimated measures and the uncertainty about the estimations [START_REF] Benavoli | Time for a Change: a Tutorial for Comparing Multiple Classifiers Through Bayesian Analysis[END_REF].

Rather than having a single probability distribution to model the data, Bayesian statistics simultaneously considers all possible distributions of a certain family and assigns a probability to each distribution. Therefore, we represent the probability of the data with a certain probability model whose parameters themselves have an associated probability distribution. Before observing any data, the probability of the parameters is represented by the prior distribution. This distribution represents our prior belief about which are the true parameters of the distribution from where the data comes. Obviously, when we have access to actual data, we have to update this belief accordingly. This is done using Bayes' rule. If we represent the data (our sample) as x and its probability (also known as likelihood function) as p(x|θ ) where θ represents the set of parameters of the distribution, then we have:

Pr(θ |x) = Pr(x|θ )Pr(θ ) Pr(x) , (1)
where p(θ ) denotes the prior distribution of the parameters, and p(θ |x) stands for the posterior probability of the parameters given the data. We can use this posterior probability to obtain relevant information about the distribution (population) for where the data comes, such as the (posterior) expected value for X . Moreover, we can also derive information about the uncertainty of the estimations. Several Bayesian methods have been proposed to compare algorithms' performance by analyzing the underlying empirical data [START_REF] Benavoli | A Bayesian Wilcoxon signed-rank test based on the Dirichlet process[END_REF][START_REF] Benavoli | Imprecise Dirichlet process with application to the hypothesis test on the probability that X ≤ Y[END_REF][START_REF] Cassio | Joint Analysis of Multiple Algorithms and Performance Measures[END_REF]. However, to the best of our knowledge, all of them focus on the pairwise comparisons of algorithms. Despite the high relevance of these methods, they come short when attempting to address the following, very commonly asked, broad comparison question:

Given a set of algorithms, and given their associated empirical data when run on a predefined set of problems: Which of those algorithms performs best for those problems?

One approach to address this question is to explicitly model the ranking of the algorithms. In [START_REF] Calvo | Bayesian Inference for Algorithm Ranking Analysis[END_REF] a first attempt to answer that question was made by proposing a Bayesian Plackett-Luce model to perform the comparative analysis of experimental results.

Bayesian Plackett-Luce model

The Plackett-Luce model is a probabilistic model defined over the space of permutations [START_REF] Luce | Individual Choice Behavior[END_REF][START_REF] Plackett | The Analysis of Permutations[END_REF], which in our context correspond to rankings of algorithms. Among the probabilistic models proposed for permutations (i.e., functions that assign a probability to each permutation) [START_REF] Irurozki | Sampling and Learning Distance-based Probability Models for Permutation Spaces[END_REF], the Plackett-Luce model collects some features that make it suitable to model rankings of algorithms. The first one is that it fulfills Luce's axiom, which states that the relative ordering of two algorithms is independent of the non-considered ones. 2The second interesting feature has to do with the interpretation of the set of parameters of the probability model. When modeling rankings of size k, the Plackett-Luce model has k parameters, one per algorithm. We will refer to these parameters as weights, w i .

If we represent rankings as σ = (σ 1 , . . . , σ k ) where σ i = j means that the j-th algorithm is ranked in the i-th position, the probability of a given ranking σ under the Plackett-Luce model is:

Pr P L (σ |w) = n i=1 w σ i n j=i w σ j . ( 2
)
In general, the weights can take any strictly positive value. However, if we restrict them such that w i = 1,3 then they can be interpreted as the probability of each algorithm being the best (the top-ranked one), i.e.,

w i = Pr(σ 1 = i).
In what follows, we will briefly review the Bayesian Plackett-Luce model, for more details, the interested reader is referred to [START_REF] Calvo | Bayesian Inference for Algorithm Ranking Analysis[END_REF].

As in any Bayesian model, we identify the three distributions mentioned above: the likelihood function, the prior distribution, and the posterior distribution. In this model, the likelihood function is calculated as the product of the probabilities of the observed rankings (of algorithms) under the Plackett-Luce model. The prior distribution of the weights is modeled with a Dirichlet distribution, which is the generalization of the Beta distribution. Briefly, the Dirichlet distribution models probability distributions over realvalued vectors that sum 1 and, thus, it is particularly suitable to model the prior distribution of the weight vector. By representing the vector of weights as w and the data (the rankings of the algorithms) as R = {σ (1) , . . . , σ (N ) }, the posterior distribution of the weights can be computed as:

Pr(w|R) = Pr(R|w)Pr(w) Pr(R) . (3)
The above equation has no closed-form, but we can easily sample the posterior distribution of weights using Markov Chain Monte Carlo methods.

The inference process is as follows. First, the performance data of the algorithms is collected, considering the scenario (the population) from which we want to draw our conclusions. The results are then transformed into rankings of algorithms (one per run) and these rankings are used to produce a sample from the posterior distribution of weights.

Once the posterior distribution of the weights has been defined, there are several aspects of the rankings' distribution represented by the model. We will focus on the probability of each algorithm being the top-ranked. There are other aspects related with the ranking of the algorithms that can be analyzed using the Bayesian Plackett-Luce model, such as the probability of an algorithm being among the two best algorithms or the expected rank for an algorithm. However, for space reasons we will limit the analyses to the two points mentioned above.

BENCHMARKING AND ANALYSIS PREPARATION

In this section we present the preliminaries for the Bayesian inference analysis -the details behind the actual benchmarking campaign (that is, the source of the generated datasets), as well as the analysis' planning and its associated practical considerations.

Black-Box Discrete Optimization Benchmarking

We build our statistical analyses upon performance data that has been obtained from an independent research thread [START_REF]Anonymous for double-blind reasons[END_REF], which is aimed at identifying suitable benchmark problems for discrete optimization heuristics and at showing that the IOHprofiler environment announced in [START_REF] Doerr | IOHprofiler: A Benchmarking and Profiling Tool for Iterative Optimization Heuristics[END_REF] is capable of handling mid-sized benchmarking requirements. Next, we list the 23 optimization problems that form the current benchmarking suite:

• OneMax (F1) and W-model extensions (F4-F10)

• LeadingOnes (F2) and W-model extensions(F11-F17)

• Harmonic (F3): linear harmonic function

• LABS: Low Autocorrelation Binary Sequences (F18) • Ising-Ring (F19) • Ising-Torus (F20) • Ising-Triangular (F21) • MIVS: Maximum Independent Vertex Set (F22) • NQP: N-Queens problem (F23)
The specified functions have been implemented altogether within the IOHexperimenter. Each function F is assessed over four problem dimensions, n ∈ {16, 64, 100, 625}, yielding altogether 92 (F, n) pairs. Each algorithm is run on 11 different instances of each of these 92 pairs, yielding 1, 012 different runs per each algorithm.

Each run is granted a budget of 100n 2 function evaluations for dimensions n ∈ {16, 64, 100} versus 5n 2 evaluations for n = 625.

The following 11 EAs were benchmarked (see [START_REF]Anonymous for double-blind reasons[END_REF] for details):

(1) gHC: a greedy (1+1)-type hill climber which flips one bit per iteration, following a deterministic order. Fips one random bit on average, but can escape local optima by flipping more than one bit. (4) fGA: The "fast GA" proposed in [START_REF] Doerr | Fast genetic algorithms[END_REF]. The number of bits flipped follow a heavy-tail power-law distribution with exponent β = 1.5. ( 5) (1 + 10) EA: The (1+10) EA with static p = 1/n. Like (1+1) EA but keeps best of ten offspring. (6) (1+10) EA r /2,2r : The two-rate (1+10) EA with self-adjusting mutation rates. Taken from [START_REF] Doerr | The (1 + λ) Evolutionary Algorithm with Self-Adjusting Mutation Rate[END_REF]. [START_REF] Demšar | Statistical Comparisons of Classifiers over Multiple Data Sets[END_REF] (1 + 10) EA norm. : (1 + 10) EA sampling the number of bits to flip from an adaptive normal distribution. Taken from [START_REF] Ye | Interpolating Local and Global Search by Controlling the Variance of Standard Bit Mutation[END_REF]. [START_REF] Doerr | Optimal Static and Self-Adjusting Parameter Choices for the (1 + (λ, λ)) Genetic Algorithm[END_REF] (1 + 10) EA var. : The (1 + 10) EA norm. with variance control, also from [START_REF] Ye | Interpolating Local and Global Search by Controlling the Variance of Standard Bit Mutation[END_REF]. [START_REF] Doerr | From black-box complexity to designing new genetic algorithms[END_REF] (1 + 10) EA log-n. The (1+10) EA with log-normal selfadaptation of the mutation rate proposed in [START_REF] Bäck | Intelligent Mutation Rate Control in Canonical Genetic Algorithms[END_REF]. [START_REF] Doerr | The (1 + λ) Evolutionary Algorithm with Self-Adjusting Mutation Rate[END_REF] (1 + (λ, λ)) GA: single-trajectory genetic algorithm using cross-over [START_REF] Doerr | From black-box complexity to designing new genetic algorithms[END_REF]. We use the variant with self-adjusting parameters analyzed in [START_REF] Doerr | Optimal Static and Self-Adjusting Parameter Choices for the (1 + (λ, λ)) Genetic Algorithm[END_REF]. [START_REF] Doerr | Fast genetic algorithms[END_REF] vGA: A (30, 30) "vanilla" GA (see, e.g., [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF]).

Analysis' Considerations and Planning

We illustrate the application of the Bayesian Plackett-Luce model by analyzing some of the performance data of the aforementioned benchmarking campaign. As described above, the numerical data encompass optimization results for 23 objective functions in 4 dimensions, each benchmarked over 11 instances by 11 competing algorithms. For each of these 23 × 4 × 11 × 11 = 11, 132 runs, the following information is recorded: whenever the algorithm has identified a solution of strictly better quality than the previous-best solution, a new data entry is created which stores the new incumbent fitness value and the number of the fitness evaluations elapsed before sampling it. 4The analyses reported in Section 4 involve performance data of all eleven algorithms. Regarding the functions, we will show how to treat them either individually or collectively (in groups). When considered in groups, we will take two approaches -either considering the results altogether (i.e., aggregating the performances across the eleven instances of each function), or taking the median of these runs. As previously noted, we use the Bayesian Plackett-Luce model to estimate the probability of each algorithm being top-ranked.

In practice, we will estimate the probability of an algorithm being the best among its competitors as its expected weight in the posterior distribution of weights. Additionally, we will analyze the uncertainty about these probabilities by estimating the 90% credible intervals for which we will use the corresponding quantiles (5% and 95%) of the posterior distribution of weights.

Finally, we address one of the main concerns with Bayesian methods, which lies in the fact that they make use of a prior distribution for the parameters. Usually, this prior belief is set as a non-informative one, to which we assign a relatively low strength. During the Bayesian inference analysis the prior belief and the information within the data are merged. The impact of the prior belief decreases as the amount of available data increases. We will illustrate and discuss the effects of the prior distribution on the final estimations in Section 4.3.

INFERENCE ANALYSES AND RESULTS

We begin our report by providing a qualitative summary of the observations obtained from a manual examination of the benchmark data. This report is followed by a rigorous, quantitative Bayesian analysis that is meant to address some of the qualitative statements with the goal to corroborate or to refute them.

As an illustration of the data visualization used to support the qualitative analysis, we depict in Fig. 1 the Expected Running Time (ERT) values for all 11 algorithms and all 23 functions in dimension n = 625. The ERT values are with respect to the best solution found by any of the (algorithm, instance) pairs. Importantly, the nature of the ERT map could substantially change upon setting alternative target values. Also note that the different instances of each problem adhere altogether to the same fitness landscape, which just undergoes isomorphic transformations (rotation, scaling, etc., see [START_REF]Anonymous for double-blind reasons[END_REF] for details). For the purpose of our analysis, we can thus regard the data for the eleven instances as eleven independent runs per each algorithm on a given problem. The following list summarizes some of the conclusions drawn in [START_REF]Anonymous for double-blind reasons[END_REF].

• The following class of algorithms exhibits similar performance profiles over the majority of problems: the (1 + (λ, λ))-GA, (1+1)-EA, (1+10)-EA var. , (1 + 10) EA log-n. , (1+10)-EA, (1+10)-EA norm. , (1+10)-EA r /2,2r , and (1+1)-fGA. These algorithms behave quite consistently, typically exhibiting fine performance. • The gHC and the vGA usually exhibit extreme performance with respect to the other algorithms. The vGA consistently suffers from poor performance over all functions, while the gHC either leads the performance on certain functions or shows deteriorated performance on others. • When counting the number of runs in which a solution of overall best target value has been identified per algorithm, the (1+10)-EA r /2,2r leads the hitting scores on the "lowdimensional" functions (n = 16, 64), with the (1+1)-fGA and

Fixed-Targets Perspective

As a first analysis, we will demonstrate algorithms' ranking by considering the results of a particular function in a certain dimension and upon setting a fixed target value φ, denoted as a triplet (F, n, φ). This subset of the data defines running results for the 11 algorithms, and in particular, it contains the recorded running time (i.e., number of evaluations) consumed by each algorithm to hit the prescribed target. Then, for each run, we rank the algorithms according to this running time -where the top-ranked algorithm is defined as the first algorithm to hit the target. Importantly, the distribution that undergoes inference is the rankings we obtain in random runs of the 11 algorithms for the selected (F, n, φ)-triplet. Fig. 2 depicts two examples of the credible intervals and expected values obtained for the probability of each algorithm being the best. These plots illustrate how the Bayesian method does not only provide us with the expected probability, but also with information about its associated uncertainty. Taking into account that we have only 11 samples to do inference, we should expect relatively high uncertainty about the probability of winning. Clearly, this is the case for certain algorithms, such as the two best algorithms in the first example (Fig. 2 (a)). However, for those algorithms that have the smallest probability of winning, 11 samples can be enough to get a quite tight interval for the probability. In general, it is hard to say how many runs are needed to get confident results. However, in Bayesian statistics we can empirically observe how the uncertainty decreases as the number or runs is increased -one of the strong points where the Bayesian inference approach shows its potential.

Another key argument in favor of using Bayesian methods is that they provide us with interpretable information. As an example, attending to the results produced by the analysis above we can state that if we run all the 11 algorithms on F17, dimension n = 625 and fix the target to φ = 625, then the estimated expected probability of RLS being the fastest algorithm hitting the target is 0.426, which is almost 15 times higher than the expected probability for gHC (that reads 0.03). Moreover, the results above estimate that, with probability 0.9, the true probability of RLS being the fastest algorithm lies in the interval [0.272, 0.585]. Statistical tests can also be used to estimate and then present uncertainty in the form of confidence intervals. Although this is not the common practice when comparing algorithms, a binomial test could be used to get intervals for the probability of being the best, but comparing this approach with the proposed Bayesian method is out of the scope of this work. Importantly, mind should be given to the different interpretation of these two approaches (see [16, p.343]).

The plots in Fig. 2 show that all the algorithms have non-zero probabilities of being the best, even if some of them never arise as best in the data. One could think that this is due to the prior distribution, which assigns equal probabilities to all the algorithms before the data is incorporated into the model. However, the reason is not rooted in the prior, which has a negligible impact, as will be demonstrated in Section 4.3. Rather, the Plackett-Luce model itself is responsible for this behavior, since its parameters (i.e., the weights associated with each algorithm), reflect the relative ordering of all the algorithms, and not only the probability of being the best.

So far we have focused on individual functions and, therefore, our conclusions are limited to the analyzed functions. In some situations, we may have a set of instances about which we want to do inference. In such cases, our sample will be the results obtained by the algorithms in a sample of instances. For each instance we can have different runs of the algorithms and, thus, there are two possible paths to follow. The common practice is to average the results of the different runs -yielding a single value for each instance and algorithm. The simplest way of averaging the runs is by taking the mean value, but due to the high risk of not having unimodal, symmetric distributions for the data, a safer approach is to take the median. An alternative way consists in taking all the runs at the same time.

Both strategies are sensible, but they are not equivalent, as the interpretation of the results changes. The difference in the interpretation is rooted in the underlying distribution behind the analyzed samples. Let us assume that the instances at hand constitute a random sample from that statistical population. 5 Then, upon taking the average value, the underlying statistical population becomes the average results obtained by the algorithms in a random instance from the population of instances. Conversely, if all values are accounted for, the underlying population becomes the result of the algorithms in random runs over random instances (again, having exactly the same number of repetitions for each instance is not a proper sample, but the results may be close enough to those obtained with a proper sample).

We will illustrate the differences between these two approaches upon comparing the algorithms' performance over two subsets of the functions at n = 625, corresponding to an "easy" subset of 5 In practice, the instances are often not randomly sampled but rather selected based on some criteria. Formally, this is not a correct way of proceeding if the goal is doing inference. Nevertheless, in many practical situations defining a statistical population of instances is not feasible and, thus, we can take the selected instances as a sample from a certain population. functions (F ∈ {1 -6, 11 -13, 15, 16}) versus the complementary subset of "non-easy" functions. Fig. 3 shows the results for the "easy" functions. Clearly, the best performing algorithm in both cases is gHC. When we analyze random runs in random instances (Fig. 3 (a)), we can see that there is a big difference between gHC and the rest. Evidently, the credible intervals are quite tight, suggesting a considerable reduction of the uncertainty. This is, mostly, because we have a larger sample than in the second case (11 samples versus 1 per function). Fig. 4 shows the results for all the runs in the "non-easy" subset of functions.

In this case we get a somewhat different picture: the estimation of the probabilities is fairly good (the credible intervals are smaller than 0.05, whereas mind should be given to the plot's scaling), but the probability values are fairly similar (simply due to overlapping). Overall, there is an uncertainty about which algorithm is best performing over this subset of functions, but not due to limitation of data, rather due to equivalence in the algorithms' performance.

Fixed-Budget Perspective

A different perspective to empirically analyze the algorithms' runs considers comparisons that are drawn per a fixed budget. The proposed Bayesian method can be used to compare the algorithms also in this scenario. In order to show a different approach to the analysis, we will demonstrate the method by tracking the evolution of the probability of each algorithm being the best as the budget increases, which is a novel perspective to the best of our knowledge. Fig. 5 presents two examples of such probabilities' evolution plots, for different cases. In Fig. 5(a), showing results for F21, it is evident that the behavior of one of the algorithms, gHC, differs from the rest. Its probability of winning quickly increases, showing that for low budgets it is clearly the best algorithm. However, although fine objective function values are quickly acquired by gHC, other algorithms obtain better values as the budget increases, and accordingly, their probability of winning increases, thus reducing the probability of gHC being the best. The picture is quite different when analyzing F22, as depicted in Fig. 5(b): the gHC exhibits steady increase in its probability of being the best, up to a probability of about 50% (a high value, taking into account that we consider altogether 11 algorithms).

Impact of the Prior Distribution

One of the main concerns about Bayesian statistics is the need to define a prior distribution for the parameters of the model. For this reason, an important aspect to investigate is the impact of the prior distribution over the final estimations. In this section we will illustrate this type of analysis using one of the previous examples. In particular, we will compare the posterior distribution of the weights obtained in the analysis of the results for function F 9 at n = 100 with a target value φ = 100 (see Fig. 2).

To this end, we have used three different types of prior distributions: (i) equal treatment of all the algorithms (uniform prior), (ii) preferred probability to the best algorithms (assigning the best algorithm higher prior probability of winning), and (iii) preferred probability to the worst-performing algorithms (deceptive prior). Distributions (ii) and (iii) make use of the empirical results to determine which algorithms have a higher prior probability of being the best. Fig. 6 shows the results of the comparison as boxplots of the posterior sample of weights. Evidently, when the empirical data is used to bias the prior distribution toward the best algorithms, the two best algorithms have slightly higher values and the opposite happens when the deceptive prior is used, as they have slightly smaller values. In any case, the actual differences are practically negligible, even if we have a relatively small sample (11 runs). These small differences are translated into small differences in our estimations, but they definitely will not alter the general comparison overview.

Although the common practice is to use non-informative priors, one nice property of Bayesian statistics is that it provides a natural way to do "incremental" analyses. That is, rather than using a noninformative prior to start the comparison of the algorithms from scratch, we can use previous analyses to set an "informed" prior, thereby testing the algorithms in a loop fashion.

Comparison with Classical Statistical Tests

Finally, as an illustration of the differences between Bayesian inference and the classical statistical test approach, we will show the results obtained with a standard statistical test workflow for the two largest budgets in Fig. 5. All the pairwise comparisons are tested using a Wilcoxon test and the p-values are corrected for multiple testing using Shaffer's method [START_REF] García | An extension on 'statistical comparisons of classifiers over multiple data sets' for all pairwise comparisons[END_REF]. The results are displayed in Fig. 7. The plots are similar to the Critical Differences plot presented in [START_REF] Demšar | Statistical Comparisons of Classifiers over Multiple Data Sets[END_REF]. Each algorithm is displayed according to its average rank in the data. Then, the algorithms which, according to the test, do not significantly differ (corrected p-value above 0.05) are linked by a horizontal line.

The results shown in these plots, although reflecting a different perspective on the data, match those observed in Fig. 5. For the two functions, the ordering of the algorithms is roughly the same in both approaches. Regarding the uncertainty, it is too high for function F21 to conclude that any of the algorithms is clearly the best. Conversely, for F22 both analyses show that gHC is clearly better than the rest of the algorithms (it has significant differences with the second one and the credible intervals in Fig. 5 do not overlap).

In the case of the Bayesian analysis we have the estimated probability of each algorithm being the best, together with a notion (in the form of credible interval) of the uncertainty about this estimation. Moreover, the Bayesian analysis can also be used to estimate the expected rank (similar, in this example, to the ones directly derived form the data; these calculations are not included herein for reasons of space).

DISCUSSION AND OUTLOOK

An important task in Computer Science in general, and in EC in particular concerns the evaluation and comparison of algorithms. Two critical steps in such an assessment comprise an appropriate design of experimentation as well as a proper statistical analysis of the performance data. Different statistical procedures and methods have been suggested in the literature to gain insight into particular aspects of algorithms' performance. In this work we have demonstrated that Bayesian inference using the Plackett-Luce model facilitates an effective analysis of algorithms' performance ranking. Since these results nicely complement classical statistical tests, we suggest to include Bayesian inference with the Plackett-Luce model as a standard statistical tool in the practical EC performance comparisons' repertoire. In this vein, we also suggest to include a Bayesian inference analysis module in IOHanalyzer, the data-analysis component of IOHprofiler.

The Bayesian Plackett-Luce model described in this work has some strong points (e.g., the ability to handle multiple algorithms), but -as with any statistical tool -also some weaknesses. It is therefore important to understand the limits of the model. An important point to note is that by aggregating performances into rankings we lose information about the magnitude of the differences. Put differently, the model is capable of estimating the probability of an algorithm being the best, but it does not quantify the performance gap with respect to the first runner-up. Developing new statistical procedures to inspect the results from other, alternative points of view, forms an important line for future research. Another important line of further research concerns a deeper analysis of the Bayesian Plackett-Luce model -not only regarding aspects of convergence of the estimators, but also in terms of identifying other useful metrics. An example for such metrics is the probability of an algorithm being among the top k algorithms, which has the potential to enrich the algorithms' comparison.

Finally, an ambitious but a necessary involves the careful definition of the most appropriate workflow for the comparisons. There is no single method that is best suited to answer all the questions we may ask, and not all the methods require the same type of data. Thus, clearly articulating the questions, devising the appropriate methods to address them, and formulating the data's requirements for their application -are tasks from which the community can greatly benefit.

( 2 )

 2 RLS: Randomized Local Search. Like gHC but flips one randomly chosen bit per iteration. (3) (1 + 1) EA: The (1+1) EA with static mutation rate p = 1/n.

Figure 1 :

 1 Figure 1: ERT values of all 11 algorithms for the 625-dimensional test suite, with respect to the best solution quality found by any of the algorithms in any of the eleven runs.

  F 9, n = 100, φ = 100

Figure 2 :

 2 Figure 2: Credible intervals (5% and 95% quantiles) and expected probability of winning (green dots) for two combinations of function, dimension and target.

Figure 3 :Figure

 3 Figure3: Credible intervals (5% and 95% quantiles) and expected probability of winning (green dots) for the "easy" subset of functions at n = 625. Plot (a) represents the results using all the runs, whereas plot (b) has been obtained with the median of the runs.

  F 21, n = 100 (b) F 22, n = 100

Figure 5 :

 5 Figure 5: Evolution of the probability of winning for different budgets. The lines represent the evolution of the expected probability of winning while the error bars represent the 90% credible intervals (5% and 95% quantiles).

Figure

  Figure Comparing the posterior distribution of weights, when obtained with different prior distributions, for the data analysis corresponding to F 9, n = 100, φ = 100.

  F 22, n = 100, Budget 200

Figure 7 :

 7 Figure7: Demonstration of the statistical test workflow applied to the results for two problem instances for a fixed budget. The plots present the average rank of each algorithm, and also display which algorithms do not exhibit statistically significant differences (linked by a horizontal line).

The careful reader should note that the term 'population' is exclusively used herein to describe the underlying population of the analysis, rather than the common use in EC referring to a set of candidate solutions.

To ease the reading, instead of using terms such as item or element to name the components of the permutations, from now on, we will refer to them as algorithms.

Note that, due to the way the probabilities are computed, we can force the sum to be any positive value without loss of generality, as any common factor for the parameters will not affect the probabilities.

More precisely, IOHexperimenter tracks much more data, but all computations in this work are solely based upon the here-described records.

ACKNOWLEDGMENTS

This work has been partially supported by the BERC 2018-2021 and ELKARTEK programs (Basque Government), and the project TIN2016-78365-R and Severo Ochoa Program SEV-2017-0718 (Spanish Ministry of Economy, Industry and Competitiveness). Research was also supported by COST Action CA15140 "Improving Applicability of Nature-Inspired Optimisation by Joining Theory and Practice (ImAppNIO)".


Carola Doerr

Johann Dreo

Pascal Kerschke

Making a Case for (Hyper-)Parameter Tuning as Benchmark Problems -A Discussion Paper

Keywords:

One of the biggest challenges in evolutionary computation concerns the selection and configuration of a best-suitable heuristic for a given problem. While in the past both of these problems have primarily been addressed by building on experts' experience, the last decade has witnessed a significant shift towards automated decision making, which capitalizes on techniques proposed in the machine learning literature.

A key success factor in automated algorithm selection and configuration are good training sets, whose performance data can be leveraged to build accurate performance prediction models. With the long-term goal to build landscape-aware parameter control mechanisms for iterative optimization heuristics, we consider in this discussion paper the question how well the 24 functions from the BBOB test bed cover the characteristics of (hyper-)parameter tuning problems. To this end, we perform a preliminary landscape analysis of two hyper-parameter selection problems, and compare their feature values with those of the BBOB functions. While we do see a good fit for one of the tuning problems, our findings also indicate that some parameter tuning problems might not be very well represented by the BBOB functions. This raises the question if one can nevertheless deduce reliable performance-prediction models for hyper-parameter tuning problems from the BBOB test bed, or whether for this specific target the BBOB benchmark should be adjusted, by adding or replacing some of its functions.

Independently of the aspect of training automated algorithm selection and configuration techniques, hyper-parameter tuning problems offer a plethora of problems which might be worthwhile to study in the context of benchmarking iterative optimization heuristics.

INTRODUCTION

Evolutionary algorithms, as well as all other commonly used search heuristics, are, strictly speaking, parametrized frameworks that the user needs to instantiate before execution. The instantiation typically requires to set a number of parameters, such as the population size, search radius, selective pressure, etc. The parametrization offers a great flexibility, which allows the user to adjust the algorithm's behavior to the needs of her problem. However, this flexibility comes at the cost of creating another optimization problem on top: parameter tuning.

Parameter tuning was very early identified as a very important problem [START_REF] Grefenstette | Optimization of Control Parameters for Genetic Algorithms[END_REF], since the performance of search heuristics is typically highly sensitive to their parameters [START_REF] Charles | The Choice of Step Length, a Crucial Factor in the Performance of Variable Metric Algorithms[END_REF]. To date, the parameter tuning problem is considered to be one of the most important challenges in evolutionary computation (and beyond) [START_REF] Bartz-Beielstein | Experimental Methods for the Analysis of Optimization Algorithms[END_REF][START_REF] Burke | Hyper-Heuristics: A Survey of the State of the Art[END_REF][START_REF] Feurer | Hyperparameter Optimization[END_REF][START_REF]Parameter Setting in Evolutionary Algorithms[END_REF].

Non-surprisingly, significant research efforts have been undertaken to assist the user in the parameter tuning problem. These works have stimulated the design of automated tools which solve the tuning problem for a given algorithm (on a given sub-problem), see [START_REF] Ansótegui | Model-Based Genetic Algorithms for Algorithm Configuration[END_REF][START_REF] Bartz-Beielstein | SPOT: A Toolbox for Interactive and Automatic Tuning in the R Environment[END_REF][START_REF] Hutter | Sequential Model-Based Optimization for General Algorithm Configuration[END_REF][START_REF] Li | Hyperband: A Novel Bandit-Based Approach to Hyperparameter Optimization[END_REF][START_REF] López-Ibáñez | The irace package: Iterated Racing for Automatic Algorithm Configuration[END_REF] for only a few examples. In a number of complementary workstreams the question of how to control good parameter values has been addressed. The underlying observation of parameter control is that -for a very significant fraction of algorithms and many problems -the range of good parameter values can change quite drastically during the optimization process. Parameter control is therefore not only aiming to identify optimal parameter values, but also to track their evolution throughout the whole run. Paradoxically, parameter control mechanisms are again parametrized, so that one might tend to conclude that not much can be gained by such approaches. Note, however, that they gain us the flexibility to adjust the first-order parameters during the run. Numerous examples are known where parameter control outperforms static choices, see [START_REF] Doerr | Theory of Parameter Control Mechanisms for Discrete Black-Box Optimization: Provable Performance Gains Through Dynamic Parameter Choices[END_REF][START_REF] Karafotias | Parameter Control in Evolutionary Algorithms: Trends and Challenges[END_REF] for recent surveys. The parameters of the control mechanism are commonly referred to as hyper-parameters, and when considering their optimization we speak of hyper-parameter tuning.

Giving the omnipresent use of search heuristics throughout all industrial branches and all data-driven scientific disciplines, we therefore see that the parameter tuning and hyper-parameter tuning problems constitute tasks for which a solid understanding is at very high demand.

Unfortunately, (hyper-)parameter tuning problems can be very complex. Examples in which a small change in the parameter value results in an exponential performance gap can be found in [START_REF] Doerr | Mutation Rate Matters Even When Optimizing Monotonic Functions[END_REF][START_REF] Lengler | A General Dichotomy of Evolutionary Algorithms on Monotone Functions[END_REF]. [START_REF] Lengler | A General Dichotomy of Evolutionary Algorithms on Monotone Functions[END_REF] shows an example for which an exponential performance gap is caused by a threshold behavior, i.e., there is an abrupt change in performance at a certain parameter value. On the other hand, it has also been demonstrated that a number of parameter tuning problems are actually much more benign (e.g., smooth, convex, and/or unimodal) than one might have expected in the first place [START_REF] Pushak | Algorithm Configuration Landscapes: -More Benign Than Expected?[END_REF]. This diversity of behaviours calls for more studies on the understanding of parameter tuning problems, apart from purely performance-oriented studies.

With this discussion paper, we are aiming at starting a focused discussion on how well (hyper-)parameter problems are represented in common benchmarks, and whether there is need and room to extend these benchmarks by (hyper-)parameter tuning problems. With this goal in mind, we provide a preliminary study in which we evaluate the fitness landscape of two hyper-parameter tuning problems taken from the evolutionary computation literature, and compare their structure to those of the 24 noiseless functions of the BBOB test bed [START_REF] Hansen | COCO: A Platform for Comparing Continuous Optimizers in a Black-Box Setting[END_REF].

Methodology: Exploratory Landscape Analysis. Exploratory landscape analysis (ELA) [START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems By Combining Exploratory Landscape Analysis and Machine Learning[END_REF][START_REF] Mersmann | Exploratory Landscape Analysis[END_REF] is a way to measure the structure of a problem f by sampling from it some (typically random or quasi-random) solutions, evaluating these, and mapping the resulting set of (x, f (x))-pairs to a vector of real numbers that are meant to capture the important characteristics of the problem. The objective is to know what makes the problem difficult, so as to adapt the algorithm's behaviour. The landscape-aware parameter tuning approach is thus to observe the problem through the estimation of its landscape features and use them to tune the algorithm's parameters to optimize its performance.

This approach led to promising results when trying to find the best suitable hyper-parameter tuning to solve a given problem instance. In this "per instance algorithm configuration" [START_REF] Leyton-Brown | Learning the Empirical Hardness of Optimization Problems: The Case of Combinatorial Auctions[END_REF] approach, a set of features is first measured on a random sample and the parameter tuning is derived from a previously learned mapping [START_REF] Hutter | Performance Prediction and Automated Tuning of Randomized and Parametric Algorithms[END_REF]. This approach necessitates to learn the landscape-parameter mapping, called an empirical performance model [START_REF] Belkhir | Feature Based Algorithm Configuration: A Case Study with Differential Evolution[END_REF], on an existing benchmark.

In such a setting, the efficiency of the landscape-aware parameter tuning strongly depends on the quality of the benchmark. More precisely, it depends on the so-called "footprint" [START_REF] Kerschke | Automated Algorithm Selection: Survey and Perspectives[END_REF][START_REF] Andrés | Performance Analysis of Continuous Black-Box Optimization Algorithms via Footprints in Instance Space[END_REF] of the benchmark instances across the features space. If the benchmark does not exhibit (some of) the landscapes of the real-world problems, then the performances of the adapted algorithm may actually be worse than a more generic counterpart.

Thus, the most urgent question addressed by this paper reduces to: Are landscapes of algorithm configuration problems similar to those observed in benchmarks that are commonly used for learning?

In this preliminary work, we use the Black-Box Optimization Benchmark (BBOB) [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF] as a reference. The noiseless functions of BBOB have been used as the learning benchmark in previous work on landscape-aware per-instance algorithm configuration [START_REF] Belkhir | Feature Based Algorithm Configuration: A Case Study with Differential Evolution[END_REF], which led to good performance improvements [START_REF] Belkhir | Per Instance Algorithm Configuration of CMA-ES with Limited Budget[END_REF]. However, it should be noted that the landscape-aware parameter tuning problem is in itself a noisy problem involving a machine learning process. As such, the recently proposed "Nevergrad" benchmark [START_REF] Rapin | Exploring the MLDA Benchmark on the Nevergrad Platform[END_REF][START_REF] Rapin | Nevergrad -A Gradient-Free Optimization Platform[END_REF] may be also considered.

Disclaimer: While we present two case studies within this discussion paper, we do not claim that they are representative. Our (1) , . . . , x (λ 1 ) } with f (x ′ ) = max{ f (x (1) ), . . . , f (x (λ 1 ) )} u.a.r.;

6 for i = 1, . . . , λ 1 = nint(λ) do x (i) ← flip ℓ (x); 7 Choose x ′ ∈ {x
8 Crossover phase:

9 for i = 1, . . . , λ 2 = nint(βλ) do y (i) ← cross c=γ /λ (x, x ′ ); 10
Choose y ∈ {x ′ , y (1) , . . . , y (λ 2 ) } with f (y) = max{ f (x ′ ), f (y (1) ), . . . , f (y (λ 2 ) )} u.a.r.;

11 Selection and update step: To investigate how well hyper-parameter tuning fits into the BBOB test bed, we describe in this section a first set of experiments, in which we compare the feature values of a hyper-parameter tuning problem with those of the 24 noiseless BBOB functions. More precisely, we consider the fitness landscape of the five-dimensional hyper-parameter tuning problem of the generalized (1 + (λ, λ)) GA presented in [START_REF] Dang | Hyper-Parameter Tuning for the (1 + (λ, λ)) GA[END_REF]. We tune its performance on the 5,000-dimensional OneMax problem. The generalized (1 + (λ, λ)) GA is summarized in Algorithm 1, we do not discuss it in great detail here, and point the interested reader to [START_REF] Dang | Hyper-Parameter Tuning for the (1 + (λ, λ)) GA[END_REF][START_REF] Doerr | Optimal Static and Self-Adjusting Parameter Choices for the (1 + (λ, λ)) Genetic Algorithm[END_REF][START_REF] Doerr | From Black-Box Complexity to Designing New Genetic Algorithms[END_REF] for an in-depth discussion of this algorithm. For the tuning problem we consider the following ranges

12 if f (y) > f (x) then x ← y; λ ← max{bλ, 1}; 13 if f (y) = f (x) then x ← y; λ ← min{Aλ, n -1}; 14 if f (y) < f (x) then λ ← min{Aλ,
: α ∈ [0.33, 10], β ∈ [1, 10], γ ∈ [0.333, 10], A ∈ [1.01, 2.5],
and b ∈ [0.4, 0.99]. For the BBOB test bed, we used the R-package smoof [START_REF] Bossek | smoof: Single-and Multi-Objective Optimization Test Functions[END_REF] and considered the first five instances per function. A detailed description of the BBOB functions can be found in [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF], and a survey on the COCO benchmarking platform, which BBOB is a part of, is available at [START_REF] Hansen | COCO: A Platform for Comparing Continuous Optimizers in a Black-Box Setting[END_REF].

For each of the 120 BBOB instances (24 functions with five instances each), as well as for the hyper-parameter tuning problem, we used a low-discrepancy point set generator to sample 3 5 = 243 points x in total at which we evaluated the fitness. More precisely, we have used a Halton point set for the generation of the configurations at which we evaluate the (1 + (λ, λ)) GA tuning problem [START_REF] De Rainville | Evolutionary Optimization of Low-Discrepancy Sequences[END_REF], and we use a Latin Hypercube sample for the BBOB evaluations (for convenience, as it is already implemented as standard sampler in flacco [START_REF] Hanster | flaccogui: Exploratory Landscape Analysis for Everyone[END_REF][START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems By Combining Exploratory Landscape Analysis and Machine Learning[END_REF], the tool that we use to compute the feature values, see below). For the hyper-parameter tuning problem we scale the [0, 1] d -distributed points from the generator to the above-specified range of the parameters α, β, γ , A, b. For each of the 243 configurations we run the self-adjusting (1 + (λ, λ)) GA eleven times, average the optimization times of the resulting eleven values, and assign this average running time as fitness to the respective configuration.

The heatmaps in Figure 1 give a first indication how the sampled configurations are distributed, and how the fitness landscape of the five-dimensional tuning problem looks like, in a global perspective (upper figure) and for a zoom into the configurations achieving an average optimization time of at most 10 7 function evaluations (bottom figure). We note here that none of the evaluated configurations performs extremely well: the best configuration has an average running time of around 43,395 evaluations, whereas in [START_REF] Dang | Hyper-Parameter Tuning for the (1 + (λ, λ)) GA[END_REF] several configurations are reported to have a much smaller average optimization time, which can go as low as 29,000 function evaluations. We also note that the plots indicate some regularity, despite the fact that the Halton point set was meant to provide a good distribution. It might be worth investigating how much our results would change depending on the sample distribution. The point sets from [START_REF] De Rainville | Evolutionary Optimization of Low-Discrepancy Sequences[END_REF] are known to achieve low L2-discrepancy, but other measures such as the star discrepancy or energy level measures might be more suitable. This, however, forms a separate research thread which we ignore in this present work.

For each of the 121 sets of 243 (x, f (x))-pairs we then compute the feature values using flacco [START_REF] Hanster | flaccogui: Exploratory Landscape Analysis for Everyone[END_REF][START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems By Combining Exploratory Landscape Analysis and Machine Learning[END_REF]. Since some feature computations are stochastic, we perform ten independent runs per data set, leaving us with 1,200 feature vectors for the BBOB functions and ten feature vectors for the hyper-parameter tuning problem.

For a first indication how similar the tuning problem is to those in the BBOB testbed, we next determine minimum and maximum for each of the 66 features (we here disregard all basic features, as well as cost-related subfeatures as they are not related to the fitness landscape itself). We then evaluate if the feature values of the tuning problem are within the range of the BBOB functions. As a result, we obtain that all but four feature values fall in the BBOB range. Figure 2 depicts the distribution of the conspicuous feature values for the 120 BBOB instances against the value of the tuning problem. All four features are based on deterministic computations, and the distribution is therefore only with respect to the 120 deterministic feature values.

The idea behind the four conspicuous features mentioned above can be summarized as follows: the cell mapping angle (cm_angle) features [START_REF] Kerschke | Cell Mapping Techniques for Exploratory Landscape Analysis[END_REF] discretize the search space into three equidistant intervals per dimension, resulting in 3 5 = 243 cells in total. For all non-empty cells, i.e., cells which contain at least one of the sampled observations, the angle (in degree) between the worst, the center and the best observation of the respective cell is computed. Afterwards, the angles of all feasible (i.e., non-empty) cells are aggregated by means of the arithmetic mean (angle.mean) and their standard deviation (angle.sd). For the nearest better clustering (nbc) features [START_REF] Kerschke | Detecting Funnel Structures by Means of Exploratory Landscape Analysis[END_REF], two distance sets are computed. The first distance set comprises the distances of all points to their respective nearest neighbors, whereas the second one contains all nearest better neighbor distances, i.e., the distances to the nearest neighbor among all observations with a better fitness value. nn_nb.mean_ratio computes the ratio between the mean distances of the two distance sets. In contrast, dist_ratio.coeff_var first computes the distance ratio per observation and afterwards aggregates these ratios by means of their coefficient of variation (arithmetic mean divided by standard deviation).

Averaged Feature Analysis

We perform a second analysis in which we average the feature values across the ten independent runs of flacco and, in case of the BBOB functions, the five instances. This approach reveals seven (feature) outliers, whose distributions are summarized in Figure 3. The features listed therein extend the previous four outliers from Figure 2 by one gradient homogeneity (cm_grad, [START_REF] Kerschke | Cell Mapping Techniques for Exploratory Landscape Analysis[END_REF]) feature and two y-distribution (ela_distr, [START_REF] Mersmann | Exploratory Landscape Analysis[END_REF]) features. For the former, we first compute the (normalized) gradients between each point in a cell and its nearest neighbor, then point each gradient towards its better end, and afterwards sum up all gradient vectors (per cell). At last, the lengths of all cumulated gradient vectors are aggregated by means of the standard deviation (sd). In contrast, the y-distribution features simply aggregate the frequency of the fitness values. More precisely, kurtosis and skewness are the kurtosis and skewness of all fitness values f (x) from the sample.

Tuning vs BBOB Functions

While we have merely looked at aggregated data in the analyses above, we next compare the feature values of the tuning problem to each of the individual BBOB functions. To this end, we compute for each of the 24 BBOB functions the vector of average feature values, and compare them to the average feature values of the tuning problem. We recall here that for the tuning problem, the averaging only concerns the randomness stemming from the feature computation, whereas for the BBOB functions we also aggregate over the five instances. Figure 4 displays the results of this comparison. While we see a seemingly good fit of the ELA feature values of the tuning problem with those of BBOB's F1, F3, F8, F13, F14 and F15 (see Figure 4b), no BBOB function could be identified for a good fit with respect to the other feature sets displayed in Figure 4a and4c.

The plots in Figure 4 also raise the interesting question if one can design an aggregated measure that captures the similarity between two or more fitness landscapes. To this end, a more solid understanding of the relevance of each feature needs to be developed. OneMax (< 10,000,000) new offspring y is created by flipping each bit of the parent x independently with mutation rate p (standard bit mutation). When y = x a new bit string is sampled until we obtain an offspring y x.

Equivalently, we can sample the mutation strength, i.e., the number of (randomly distributed) bits to be flipped, from the conditional binomial distribution Bin >0 (n, p)(k). This distribution evaluates to 0 for k = 0 and to Bin(n, p)(k)/(1

-(1 -p) n ) = n k p k (1 -p) n-k /(1 - (1 -p) n ) for 1 ≤ k ≤ n.
This is the description used in Algorithm 2 and explains the subscript " >0 " in the name of the algorithm. Once an offspring y x is sampled, the algorithm evaluates the fitness f (y) and replaces the parent x by y if and only if f (y) ≥ f (x). In this case the mutation rate p is increased to Ap, where A > 1 is one of the three hyper-parameters of this algorithm. When f (y) < f (x) the offspring y is discarded and the mutation rate lowered to bp, where b < 1 is the second hyper-parameter. The search process continues until a user-defined stopping criterion is met.

The self-adjusting (1 + 1) EA >0 has thus three hyper-parameters: the update strengths A and b and the initial mutation rate p 0 ∈ (0, 1). It has been demonstrated in [START_REF] Doerr | Sensitivity of Parameter Control Mechanisms with Respect to Their Initialization[END_REF] that the influence of the initial mutation rate is negligibly small when considering "easy" optimization problems such as OneMax and LeadingOnes. We therefore fix in the following p 0 = 1/n and concentrate on the twodimensional tuning problem A ∈ (1, 6] and b ∈ (0, 1). This is the tuning problem that was studied in [START_REF] Doerr | On the Effectiveness of Simple Success-Based Parameter Selection Mechanisms for Two Classical Discrete Black-Box Optimization Benchmark Problems[END_REF], with a purely performanceoriented mindset. Here, instead, we are rather interested in the characteristics of its fitness landscape.

We consider as test case the (1 + 1) EA >0 on the 250-dimensional LeadingOnes problem Lo : {0, 1} n → [0..n], x → max{i ∈ [0..n] | ∀j ∈ [i] : x j = 1}. In this section, all numerical evaluations for the tuning problem are built upon the data from [START_REF] Doerr | On the Effectiveness of Simple Success-Based Parameter Selection Mechanisms for Two Classical Discrete Black-Box Optimization Benchmark Problems[END_REF]. This data is as follows: For each configuration x = (A, b) with A ∈ {1 + 0.1k | 1 ≤ k ≤ 50} and b ∈ {0.02k | 1 ≤ k < 50} we assign as f (x)-value the average running time of 101 independent runs of the self-adjusting (1 + 1) EA >0 with update strength A and b. We note that this results in a total number of 2,450 (x, f (x))-pairs for the tuning problem. For the BBOB data, however we build our computations on 50D = 100 points only (for the simple reason that we had this data readily available -an update with a similar-sized sample will be done for the camera-ready version, should this paper be accepted).

Alternative Parameter Representations

We briefly note that a recent work [START_REF] Doerr | Self-Adjusting Mutation Rates with Provably Optimal Success Rules[END_REF] shows that a better parametrization of A and b would require that b = (1/A) 1/(s-1) , so that A is the update strength and s the success rule of the algorithm. This setting generalizes the well-known one-fifth success rule from evolution strategies [START_REF] Devroye | The Compound Random Search[END_REF][START_REF] Rechenberg | Evolutionsstrategie[END_REF][START_REF] Schumer | Adaptive Step Size Random Search[END_REF]. It was proven in [START_REF] Doerr | Self-Adjusting Mutation Rates with Provably Optimal Success Rules[END_REF] that the optimal update strength F for the (1+1) EA >0 on the LeadingOnes problem satisfies F = 1 + o(1), while the optimal success rule is around 1.285. With this setting, the expected optimization time on the n-dimensional LeadingOnes function is around 0.404n 2 (to be very precise, the latter is based on a numerical evaluation of an otherwise rigorous mathematically proven bound). Since such a parametrization requires a very deep understanding of the underlying optimization process, and since we want to build our analysis upon the data provided by [START_REF] Doerr | On the Effectiveness of Simple Success-Based Parameter Selection Mechanisms for Two Classical Discrete Black-Box Optimization Benchmark Problems[END_REF], we ignore this alternative parametrization in the following, and consider the tuning of A and b as used in Algorithm 2. We mention this alternative representation because we consider it an interesting question to investigate the influence of such representations on the fitness landscapes and on the difficulty of the underlying tuning problem. We plan to address this question in future work.

Heatmaps

To get a first impression of the global fitness landscape, we first plot again some heatmaps of the tuning problem, cf. Figure 5. The left figure displays the full data set, whereas the ones in the middle and on the right zoom into those configurations which obtain an average running time of at most 50,000 and 30,000, respectively. For comparison, the best configurations achieve a running time slightly above 25,000 ≈ 0.403n 2 . We observe that the landscape appears to be very flat when looking at the largest resolution. However, some structure can be observed when zooming into the more performant runs. In principle, it might be possible to extend the results from [START_REF] Doerr | Self-Adjusting Mutation Rates with Provably Optimal Success Rules[END_REF] to formulate the dependence of the average running time of the self-adjusting (1 + 1) EA >0 on the LeadingOnes problem. Based on those results, a rather smooth fitness landscape can be expected. The heatmap structure may indicate that this tuning problem is not too different from the BBOB functions, a question that we will evaluate in more detail in the next subsection.

Feature Value Comparison

Similar as in Section 2 we again first look for feature values of the tuning problem that do not fall into the range of the 120 BBOB instances. This is the case for all but three features, whose distributions are illustrated in Figure 6.

The first feature (dist_ctr2worst.mean) again exploits information based on the 3 5 = 243 cells of the discretized search space. For all (non-empty) cells, it computes the distance between the cell's worst point and center. Afterwards the computed distances are averaged across the non-empty cells. The second feature (nb_fitness.cor) measures the correlation between an observation's fitness value and its "indegree" (i.e., the number of points from the sample for which the current observation constitutes the nearest better neighbor). At last, nn_nb.sd_ratio is the ratio between the standard deviations of the nearest neighbor and nearest better neighbor distance sets.

Detailed Comparison

Finally, we investigate how the feature values of the parameter tuning problem compare to those of the BBOB functions. Based on the data given in Figure 7, it seems that this tuning problem is indeed, as conjectured in the beginning of this section, better represented by the BBOB functions than the problem considered in Section 2. In particular for functions F3 and F15, visual comparison leads us to observe a slightly better fit across various (but not all) feature values.

CONCLUSIONS AND OUTLOOK

With this work we aim to trigger a discussion on the suitability of hyper-parameter tuning as benchmarking problems for optimization algorithms. Using two examples stemming from the theory of EA literature we have given some indication that the fitness landscapes described by hyper-parameter tuning may currently not be very well reflected in the 24 BBOB functions. However, we cannot rule out the possibility that performance data from the BBOB test bed can nevertheless provide discriminating information about the algorithms' performances, in particular by combining data across several benchmark functions. For a proper scientific evaluation, we therefore suggest a more detailed, double-sided study on:

(1) whether BBOB performance data can nevertheless be used as training sets to build reasonable performance models for (hyper-)parameter tuning problems, and/or (2) which parameter tuning problems are particularly appropriate for benchmarking (general-purpose and/or problemspecific) solvers

We recall that an answer to the second question must not ignore the fact that for a reasonable experimental setting the function evaluations of the tuning problem should not be too time-consuming (unless one is interested in the very expensive optimization case, in which one deals with very small evaluation budgets). Finding suitable tuning problems that are fast to evaluate, but at the same time also representative for "real-world" tuning problems might be a challenging task.

Going forward, we plan on addressing both problems listed above by investigating parameter tuning problems of different types, e.g., exhibiting different dimensionalities, different underlying problems, etc. In concrete terms, a candidate tuning problem that we are particularly keen on understanding is the (ideally automated) configuration of the hyper-parameters of the CMA-ES. We furthermore plan on extending the tuning of its hyper-parameters to setting different CMA-ES modules that have been suggested in the research literature. A suitable framework for such a study is provided by the modular CMA-ES framework available from [START_REF] Sander Van Rijn | Evolving the Structure of Evolution Strategies[END_REF].

One of our long-term objectives is the development of landscapeaware techniques to control the selection and the configuration of algorithms, which leverage existing experimental data and appropriate tools from the machine learning literature to select algorithms and parameters on the fly. Besides building a suitable training set, such an approach also requires to build efficient predictors for the fitness landscapes. That is, we need to design mechanisms that require only a very limited amount of additional function evaluations for approximating the feature values in the optimization problem's current state. This question comprises three main difficulties: (1) predicting the feature values based on a small number of samples (see [START_REF] Belkhir | Surrogate Assisted Feature Computation for Continuous Problems[END_REF][START_REF] Kerschke | Low-Budget Exploratory Landscape Analysis on Multiple Peaks Models[END_REF] for a discussion and first approaches), (2) selecting a subset of features that is large enough for landscape discrimination, while being small enough to reduce the learning complexity [START_REF] Kerschke | Low-Budget Exploratory Landscape Analysis on Multiple Peaks Models[END_REF][START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems By Combining Exploratory Landscape Analysis and Machine Learning[END_REF], and

(3) the refinement of feature value computations that are designed to capture the structure of the global optimization problem to a local perspective, which measures how the relevant part of the search space currently "seen by the algorithm" looks like.

PERSPECTIVES

Tuning problems might help to bridge the gap between continuous and discrete evolutionary computation. Independently from the type of sub-problem the search heuristic is solving, the most generic hyper-parameter setting problems involve several parameter types: continuous (e.g., mutation rate), integer (e.g., population size), or qualitative parameters (e.g., choice of modules). Understanding discrete and continuous optimization is therefore required for most real-world applications. We also want to outline the interplay between machine learning and black-box optimization: Recent works on landscape-aware heuristics have shown that machine learning can be used to learn a feature-parameter mapping [START_REF] Belkhir | Feature Based Algorithm Configuration: A Case Study with Differential Evolution[END_REF][START_REF] Belkhir | Per Instance Algorithm Configuration of CMA-ES with Limited Budget[END_REF][START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems By Combining Exploratory Landscape Analysis and Machine Learning[END_REF]. That alone forms an interesting use case for machine learning practitioners, as the corresponding data set is cleanly generated on-demand and at almost any scale, by sampling a well-defined decision space.

In both domains, the hyper-parameter setting problems are very similar and search heuristics are already applied to solve them both. However, while the landscape-aware approach has already been proven to be useful in the black-box optimization context, we are not aware if similar concepts are in use within the machine learning communities.

Algorithm 1 : 4 Mutation phase: 5

 145 The self-adjusting (1 + (λ, λ)) GA variant with five hyper-parameters α, β, γ , A, b 1 Initialization: Sample x ∈ {0, 1} n u.a.r.; 2 Initialize λ ← 1; 3 Optimization: for t = 1, 2, 3, . . . do Sample ℓ from Bin >0 (n, p = αλ/n);

Figure 1 :

 1 Figure 1: Heatmaps depicting the pairwise hyper-parameter landscapes of the "GA on OneMax" tuning problem for two precision values: all configurations (top) and configurations with less than 10 7 function evaluations (bottom).

Figure 2 :Figure 3 :

 23 Figure 2: Features based on 120 BBOB instances (shown as violin plots) in comparison to the respective feature values based on the hyper-parameter landscape of the "(1 + (λ, λ)) GA on OneMax" tuning problem (red horizontal line), for the four features for which the tuning problem does not fall into the range of the BBOB instances.

Algorithm 2 : 5 y 8 x 10 p

 25810 The (1 + 1) EA >0 with update strengths A and b and initial mutation rate p 0 ∈ [1/n 2 , 1/2] for the maximization of a pseudo-Boolean function f : {0, 1} n → R 1 Initialization: Sample x ∈ {0, 1} n uniformly at random and compute f (x); 2 Set p = p 0 ; 3 Optimization: for t = 1, 2, 3, . . . do 4 Sample ℓ from Bin >0 (n, p); ← flip ℓ (x); 6 evaluate f (y); 7 if f (y) ≥ f (x) then ← y and p ← min{Ap, 1/2} 9 else ← max{bp, 1/n 2 }

  cell mapping, information content, and nearest better clustering features

Figure 4 :

 4 Figure 4: Comparison of selected feature values of the tuning problem (bold black line) with those of the 24 BBOB functions (slim red lines). All values have been normalized to the interval [0, 1]

Figure 5 :

 5 Figure 5: Heatmaps depicting the hyper-parameter landscapes of LeadingOnes for different precision values: all configurations (left), as well as configurations with < 50k (middle) and < 30k evaluations (right).

Figure 6 :

 6 Figure 6: Features based on the 24 BBOB functions (shown as violin plots) in comparison to the respective feature values based on the landscape of the LeadingOnes hyperparameter tuning problem (red horizontal line), for the three features for which the tuning problem does not fall into the range of the BBOB functions.

  cell mapping, information content, and nearest better clustering features

Figure 7 :

 7 Figure 7: Comparison of selected feature values of the "(1 + 1) EA >0 on 250d-LeadingOnes" tuning problem (bold black line) with those of selected or all 24 BBOB functions (slim red lines). All values have been normalized to the interval [0, 1]

SECOND CASE STUDY: THE (1+1) EA WITH SELF-ADJUSTING MUTATION RATESWe add a second case study to our findings from Section 2, in which we consider the algorithm proposed in[START_REF] Doerr | On the Effectiveness of Simple Success-Based Parameter Selection Mechanisms for Two Classical Discrete Black-Box Optimization Benchmark Problems[END_REF], a (1+1) evolutionary algorithm (EA) with self-adjusting mutation rates. Its pseudocode can be found in Algorithm 2. The (1 + 1) EA >0 works as follows: after a random initialization, the algorithm always maintains a best-so-far solution x, breaking ties towards the last evaluated search point of current-best fitness value. In each iteration one

ACKNOWLEDGMENTS

Our work is supported by the European Cooperation in Science and Technology through COST Action CA15140. C. Doerr acknowledges support from the Paris Ile-de-France Region. P. Kerschke is supported by the European Research Center for Information Systems (ERCIS).

We thankfully acknowledge the research network COnfiguration and SElection of ALgorithms (COSEAL) and the participants of the Dagstuhl seminar 16412 on Automated Algorithm Selection and Configuration for numerous discussions on the subject of this work.