Mohamed El Yafrani

Marcella Scoczynski

Inkyung Sung

Markus Wagner

Carola Doerr

Peter Nielsen

Mohamed El

MATE: A Model-based Algorithm Tuning Engine -A proof of

Keywords: Parameter tuning, Model-based tuning, Genetic programming

come   L'archive ouverte pluridisciplinaire

Motivation

The performance of many algorithms is highly dependent on tuned parameter configurations made with regards to the user's preferences or performance criteria [START_REF] Belkhir | Feature based algorithm configuration: A case study with differential evolution[END_REF], such as the 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, and so on. This configuration task can be considered as a second layer optimisation problem [START_REF] Hoos | Programming by optimization[END_REF] relevant in the fields arXiv:2004.12750v2 [cs.NE] 15 Feb 2021 of optimisation, machine learning and AI in general. It is a field of study that is increasingly critical as the prevalence of the application of such methods is expanded. Over the years, a range of automatic parameter tuners have been proposed, thus leaving the configuration to a computer rather than manually searching for performance-optimised settings across a set of problem instances. These tuning environments can save time and achieve better results [START_REF] Ansótegui | A gender-based genetic algorithm for the automatic configuration of algorithms[END_REF].

Among such automated algorithm configuration (AAC) tools, we cite GGA [START_REF] Ansótegui | A gender-based genetic algorithm for the automatic configuration of algorithms[END_REF], ParamILS [START_REF] Hutter | Paramils: an automatic algorithm configuration framework[END_REF], SPOT [START_REF] Bartz-Beielstein | SPOT: A toolbox for interactive and automatic tuning in the R environment[END_REF] and irace [START_REF] López-Ibáñez | The irace package: Iterated racing for automatic algorithm configuration[END_REF]. These methods have been successfully applied to (pre-tuned) state-of-the-art solvers of various problem domains, such as mixed integer programming [START_REF] Hutter | Automated configuration of mixed integer programming solvers[END_REF], AI planning [START_REF] Fawcett | Fd-autotune: Domain-specific configuration using fast downward[END_REF], machine learning [START_REF] Snoek | Practical bayesian optimization of machine learning algorithms[END_REF], or propositional satisfiability solving [START_REF] Hutter | The configurable SAT solver challenge (CSSC)[END_REF]. Figure 1 illustrates the abstract standard architecture adopted by these tools. However, the outcomes of these tools are static (or feature-independent), which means an algorithm configuration derived by any of these tools is not changed depending on an instance of a target optimisation problem. This leads to a significant issue as theoretical and empirical studies on various algorithms and problems have shown that parameters of an algorithm are highly dependent on features of a specific instance of a target problem [START_REF] Doerr | Fast genetic algorithms[END_REF] such as the problem size [START_REF] Böttcher | Optimal fixed and adaptive mutation rates for the LeadingOnes problem[END_REF][START_REF] Witt | Tight bounds on the optimization time of a randomized search heuristic on linear functions[END_REF].

A possible solution to this issue is to cluster problem instances into multiple sub-groups by their size (and other potential features), then use curve fitting to map features to parameters [START_REF] Mascia | Tuning algorithms for tackling large instances: an experimental protocol[END_REF][START_REF] El Yafrani | Efficiently solving the traveling thief problem using hill climbing and simulated annealing[END_REF]. A similar approach is also found in [START_REF] Liefooghe | Towards landscapeaware automatic algorithm configuration: preliminary experiments on neutral and rugged landscapes[END_REF] that first partitions problem instances based the values of their landscape features and selects an appropriate configuration of a new problem instance based on its closeness to the partitions. However, the former approach does not scale well to multiple features and parameters, and the latter faces over-fitting issues due to the nature of the partitioning approach, making it difficult to assign an unseen instance to a specific group. Some works have incorporated problem features in the parameter tuning process. SMAC [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF] and PIAC [START_REF] Leyton-Brown | Learning the empirical hardness of optimization problems: The case of combinatorial auctions[END_REF] are examples of model-based tools that consider instance features to define parameter values by applying machine learning techniques to build the model. However, an issue of these approaches is the low explainability of the outcome. For instance, while machine learning techniques such as random forest and neural networks can be used to map the parameters to problem features with a high accuracy, they are considered as black-boxes, i.e., the outcome is virtually impossible to understand or interpret. Explainability is an important concept, as not only it allows us to understand the relationships between input and output [START_REF] Rai | Explainable ai: from black box to glass box[END_REF], but in the context of parameter tuning, it can provide an outcome that can be used to inspire fundamental research [START_REF] Friedrich | Heavy-tailed mutation operators in single-objective combinatorial optimization[END_REF][START_REF] Friedrich | Escaping large deceptive basins of attraction with heavy-tailed mutation operators[END_REF].

To tackle these issues, we propose an offline algorithm tuning approach that extracts relationships between problem features and algorithm parameters using a genetic programming algorithm framework. We will refer to this approach as MATE, which stands for Model-based Algorithm Tuning Engine. The main contributions in this work are as follows:

1. We formulate the model-based parameter tuning problem as a symbolic regression problem, where knowledge about the problem is taken into account in the form of problem features; 2. We implement an efficient Genetic Programming (GP) algorithm that configures parameters in terms of problem features; and 3. In our empirical investigation, we rediscover asymptotically-correct theoretical results for two algorithms (1+1-EA and RLS) and four problems (OneMax, LeadingOnes, BinValue, and Jump). In these experiments, MATE shows its potential in algorithm parameter configuration to produce models based on instance features.

Background

Several methods have tried to tackle the dependence between the problem features and the algorithm parameters. The Per Instance Algorithm Configuration (PIAC) [START_REF] Leyton-Brown | Learning the empirical hardness of optimization problems: The case of combinatorial auctions[END_REF], for example, can learn a mapping between features and best parameter configuration, building an Empirical Performance Model (EPM) that predicts the performance of the algorithm for sample (instance, algorithm/configuration) pairs. PIAC methodology has been applied to several combinatorial problems [START_REF] Hutter | Performance prediction and automated tuning of randomized and parametric algorithms[END_REF][START_REF] Xu | SATzilla: portfolio-based algorithm selection for SAT[END_REF][START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF] and continuous domains [START_REF] Belkhir | Per Instance Algorithm Configuration of CMA-ES with Limited Budget[END_REF]. Sequential Model-based Algorithm Configuration (SMAC) [START_REF] Hutter | Sequential model-based optimization for general algorithm configuration[END_REF] is also an automated algorithm configuration tool which considers a model, usually a random forest, to design the relationship between a performance metric (e.g. the algorithm runtime) and algorithm parameter values. SMAC can also include problem features within the tuning process as a subset of input variables.

Table 1 presents a summary for some state-of-the-art methods including the approach proposed in this paper. The term 'feature-independent' means that the corresponding approach does not consider instance features. 'Model-based' approaches use a trained model (e.g. machine learning, regression, etc.) to design parameter configurations. Model-free approaches generally rely on an experimental design methodology or optimisation method to find parameter settings of an algorithm that optimise a cost metric on a given instance set.

The main differences between MATE and the other related approaches are:

1. A transparent machine learning method (GP) is utilised to enable humanreadable configurations (in contrast to, e.g., random forests, neural networks, etc.). 2. The training phase is done on one specific algorithm and one specific problem in our approach -the model is less instance-focused but more problem-domain focused by abstracting via the use of features. For example, the AAC experiments behind [START_REF] Friedrich | Heavy-tailed mutation operators in single-objective combinatorial optimization[END_REF][START_REF] Friedrich | Escaping large deceptive basins of attraction with heavy-tailed mutation operators[END_REF] have guided the creation of new heavy-tailed mutation operators that were beating the state-of-the-art. Similarly, the AAC and PIAC experiments in [START_REF] Treude | Predicting good configurations for github and stack overflow topic models[END_REF] showed model dependencies on easily-deducible instance features.

Lastly, our present paper is much aligned with the recently founded research field "Data Mining Algorithms Using/Used-by Optimisers (DUO)" [START_REF] Agrawal | Better software analytics via "duo": Data mining algorithms using/used-by optimizers[END_REF]. There, data miners can generate models explored by optimisers; and optimisers can adjust the control parameters of a data miner.

The MATE Framework

Problem Formulation and Notation

Let us denote an optimisation problem by B whose instances are characterised by the problem-specific features F = {f 1 , . . . , f M }. A target algorithm A with its parameters P = {p 1 , . . . , p N } is given to address the problem B. A set of instances I = {i 1 , . . . , i L } of the problem B and a L × M matrix V, whose element value v i,j represents the jth feature value of the ith problem instance, are given.

Under this setting, we define the model-based parameter tuning problem as the problem of deriving a list of mappings M = {m 1 , . . . , m N } where each mapping m j : R M → R, which we will refer to as a parameter expression, returns a value for the parameter p j given feature values of an instance of the problem B. Specifically, the objective of the problem is to find a parameter expression set M * , such that the performance of the algorithm A across all the problem instances in I is optimised.

Architecture Overview

In this section, we introduce our approach for parameter tuning based on the problem features. Figure 2 illustrates the architecture of the MATE tuning engine. In contrast to static methods, we consider the features of the problem. These feature are to be used in the training phase in addition to the instances, the target algorithm and the parameter specifications. Once the training is finished, the model can be used on unseen instances to find the parameters of the algorithm in terms of the problem feature values of the instance. For example, a desired outcome of applying the MATE framework can be:

-Mutation probability of an evolutionary algorithm in terms of the problem size; -Perturbation strength in an iterated local search algorithm in terms of the ruggedness of the instance and the problem size; and -Population size of an evolutionary algorithm in terms of the problem size.

Note that all the examples include the problem size as a problem feature. In both theory and practice, the problem size is among the most important problem features, and it is usually known prior to the optimisation, without any need for a pre-processing step. More importantly, an extensive number of theoretical studies showed that the optimal choice of parameters is usually expressed in terms of the problem size (see, e.g. [START_REF] Böttcher | Optimal fixed and adaptive mutation rates for the LeadingOnes problem[END_REF][START_REF] Doerr | Fast genetic algorithms[END_REF][START_REF] Witt | Tight bounds on the optimization time of a randomized search heuristic on linear functions[END_REF]).

The Tuning Algorithm

We use a tree-based Genetic Programming system as the tuning algorithm. It starts with a random population of trees, where each tree represents a potential parameter expression. Without loss of generality, we assume that the target problem is always a maximisation problem 5 .

The Score Function and Bias Reduction The score function is expressed as the weighted sum of the obtained objective values on each instance in the training set I. Using the notations previously introduced, the score is defined in Equation (1):

S(t) = 1 L Σ i∈I z A (m 1 (v i,1 , . . . , v i,M ), . . . , m N (v i,1 , . . . , v i,M ), i) R i (1)
where:

-S(.) is the GP score function, z A (ϕ 1 , . . . , ϕ N , i) is a function measuring the goodness of applying the algorithm A with the parameter values ϕ 1 , . . . , ϕ N to instance i, -R i is the best known objective value for instance i.

The weights are used as a form of normalisation to reduce the bias some instances might induce. A solution to address this issue would be to use the optimal value or a tight upper bound. However, since we assume the such values are unknown (the problem itself can be unknown), we use the best known objective value (R i ) as a reference instead. In order to always ensure that score is well contained, the reference values are constantly updated whenever possible during the tuning process.

Replacement Strategy -Statistical Significance and Bloat Control

As the target algorithm can be stochastic, it is mandatory to perform multiple runs to ensure statistical significance (refer to Table 3). Thus, the replacement of trees is done based on the Wilcoxon rank-sum test.

Another aspect to take into account during the replacement process is bloat control. In our implementation, we use a simple bloat minimisation method based on the size of tree (number of nodes).

Given a newly generated tree (Y ), we compare it against each tree (X) in the current population starting from the ones with the lowest scores using the following rules:

-If Y is deemed to be significantly better than X (using the Wilcoxon test).

then we replace X with Y irrespective of the sizes. -If there is no statistical significance between X and Y , but Y has a smaller size than X, then we replace X with Y . -Otherwise, we do not perform the replacement.

Experimental Setting

To evaluate our framework, we consider two target algorithms, the (1+1) EA(µ) and RLS(k). The (1+1) EA(µ) is a simple hill-climber which uses standard bit mutation with mutation rate µ. RLS(k) differs from the (1+1) EA(µ) only in that it uses the mutation operator that always flips k uniformly chosen, pairwise different bits. That is, the mutation strength k is deterministic in RLS, whereas it is binomially distributed in case of the (1+1) EA(µ), Bin(n, µ), where n is the number of bits.

We use MATE to configure the two algorithms for the four different problems with different time budgets as summarised in Table 2. In the table, the features of the problems used to tune the algorithm parameters and the different feature values chosen to generate problem instances of the problems are also presented. These problems have been chosen because they are among the best studied benchmark problems in the theory of evolutionary algorithms [START_REF] Doerr | Theory of Evolutionary Computation-Recent Developments in Discrete Optimization[END_REF]. The details of our GP implementation for the experiments are presented in Table 3. Based on Table 3 and the set of features, our GP method uses a minimalistic set of 6 terminals at most: m, n and {1, 2, -1, -2}.

It is worth noting that we are focusing in this paper on tuning algorithms with a single parameter. This is done to deliver a first prototype that is validated on algorithms and problems extensively studied by the EA theory community. An extension to tuning several algorithm parameters forms an important direction for future work.

For example, given a budget of (1 + o(1))en ln(n), it is known that the (1+1)EA(1/n) optimises the OneMax function as well as any other linear functions with a decent probability. It is also known that the 1/n is asymptotically optimal [START_REF] Lengler | Fixed budget performance of the (1+1) EA on linear functions[END_REF]. Note, though, that such fixed-budget results are still very sparse [START_REF] Jansen | Analysing stochastic search heuristics operating on a fixed budget[END_REF], since the theory of EA community largely focuses on expected optimisation times. Since these can nevertheless give some insight into the optimal parameter settings, we note the following:

-OneMax and BinValue: the (1+1)EA(1/n) optimises every linear function in expected time en ln(n), and no parameter configuration has smaller expected running time, apart from possible lower order terms [START_REF] Witt | Tight bounds on the optimization time of a randomized search heuristic on linear functions[END_REF]. For RLS, it is not difficult to see that k = 1 yields an expected optimisation time of (1 + o(1))n ln(n), and that this is the optimal (static) mutation strength; -LeadingOnes: on average, RLS(1) needs n 2 /2 steps to optimise LeadingOnes.

This choise also minimises the expected optimisation time. For the (1+1) EA, µ ≈ 1.59/n minimises the expected optimisation time, which is around 0.77n 2 for this setting [START_REF] Böttcher | Optimal fixed and adaptive mutation rates for the LeadingOnes problem[END_REF]. The standard mutation rate µ = 1/n requires 0.86n 2 evaluations, on average, to locate the optimum, of the LeadingOnes function.

For LeadingOnes, it is known that the optimal parameter setting drastically depends on the available budget. This can be inferred from the proofs in [START_REF] Böttcher | Optimal fixed and adaptive mutation rates for the LeadingOnes problem[END_REF][START_REF] Doerr | Analyzing randomized search heuristics via stochastic domination[END_REF]; and -Jump: mutation rate m/n minimises the expected optimisation time of the (1+1) EA on Jump(m, n), which is nevertheless Θ((e/m) m n m ) [START_REF] Doerr | Fast genetic algorithms[END_REF].

Performance Analysis

Training Phase The experimental study is conducted by running MATE ten times on each algorithm, problem and budget combination (refer to Table 4 for the list of budgets). This results in an elite population of 20 individuals for each setting, from which we select the top 5 expressions in terms of the score. These results are then merged and the 3 most frequent expressions are selected. For instance, the expression 2/n for OneMax with 0.5enln(n) appears 92 times over the 200 individuals (population size (20) × runs ( 10)).

In the current implementation, expression types (integers and non-integers) are not taken into account during the evolution. Therefore, the resulting expressions are converted into integers in the case of RLS by merging all real numbers r using r (e.g. k = 3/2 will be replaced by k = 1). On the other hand, expressions are simplified for EA by eliminating additive constants (e.g. µ = 1/(n + 1) is replaced by µ = 1/n).

Evaluation Phase I To assess the performance of MATE, we evaluate for each problem-budget combination each of the top 3 most frequent expressions, by running them 100 independent times on each training dimension. We then normalise the outputs as in Equation [START_REF] Agrawal | Better software analytics via "duo": Data mining algorithms using/used-by optimizers[END_REF]. The results are shown in the box plots in Table 4.

Comparison amongst the top 3 configurations. When comparing the top 3 ranked configurations, we observe the following from Table 4 while we compare medians.

-OneMax: For (1+1) EA, µ = 1/n, which ranked second for budgets 0.5en ln n and en ln n and first for budget 2en ln n performs better than µ = 1/2 * n; while for RLS, the expression k = 1 appears at least on 94%, providing the best results; -BinValue: µ = 1/n represents 18% on en ln n for (1+1) EA experiments, and a similar performance with µ = 2/n and µ = 3/n; while on 0.5en ln n case the µ = 1/n expression provides better results than µ = 1/2 and µ = 1/3; on the same way the expression k = 1 corresponds to 60% of the cases on RLS with the budget of 2n ln n with a better performance than k = 2 and k = n; -LeadingOnes: µ = 1/n is the most frequent expression among all considered budgets on (1+1) EA and µ = 2/n presents the best performance amongst the top 3 expressions for all budget cases; k = 1 represents 88% on RLS cases with 0.75n 2 iterations and performs better than k = 2 and k = 3 for both considered budgets. -Jump: µ = 2/n and µ = m/n present similar results for both budget cases; µ = 1/n appears on 36% and 68% of the cases on (1+1) EA on the considered budgets respectively, and performs worse than the other two µ configurations; for RLS experiments k = m is the most frequent expression and performs better than k = 2 * m and k = 3.

Comparison of top 3 configurations against other parameter settings. For a fair assessment of our results, we add to this comparison some expressions that were not ranked in the top 3. These are µ = i/n with i ∈ {1, 3/2, 2, 5/2, 3, 4} for (1+1) EA(µ) for OneMax and LeadingOnes. For readability purposes, the top 3 expressions are complemented with 3 of these additional expressions in the same order they are shown. We can observe in Table 4 that these additional expressions present low frequencies, µ = 3/n being the highest case with 12% with the budget en ln n, while expressions µ = 3/(2n) and µ = 5/(2n) are the lowest cases among the considered budgets. Note that the frequencies do not necessarily sum up to 100% as other expressions not reported here might have occurred.

Comparison with theoretical results. As we have mentioned in the beginning of this section, one should be careful when comparing theoretical results that have been derived either in terms of running time or in terms of asymptotic convergence analysis, as typically done in runtime analysis. It is well known that optimal parameter settings for concrete (typically, comparatively small) dimensions can be different from the asymptotically optimal ones [START_REF] Buskulic | Maximizing drift is not optimal for solving onemax[END_REF][START_REF] Chicano | Fitness probability distribution of bit-flip mutation[END_REF]. We nevertheless see that the configurations that minimise the expected running times (again, in the classical, asymptotic sense) also show up in the top 3 ranked configurations. In Table 4, we highlight the asymptotically optimal best possible running time

2/n [46%] 1/n [32%] 1/(2*n) [10%] 3/n [4%] 5/(2*n) [2%] 3/(2*n) [0%] n ln(n)* 0.85 0.90 0.95 1.00 1 [98%] 3 [2%] 2 [0%]
en ln(n)* 0.95 0.96 0.97 0.98 0.99 1.00

2/n [46%] 1/n [26%] 1/(2*n) [14%] 3/n [12%] 5/(2*n) [2%] 3/(2*n) [0%]
2n ln(n)** 0.90 0.92 0.94 0.96 0.98 1.00

1 [94%] 3 [4%] 2 [2%]
2en ln(n)** 0.975 0.980 0.985 0.990 0.995 1.000 2n ln(n)** 0.980 0.985 0.990 0.995 1.000

1/n [44%] 2/n [26%] 1/(2*n) [12%] 3/n [8%] 3/(2*n) [8%] 5/(2*n) [0%] BinValue 0.
1 [60%] 2 [32%] n [6%]
LeadingOnes 0.5n 2 0.6 0.7 0.8 0.9

1/n [52%] 2/n [28%] 4/n [20%] 3/n [0%] 5/(2*n) [0%] 3/(2*n) [0%] 0.5n 2 * 0.7 0.8 0.9 1.0 2 [70%] 1 [26%] 3 [4%]
0.8n 2 0.70 0.75 0.80 0.85 0.90 0.95 1.00

1/n [62%] 2/n [18%] [16%] 4/n [2%] 5/(2*n) [0%] 3/(2*n) [0%]
0.75n 2 ** 0.75 0.80 0.85 0.90 0.95 1.00

1 [88%] 2 [12%] 3 [0%]
0.9n 2 ** 0.75 0.80 0.85 0.90 0.95 1.00

1/n [48%] 2/n [28%] 3/n [18%] 5/(2*n) [4%] 4/n [0%] 3/(2*n) [0%]
Jump n m 0.875 0.900 0.925 0.950 0.975 1.000

1/n [36%] 2/n [32%] m/n [12%]
n m 0.850 0.875 0.900 0.925 0.950 0.975 1.000

m [34%] 2*m [24%] 3 [20%]
en m ** 0.92 0.94 0.96 0.98 1.00

1/n [68%] 2/n [22%] m/n [6%]
2n m 0.850 0.875 0.900 0.925 0.950 0.975 1.000

m [42%] 2*m [20%]
3 [18%] † The y-axis show the best found expressions with its frequency between square brackets, and the x-axis represents the normalised fitness.

by an asterisk *. Budgets exceeding this bound are marked by two asterisks **. As for the individual problems, we note the following:

-OneMax: It is interesting to note here that the performance is not monotonic in k, i.e., k = 2 performs worse than k = 1 and k = 3. This is caused by a phenomenon described in [START_REF] Doerr | Optimal parameter choices via precise black-box analysis[END_REF]Section 4.3.1], which states that, regardless of the starting point, the expected progress is always maximised by an uneven mutation strength. MATE correctly identifies this and suggests uneven mutation strengths in almost all cases. -BinValue: We observe that it is very difficult here to distinguish the performance of the different configurations. This is in the nature of BinValues, as setting the first bit correctly already ensures 50% of the optimal fitness values. We very drastically see this effect in the recommendation to use k = n for the RLS cases. With this configuration, the algorithm evaluates only two points: the random initial point x and its pairwise complement x, regardless of the budget. As can be seen in Table 4, the performance of this simple strategy is quite efficient, and hard to beat -LeadingOnes: As mentioned earlier, for the (1+1) EA, the optimal mutation rate in terms of minimising the expected running time is around µ = 1.59/n. We see that µ = 3/(2n), which did not show in the top 3 ranked configurations performs better than any of the suggestions by MATE. -Jump: as discussed, mutation rate µ = m/n minimises the expected optimisation time. MATE recognises it as a good configuration in some of the runs. However, we see that µ = 2/n, which equals µ = m/n for 5 out of our 12 training sets, shows comparable performance, and in the en m budget case even slightly better performance.

Evaluation Phase II To properly assess the performance of MATE, we conducted experiments for OneMax and LeadingOnes instances of larger sizes that were not considered in the training phase. The goal of this experiment is to empirically demonstrate that our approach generalises well for large and unseen instances. These results are presented in Table 5 where 100 runs were performed for OneMax with n ∈ {1000, 2000, 5000} and LeadingOnes with n ∈ {750, 1000}.

We can observe the following:

-There is less overlap amongst the confidence intervals especially for smaller budgets, which means there is a higher level of separability amongst the performances of the different expressions. -By comparing these results with the ones from Table 4, we can observe that the results of the top 3 expressions on large instances are statistically better in the majority of cases. -OneMax: For (1+1) EA, in contrast to the results in Table 4 where µ = 1/n and µ = 3/(2n) show a similar performance, here µ = 1/n performs better than the other expressions. For RLS, the best performing expression is k = 1, which was ranked first. -LeadingOnes: For (1+1) EA the best expressions are µ = 2/n, which was ranked second, and µ = 3/(2n), which was not ranked among the top 3 expressions. For RLS, k = 1, ranked first and second, is the best performing expression.

Comparative study

Herein, we compare the performance of MATE with irace and SMAC. The goal is to investigate the sensitivity of the obtained parameters on unseen instances. Table 6 shows that MATE significantly outperforms irace and SMAC for (1+1) EA. On the other hand, the three methods show a similar performance on RLS. This is due to the fact that the parameter µ in (1+1) EA is highly sensitive to the problem feature n. In contrast, the parameter k in RLS is independent from n and its best value (k = 1) was identified by the three methods for both OneMax and LeadingOnes.

Conclusions and Future Directions

With this article, we have presented MATE as a model-based algorithm tuning engine: its human-readable models map instance features to algorithm parame-ters. Our experiments showed that MATE can find known asymptotic relationships between the feature values and algorithm parameters. We also compared the performance of MATE with iRace and SMAC investigating the sensitivity of the obtained parameters on unseen instances of larger size. With its scalable models, MATE performed best. It is worth noting that MATE can be a useful guideline tool for theory researchers due to its white-box nature, similarly to how results in [START_REF] Doerr | Simple on-the-fly parameter selection mechanisms for two classical discrete black-box optimization benchmark problems[END_REF] inspired the analysis of a generalised one-fifth success rule in [START_REF] Doerr | Self-adjusting mutation rates with provably optimal success rules[END_REF]. But MATE can also be extended to be used as a practical toolbox for feature-based algorithm configuration.

In the future, we intend to explore, among other, the following three avenues. First, the design of MATE itself will be subject to extensions, e.g. to better handle performance differences between instances via ranks or racing. Second, while our proof-of-concept study here was motivated by theoretical insights, we will investigate more realistic problems for which instance features are readily available, such as the travelling salesperson problem and the assignment problem. Third, we will investigate approaches to extend MATE to handle multiple parameters to demonstrate its ability to tune more sophisticated algorithms.

Fig. 1 :

 1 Fig. 1: Standard architecture of tuning frameworks.

Fig. 2 :

 2 Fig. 2: Architecture of the proposed MATE framework

Table 1 :

 1 Summary of the state-of-the-art related works

	Approach Name	Algorithm	Characteristics		Ref.
	GGA	Genetic Algorithm	Feature-independent, model-free [2]
	ParamILS	Iterated Local Search	Feature-independent, model-free [23]
	irace	Racing procedure	Feature-independent, model-free [30]
	SPOT	classical regression, tree-based, ran-dom forest and Gaussian process	Feature-independent, based	model-	[3]
	PIAC	Regression methods	Feature-dependent, model-based [28]
	SMAC	Random Forest	Feature-dependent, model-based [22]
	MATE	Genetic Programming	Feature-dependent, based, explainable	model-

Table 2 :

 2 Summary of problems

	Problem	Features	Training set
	OneMax(n)	n: number of bits	n = 10, 20, 50, 100, 200, 500
	BinValue(n)	n: number of bits	n = 10, 20, 50, 100, 200, 500
	LeadingOnes(n)	n: number of bits	n = 10, 20, 50, 100, 200, 500
			(m, n) = (2, 10), (3, 10), (4, 10), (5, 10),
		m: width of region with bad	(2, 20), (3, 20), (4, 20),
	Jump(m, n)	fitness values	(2, 50), (3, 50),
		n: number of bits	(2, 100), (3, 100),
			(2, 200)

Table 3 :

 3 MATE setup

	Attribute/Parameter	Value/Content
	Terminals	{1, 2, -1, -2} F
	Functions	Arithmetic operators
	Number of GP generations	100
	Population size	20
	Tournament size	5
	Replacement rate	< 75%
	Initialisation	grow (50%) and full (50%) methods
	Mutation operator	random mutations
	Mutation probability	0.2
	Crossover operator	sub-tree gluing
	Crossover rate	80%
	Number of independent runs of target algorithm	10
	p-value for the Wilcoxon ranksum test	0.02
	4 Computational Study

Table 4 :

 4 Results for 20 settings.

		1+1-EA			RLS
		Budget	Result		Budget	Result
	OneMax	0.5en ln(n)	0.85	0.90	0.95

Table 5 :

 5 Results for larger OneMax and LeadingOnes instances

		1+1-EA							RLS
		Budget	Result						Budget	Result
	OneMax	0.5en ln(n) en ln(n)*	2/n [46%] 1/n [32%] 1/(2*n) [10%] 3/n [4%] 5/(2*n) [2%] 3/(2*n) [0%] 2/n [46%] 1/n [26%] 1/(2*n) [14%] 3/n [12%] 5/(2*n) [2%] 3/(2*n) [0%]	0.975 0.980 0.985 0.990 0.995 0.994 0.995 0.996 0.997 0.998 0.999 1.000	n ln(n)* 2n ln(n)**	1 [98%] 3 [2%] 2 [0%] 1 [94%] 3 [4%] 2 [2%]	0.94 0.97	0.96 0.98	0.98 0.99	1.00 1.00
		2en ln(n)**	1/n [44%] 2/n [26%] 1/(2*n) [12%] 3/(2*n) [8%] 3/n [8%] 5/(2*n) [0%]	0.994	0.996		0.998	1.000
	LeadingOnes	0.5n 2 0.8n 2 **	1/n [52%] 2/n [28%] 4/n [20%] 3/n [0%] 5/(2*n) [0%] 3/(2*n) [0%] 1/n [62%] 2/n [18%] 3/n [16%] 4/n [2%] 3/(2*n) [0%] 5/(2*n) [0%]	0.65 0.70 0.75 0.80 0.85 0.80 0.85 0.90 0.95 1.00	0.5n 2 * 0.75n 2 **	2 [70%] 1 [26%] 3 [4%] 0.70 0.75 0.80 0.85 0.90 0.95 1.00 0.80 0.85 0.90 0.95 1.00 1 [88%] 2 [12%] 3 [0%]
		0.9n 2 **	1/n [48%] 2/n [28%] 3/n [18%] 5/(2*n) [4%] 4/n [0%] 3/(2*n) [0%]	0.80	0.85	0.90	0.95	1.00

Table 6 :

 6 For a fair comparison, we run irace and SMAC with 2000 maximum experiments (which we believe is equivalent to the 100 GP generations with a population size of 20 individuals in MATE) considering the training instances presented in Table2. We report the best elite parameter values returned by irace (2 candidates), SMAC (1 candidate) and MATE (most frequent expressions) in the columns µ and k in Table6, while the score (Eq. 1) is shown in column Score with the standard deviation as a subscript. These parameter values are then applied over 100 runs performed for OneMax with n ∈ {1000, 2000, 5000} and LeadingOnes with n ∈ {750, 1000}. Results for MATE, irace and SMAC for OneMax and LeadingOnes instances. .025 0.430 0.030.002 0.024 0.290.007 0.5n2 2 0.860.014 1 0.980.026 1 0.980.02

					1+1-EA			RLS
				MATE		irace		SMAC	MATE	irace	SMAC
		Budget µ Score	µ	Score	µ	Score	Budget k Score	k Score	k Score
		enln(n) 2	2 n	0.990.001 0.258 0.570.002 0.009 0.80.003 nln(n) 1 10	1 10	1 10
			1 n	0.990.001 0.216 0.580.002			3 0.960.002
			1 2n 0.980.002					2 0.940.003
	OneMax	enln(n) 1 n 2 n 1 2n 10 10 10	0.009 0.820.002 0.016 0.760.003 2nln(n) 1 10 0.013 0.790.003 3 0.980.001 2 0.970.002	1 10	1 10
		2enln(n) 2 n	10	0.594 0.540.002 0.008 0.860.002
			1 n	10	0.589 0.540.002
			1 2n 10
		0.5n 2	1 n 0.702 n 0.80.021 0.409 0.030.002			1 0.970.027 5 0.610.01
	LeadingOnes	0.8n 2 0.9n 2	4 n 1 n 2 n 3 n 1 n	0.710.015 0.950.023 0.255 0.050.002 0.005 0.830.017 0.75n 2 1 10 3 0.750.013 0.990.012 0.258 0.050.002 2 0.950.009 0.910.018 3 0.820.01 0.990.011 0.158 0.070.003 0.006 0.750.013	1 10	1 10
			2 n	10	0.153 0.070.006
			3 n	0.950.014

The current MATE implementation is publicly available at https://gitlab.com/ yafrani/mate

Acknowledgements

M. Martins acknowledges CNPq (Brazil Government). M. Wagner acknowledges the ARC Discovery Early Career Researcher Award DE160100850. C. Doerr acknowledges support from the Paris Ile-de-France Region. Experiments were performed on the AAU's CLAUDIA compute cloud platform.


Maxim Buzdalov

Carola Doerr

Optimal Static Mutation Strength Distributions for the (1 + 𝜆) Evolutionary Algorithm on OneMax

Most evolutionary algorithms have parameters, which allow a great flexibility in controlling their behavior and adapting them to new problems. To achieve the best performance, it is often needed to control some of the parameters during optimization, which gave rise to various parameter control methods. In recent works, however, similar advantages have been shown, and even proven, for sampling parameter values from certain, often heavy-tailed, fixed distributions. This produced a family of algorithms currently known as "fast evolution strategies" and "fast genetic algorithms".

However, only little is known so far about the influence of these distributions on the performance of evolutionary algorithms, and about the relationships between (dynamic) parameter control and (static) parameter sampling. We contribute to the body of knowledge by presenting, for the first time, an algorithm that computes the optimal static distributions, which describe the mutation operator used in the well-known simple (1 + 𝜆) evolutionary algorithm on a classic benchmark problem OneMax. We show that, for large enough population sizes, such optimal distributions may be surprisingly complicated and counter-intuitive. We investigate certain properties of these distributions, and also evaluate the performance regrets of the (1 + 𝜆) evolutionary algorithm using commonly used mutation distributions.

INTRODUCTION

Many evolutionary algorithms are composed of operators that are used within in several algorithms. In the context of bit string representations, a very popular example of such an operator is standard bit mutation (SBM), the variation operator that takes as input a point 𝑥 (the "parent"), modifies it by changing the entry in each position with some positive probability 𝑝 (independently of all other decisions), and outputs the so-obtained "offspring" 𝑦. SBM is a global search operator, since the probability that it generates a particular point 𝑦 is positive, regardless of the input 𝑥. Another common mutation operator is the flip 1 operator, which creates the offspring by changing the entry in exactly one uniformly selected position. This search operator is a local one. Both mutation operators are unbiased in the sense proposed by Lehre and Witt [START_REF] Per | Black-Box Search by Unbiased Variation[END_REF], i.e., they treat all positions and all possible values identically.

By a characterization derived in [START_REF] Doerr | Optimal parameter choices via precise black-box analysis[END_REF], unary unbiased operators are exactly the ones that can be described by a distribution 𝐷 over the integers {0, 1, . . . , 𝑛}. To apply them, one first samples a mutation strength ℓ from this distribution and then changes the input by applying the flip ℓ operator, which flips the entry in ℓ uniformly chosen, pairwise different positions. That SBM is unbiased now easily follows from the observation that it can be exactly characterized by the binomial distribution Bin(𝑛, 𝑝) with 𝑛 trials (one for each position) and success probability 𝑝 (success=bit flip). Likewise, flip 1 is the operator associated to the 1-point distribution assigning all probability mass to mutation strength 1.

Our Contribution: We want to analyze in this work how common mutation operators, such as the ones mentioned above or the fast mutation operator suggested in [START_REF] Doerr | Fast genetic algorithms[END_REF], compare to an optimal one. To analyze this question, we introduce the (1 + 𝜆) UUSD-EA, the family of all (1 +𝜆)-type mutation-only algorithms whose mutation operator can be defined via a unary unbiased satic distribution (which may hence depend on the problem dimension 𝑛 and the offspring population size 𝜆, but which may not change during the run). This family comprises all (1 + 𝜆) evolutionary algorithms (EAs), their Randomized Local Search (RLS) counterparts (which use the flip 1 operator instead of SBM), the fastGAs, the normalized EAs [START_REF] Ye | Interpolating Local and Global Search by Controlling the Variance of Standard Bit Mutation[END_REF], etc.

We numerically compute for various combinations of 𝑛 and 𝜆 the minimal expected runtime that can be obtained by any (1+𝜆) UUSD-EA on OneMax in dimension 𝑛, and we compare these runtimes to that of the classically studies (1 + 𝜆)-type algorithms mentioned above. Our lower bounds are constructive, in that we also derive the distributions associated to these optimal (1 + 𝜆) UUSD-EAs. This allows us to study the properties of these optimal distributions. Approach: Since the optimal distributions cannot be determined by exact analytical approaches, we apply the separable CMA-ES [START_REF] Ros | A simple modification in CMA-ES achieving linear time and space complexity[END_REF] to identify them. The CMA-ES shows very good performance, and provides us with distributions that are optimal up to the last few digits in the available machine precision. This is also one of the cases where the separable version of CMA-ES is not only computationally faster, but also yields results of same or better quality than the default CMA-ES implementations. Main result: Our numerical results show, among other things, that the flip 1 operator is optimal when 𝜆 is small, whereas the conditional binomial distribution Bin >0 (𝑛, 1/2) is optimal when 𝜆 is very large -with this distribution, the (1 +𝜆) UUSD-EA performs a uniform random search. The optimal distributions in the middle regime show a rather complex behavior, for which we can identify a few patterns, but for which we also observe a few phenomena that may look counter-intuitive at first glance. For instance, the first mutation strength different from one-bit flip that gets nonzero probability when 𝜆 grows appears to be either 𝑛 or 𝑛 -1.

Relationship to black-box complexity and to parameter control: Our work can be seen as a continuation of black-box complexity theory for 𝜆-parallel [START_REF] Badkobeh | Unbiased Black-Box Complexity of Parallel Search[END_REF] elitist [START_REF] Doerr | Introducing Elitist Black-Box Models: When Does Elitist Behavior Weaken the Performance of Evolutionary Algorithms?[END_REF] unary unbiased [START_REF] Per | Black-Box Search by Unbiased Variation[END_REF] black-box algorithms with static configuration. A key advantage of lower bounds such as ours is that it allows to rigorously quantify the impact of individual decisions. In our case, the driving motivation behind our analyses is a rigorous quantification of the difference between static and dynamic algorithm configurations. Put differently, we aim at quantifying the gap between algorithms using parameter control and such that do not. In contrast to classical runtime and black-box complexity results, our work focuses on an exact numeric evaluation of this gap for concrete problem dimensions.

Impact: While the results of our concrete analysis may mostly appeal to theoreticians, our work invites to take a different view on mutation operators by defining them via distributions over the possible mutation strengths. This alternative view makes it substantially easier to generalize concepts such as SBM, flip 1 , etc., an advantage that can be leveraged, for example, for designing smooth convergence from global to local search behavior. But the design principle can also lead to performance gains in the static setting. A first strong result in this context is the fastGA [START_REF] Doerr | Fast genetic algorithms[END_REF], which has become the new state-of-the-art in several applications [START_REF] Meunier | Black-Box Optimization Revisited: Improving Algorithm Selection Wizards through Massive Benchmarking[END_REF]. Our work indicates that there is quite some untapped potential in the design of new mutation operators, and we hope that our work inspires new work in this direction.

Based on the findings of our work, we formulate to conjectures:

Conjecture 1: We conjecture that for each 𝑛 there exists a threshold 𝜆 1 (𝑛) such that for all 𝜆 ≤ 𝜆 1 (𝑛) the 1-point distribution is optimal. Our guess on the particular dependency is 𝜆 1 (𝑛) = Θ(log 𝑛). However, we currently do not have a formal proof for this.

Conjecture 2:

We conjecture that for each 𝑛 and each arbitrarily small 𝜀 there exists another threshold 𝜆 𝑏 (𝑛, 𝜀) such that, for all 𝜆 ≥ 𝜆 𝑏 (𝑛, 𝜀), the optimal distribution is closer than 𝜀, in any imaginable measure -such as the maximum of differences between probabilities across all mutation strengths -to the conditional binomial distribution Bin >0 (𝑛, 1/2). This is equivalent to stating that the uniform random search is arbitrarily close to being optimal for large enough 𝜆. Again, we do not yet have a formal proof.

Related Work:

Our study continues the recent works [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF] and [START_REF] Buzdalov | Optimal Mutation Rates for the (1 + 𝜆) EA on OneMax[END_REF], which provide optimal dynamic configurations for (1+1) and (1+𝜆)type algorithms, respectively. While their works are restricted to specific mutation operators (variants of SMB and the 𝑘-bit flips flip 𝑘 ), we study in this work in a generalization to arbitrary unary unbiased variation operators. In contrast to [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF][START_REF] Buzdalov | Optimal Mutation Rates for the (1 + 𝜆) EA on OneMax[END_REF] we focus on static configurations, with the idea to build a rigorous baseline against which we can compare dynamic parameter control methods.

For 𝜆 = 1, the work [START_REF] Doerr | Optimal parameter choices via precise black-box analysis[END_REF] quantifies the asymptotic advantage of the best unary unbiased algorithm with dynamic distributions against the best static one (RLS). To extend this work to (1 + 𝜆)-type algorithms, a rigorous bound on the best static case is needed -a baseline that we provide in this work for various combinations of 𝑛 and 𝜆, with the hope that the insights generated by our examples can be leveraged to rigorously prove certain characteristics of the optimal static unary unbiased operators.

For 𝜆 > 1, related works can be found in the context of the parallel black-box complexity model [START_REF] Badkobeh | Unbiased Black-Box Complexity of Parallel Search[END_REF][START_REF] Badkobeh | Black-box Complexity of Parallel Search with Distributed Populations[END_REF][START_REF] Per | Parallel Black-Box Complexity with Tail Bounds[END_REF] and for the (1 + 𝜆) EA [START_REF] Gießen | The Interplay of Population Size and Mutation Probability in the (1 + 𝜆) EA on OneMax[END_REF][START_REF] Gießen | Optimal Mutation Rates for the (1+𝜆) EA on OneMax Through Asymptotically Tight Drift Analysis[END_REF]. All these works, however, are either less interested in exact runtime bounds (and focus on asymptotic runtime guarantees instead) or they concern specific mutation operators only. For generalized mutation operators, concrete examples can be found in the mentioned works [START_REF] Doerr | Fast genetic algorithms[END_REF][START_REF] Ye | Interpolating Local and Global Search by Controlling the Variance of Standard Bit Mutation[END_REF]. We are not aware, however, of previous works explicitly studying optimal unary mutation operators.

Availability of Code and Data: All project source code and data are available for public use at [START_REF] Anonymized | Code and data repository for this paper[END_REF].

FROM MUTATION OPERATORS TO MUTATION STRENGTH DISTRIBUTIONS

We are concerned in this work with a generalizing view on unary unbiased variation operators, often referred to in the evolutionary computation context as "mutation". In a nutshell, a mutation operator takes as input a search point 𝑥 ∈ S, S denoting the search space, and creates from it an offspring 𝑦 ∈ S. More formally, a mutation operator is a family (𝐷 (𝑥)) 𝑥 ∈S of unary distributions over the search space S. When fed with an input 𝑥, a new search point 𝑦 is sampled from 𝐷 (𝑥). One of the most common mutation operators is standard bit mutation (SBM). In the context of pseudo-Boolean optimization (i.e., the maximization of a function 𝑓 : {0, 1} 𝑛 → R) -which is the setting that we assume for the remainder of the paper -SBM is often explained as follows: to obtain an offspring 𝑦 from 𝑥, we first create a copy of 𝑥 and then we decide for each bit position 𝑖 ∈ [𝑛] := {1, . . . , 𝑛} whether the entry shall be updated to 1 -𝑥 𝑖 ("bit flip") or whether the current entry is maintained.The bit flip decisions are made independently of each other. The probability 𝑝 ∈ (0, 1] to flip an entry is referred to as the mutation rate. SBM is a prime example of a unary unbiased mutation operator in the sense proposed by Lehre and Witt in [START_REF] Per | Black-Box Search by Unbiased Variation[END_REF]. By a characterization proven in [13, Lemma 1], this class subsumes all variation operators that are fully described by a distribution 𝐷 over the possible mutation strengths ℓ ∈ [0..𝑛] := [𝑛] ∪ {0}. When applied to a search point 𝑥, the operator first samples a mutation strength ℓ ∈ [0..𝑛] from its operator-specific distribution 𝐷 and then creates the offspring 𝑦 by flipping the entries in ℓ pairwise different, uniformly selected entries.

It is not difficult to see that the operator-specific distribution of SBM is the binomial distribution Bin(𝑛, 𝑝) with 𝑛 trials and success probability 𝑝. Another common operator is the 1-bit-flip operator 1pt, which is used, for example, within Randomized Local Search (RLS). 1pt creates the offspring by flipping exactly one uniformly chosen bit. Its operator-specific distribution over the mutation strengths [0..𝑛] is hence the 1-point distribution that assigns all probability mass to ℓ = 1. Likewise, the 𝑘pt operator flips 𝑘 pairwise different, uniformly selected bits, and its operator-specific distribution is the 1-point distribution giving all probability mass to ℓ = 𝑘. Other common unary unbiased mutation operators include the "shift" SBM, SBM 0→1 , which is similar to SBM but which assign all probability weight from ℓ = 0 to ℓ = 1 and the "resampling" SBM, SBM >0 , which assigns the probability weight of sampling mutation strength 0 proportionally to all positive mutation strengths 1 ≤ ℓ ≤ 𝑛 by assigning to each of these values probability Bin(𝑛, 𝑝)(ℓ)/(1-(1-𝑝) 𝑛 ); see [START_REF] Carvalho | A Simple Proof for the Usefulness of Crossover in Black-Box Optimization[END_REF][START_REF] Carvalho | Towards a More Practice-Aware Runtime Analysis of Evolutionary Algorithms[END_REF] for a discussion and motivation of these two latter variants. Motivated by the observation that infrequent large "jumps" can be beneficial in evolutionary algorithm behavior, B. Doerr et al. introduced in [START_REF] Doerr | Fast genetic algorithms[END_REF] the fast genetic algorithm (GA), which samples the mutation strength from the heavy-tailed, power-law distribution

P [ℓ = 𝑘] = (𝐶 𝛽 𝑛/2 ) -1 𝑘 -𝛽 with 𝐶 𝛽 𝑛/2 = 𝑛/2
𝑖=1 𝑖 -𝛽 and 𝛽 being some constant, often set as 𝛽 = 1.5. Finally, in [START_REF] Ye | Interpolating Local and Global Search by Controlling the Variance of Standard Bit Mutation[END_REF] a normalized mutation operator was suggested, which samples the mutation strength from a normal distribution N (𝜇, 𝜎 2 ). evaluate 𝑓 (𝑦 (𝑖) );

In contrast to the examples discussed above, this operators allows

7 𝑦 ← select arg max{𝑓 (𝑦 (𝑖) ) | 𝑖 ∈ [𝜆]} ; 8 if 𝑓 (𝑦) ≥ 𝑓 (𝑥) then 𝑥 ← 𝑦;
to scale the mean and the variance of the distribution independently of each other.

The characterization from [13, Lemma 1], which identifies mutation operators via their distributions over the set [0..𝑛], is classically only used to verify that a certain operator is unbiased. We use it here the other way around, by asking ourselves how different optimal mutation operators are from those commonly studied in the evolutionary computation community.

The (1 + 𝜆) UUSD-EA. We study this question in the context of the (1 + 𝜆) unary unbiased static distribution EA (UUSD-EA), which is given by Algorithm 1. The (1 + 𝜆) UUSD-EA is initialized uniformly at random. In each iteration, it samples 𝜆 points, which are all sampled from the same unary unbiased mutation operator. We denote the distribution from which the mutation strength is sampled by 𝐷 (𝑛, 𝜆) to indicate that it may depend on 𝑛 and 𝜆, but not on any information accumulated during the run of the algorithm. That is, the (1 + 𝜆) UUSD-EA allows only static mutation operators. For creating the 𝜆 search points that shall be evaluated in the current iteration, the (1 + 𝜆) UUSD-EA samples for each one of them a mutation strength ℓ (𝑖) and then creates the 𝑖-th "offspring" by applying the flip ℓ (𝑖) operator, which flips ℓ (𝑖) pairwise different, uniformly chosen bits in the input 𝑥. The best of these offspring replaces 𝑥 if it is at least as good as it. It is irrelevant for the context of our work how the ties are broken in line 8, as our results apply to all tie-breaking rules.

OneMax: We focus on OneMax, i.e., our goal is to determine the optimal static unary unbiased distributions for maximizing the function 𝑛 𝑖=1 𝑥 𝑖 . In the context of our work, this problem is equivalent to that of minimizing the Hamming distance to an arbitrary bit string 𝑧 ∈ {0, 1} 𝑛 [START_REF] Per | Black-Box Search by Unbiased Variation[END_REF].

Expected Runtimes: As common in literature on theory of evolutionary computation, we understand as an optimal distribution the one that minimizes the expected runtime, which we measure here in terms of generations. Since the offspring population size 𝜆 is fixed, this does not impact our results: a parallel runtime of 𝑇 generations corresponds to a runtime of exactly (𝑇 -1)𝜆 + 1 function evaluations. As runtime we therefore denote the number of iterations that the algorithm performs until it evaluates an optimal solution for the first time.

(Non-)Uniqueness of the Optimal Distributions: We note that we do not have any guarantee at the moment that the optimal distributions are unique. In fact, we observe that for any fixed problem dimension 𝑛, there is a certain threshold 𝜆(𝑛) before which small differences in the distributions cause measurable effects on the expected running times, so that the optimal distribution seems to be unique from the point of view of our computations. After the threshold, however, differences between the distributions have no measurable effect on the expected runtime, so that our algorithm may consider different ones as optimal, and may hence not always return the same distribution.

Notation: For combinations of 𝑛 and 𝜆 for which the optimal distributions are unique, we denote by 𝑃 * (𝑘 | 𝑛, 𝜆) the probability that this distribution assigns to flipping 𝑘 bits, 𝑘 ∈ [𝑛]. To ease the reading, we sometimes use the same notation also for those distributions which yield expected running times that deviates only negligibly from what appears to be the true optimum.

ALGORITHM FOR COMPUTING THE OPTIMAL DISTRIBUTIONS

Our algorithm to compute the optimal static unary unbiased distributions is based on the dynamic programming approach from [START_REF] Buzdalov | Optimal Mutation Rates for the (1 + 𝜆) EA on OneMax[END_REF]. We start the description by explaining the principles of dynamic programming that are used to compute the expected running time 𝑇 𝑓 of the (1+𝜆) UUSD-EA, measured in iterations, assuming it starts at fitness 𝑓 and the values 𝑇 𝑓 ′ are known for all 𝑓 ′ > 𝑓 . We note then that a (practical) analytic solution of the problem of finding an optimal distribution is quite unlikely to exist even for small values of 𝑛, and instead give a black-box minimization scheme with the use of a separable CMA-ES [START_REF] Ros | A simple modification in CMA-ES achieving linear time and space complexity[END_REF], a simplified and more computationally efficient version of the well-known continuous optimizer [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]. We complete with an investigation of convergence properties, which allows us to say, with great confidence, that separable CMA-ES finds a globally optimal distribution in a constant fraction of runs.

Dynamic Programming on Expected Times

We first explain how to compute 𝑇 𝑓 for a given distribution 𝐷 =

(𝐷 𝑘 ) 𝑘 ∈ [0..𝑛] .
We begin with computing the probabilities 𝑆 𝑘,𝑔 of sampling an offspring with fitness 𝑔 by flipping exactly 𝑘 bits chosen uniformly without replacement in a solution of fitness 𝑓 . Note that these quantities depend only on the current fitness 𝑓 and the problem properties, that is, they depend neither on 𝜆 nor on

𝐷. It holds from simple combinatorics that 𝑆 𝑘,𝑔 = 𝑛-𝑓 𝑔-𝑓 +𝑖 𝑓 𝑖 / 𝑛
𝑘 , where we assume this probability to be zero if one of the binomial coefficient arguments are out of bounds.

Next we compute the probabilities (𝑄

𝑔 ) 𝑔=𝑓 +1,𝑓 +2,...,𝑛 of sampling a single offspring of fitness 𝑔. For 𝑔 > 𝑓 they are derived from 𝑆 𝑘,𝑔 by using the distribution parameters (𝐷 𝑘 ) 𝑘 ∈ [0..𝑛] as follows:

𝑄 (1) 𝑔 = ∑ 𝑛 𝑘=1 𝐷 𝑘 𝑆 𝑘,𝑔 .
As the (1 + 𝜆) UUSD-EA is an elitist algorithm, and the behavior of the algorithm with different parents having the same fitness value is the same, with the remaining probability the algorithm remains in the same state, which we capture as 𝑄 The probability of each possible fitness improvement after sampling all the 𝜆 offspring is then computed as follows:

𝑄 (𝜆) 𝑔 = ∑ 𝑔 𝑖=𝑓 𝑄 (1) 𝑖 𝜆 - ∑ 𝑔-1 𝑖=𝑓 𝑄 (1) 𝑖 𝜆 . (1)
In simple words, the new fitness is 𝑔 if all offspring have the fitness in [𝑓 ..𝑔], but not all of them have the fitness in [𝑓 ..𝑔 -1], counting all offspring with fitness smaller than 𝑓 towards 𝑓 . Finally, the expected time to reach the optimum from the fitness 𝑓 is computed using the following expression:

𝑇 𝑓 = 1 1 -𝑄 (𝜆) 𝑓 1 + ∑ 𝑛 𝑔=𝑓 +1 𝑄 (𝜆) 𝑔 𝑇 𝑔 , (2)
by standard arguments as detailed in [START_REF] Buskulic | Maximizing drift is not optimal for solving OneMax[END_REF][START_REF] Buzdalov | Optimal Mutation Rates for the (1 + 𝜆) EA on OneMax[END_REF].

The time and memory complexities of such a step is 𝑂 (𝑛 2 ) for computing 𝑆 𝑘,𝑔 and 𝑂 (𝑛) for the other stages. As there are 𝑛 different fitness values for 𝑓 , the time complexity of computing all the expected running times for one set of algorithm parameters (𝜆, 𝐷) is 𝑂 (𝑛 3 ), whereas the memory complexity is still 𝑂 (𝑛 2 ) as the 𝑆 𝑘,𝑔 matrix may be discarded once 𝑓 changes. However, since 𝑆 𝑘,𝑔 depend only on 𝑛 and 𝑓 , one may evaluate up to 𝑂 (𝑛) different combinations of algorithm parameters (𝜆, 𝐷) in a single run by merging the activities corresponding to each 𝑓 , which preserves the total time complexity of 𝑂 (𝑛 3 ) and the memory complexity of 𝑂 (𝑛 2 ). This feature will turn later to be beneficial when a population-based optimizer is applied to find the best possible distribution 𝐷.

Optimization with Separable CMA-ES

When solving the problem of finding an optimal distribution 𝐷, one could choose to express each of 𝑇 𝑓 as a function of 𝑛 + 1 distribution parameters (𝐷 𝑘 ) 𝑘 ∈ [0..𝑛] , promote such expressions through dynamic programming, take their weighed sum for the expected running time from a random point and perform analytical optimization using standard analysis approaches. However, the presence of 𝑄 (𝜆) 𝑓 in the denominator in (2) makes the resulting expression nonlinear even for 𝜆 = 1, and having a 𝜆 in the exponent in [START_REF] Anonymized | Code and data repository for this paper[END_REF] makes it even harder. The resulting expression appears to be a ratio of polynomials of degree Θ(𝑛𝜆) with 𝐷 𝑘 as variables, which makes it infeasible to perform an exact analytical minimization even for small problem sizes.

Such large degree of the polynomials also effectively prevents gradient-based optimization, since the exact numeric computation of the derivatives -although possible -would require considerable computational resources, which furthermore significantly increase with 𝜆. At the same time, evaluating the expected runtime of the (1 + 𝜆) UUSD-EA for a given distribution does not depend on 𝜆, assuming the computations are done in machine precision. This makes it possible to apply black-box optimization techniques to identify distributions which minimize the expected runtime. From the large set of possible black-box optimization techniques, we chose the CMA-ES family of algorithms [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] that are well suitable for such kind of optimization. The application of CMA-ES requires a few clarifications, which we list below. programming platform was due to the computational complexity of the fitness evaluation. On one hand, it is too expensive to be implemented in Python or Matlab (the languages with reference implementations of CMA-ES maintained by the authors of this algorithm). On the other hand, the costs of inter-process communication are not negligible compared to fitness evaluation, so evaluation of fitness shall happen in the same process as the optimizer.

• To search the distributions, we tune the CMA-ES to respect box constraints (each variable is in [0, 1]) and, before fitness evaluation, we normalize variable values in the Lamarckian sense, that is, without updating the individuals themselves, which would clash with the assumptions made by CMA-ES.

• As we perform distribution optimization for a noiseless problem, we can safely set 𝐷 0 to zero, which restricts the search to the 𝑛-dimensional cube [0, 1] 𝑛 .

We had to modify the implementation of CMA-ES, as the particular class hierarchy in Apache Commons Math does not allow to evaluate the whole population in a single call, which we benefit from. Besides, after some experimentation, we switched to the separable CMA-ES [START_REF] Ros | A simple modification in CMA-ES achieving linear time and space complexity[END_REF], which is more efficient in terms of computational costs, but also produced better results (likely due to less operations that caused precision loss). We also vastly optimized its implementation, which resulted in an 8× reduction of the wall-clock running time.

The configurable parameters of CMA-ES were: population size 10, initial step size 1.0, random initial guess, computational budget 100𝑛 2 . The optimizer, however, did not reach the computational budget, as, in all runs, it converged to a single point and terminates at one of the degeneration criteria much earlier than that.

Convergence Analysis

The preliminary experiments showed that CMA-ES typically converges relatively quickly to nearly the same value in most of runs. We observed that in a constant fraction of runs the algorithm optimized the distribution to produce the same expected running time up to precision of 10 -12 and better. Fig. 2 shows examples for some values of 𝜆 for 𝑛 = 3 and 𝑛 = 100. Note that there is no curve for 𝜆 = 1 in the 𝑛 = 100 plot since in this case all the values were identical. Since the delivered precision is very close to the precision available for the 64-bit floating point machine numbers, we assume that CMA-ES reaches the global optimum of the problem in most of the cases. For further analyses we selected only the runs that are at most 1 + 10 -9 times greater than the obtained minimum.

We have also performed the robustness analysis for the distributions produced by the optimizer. For each 𝑛, 𝜆, and 𝑘, the standard deviation of the values 𝑃 * (𝑘 | 𝑛, 𝜆) was computed across all the good enough distributions. The results are presented in Fig. 3, where an intriguing picture appears. First, small 𝜆 produce very small (much less than 10 -10 ) maximal standard deviations, which then jump to the region of [10 -8 ; 10 -5 ] and remain there until 𝜆 reaches a certain threshold. Above that threshold, the maximal standard deviations experience some sort of phase shift and raise to very high values reaching 0.1 and above. Note that, by our selection of the data that is considered in this computation, the (1 + 𝜆) UUSD-EA still shows nearly identical expected running times on such different distributions and values of 𝜆. We show later in the next section that this is, in fact, an expected behavior that corresponds to situations when there is a single global optimum, but a number of different distributions coincide in its running time expectation with that global optimum up to the machine precision. Fig. 4 displays the average wall-clock times required to find the optimal distribution for all available 𝑛 and few values of 𝜆. The available data suggests that the time complexity scales polynomially with the degree of 4 + 𝜀 for some small 𝜀, which, together with the earlier cubic runtime bound for the evaluation of a given distribution, suggests that the CMA-ES requires 𝑂 (𝑛 1+𝜀 ) iterations before hitting one of its termination criteria.

OPTIMAL DISTRIBUTIONS

We now take a closer look at the distributions returned by the algorithm from Sec. 3. As mentioned, the data is available at [START_REF] Anonymized | Code and data repository for this paper[END_REF], and we present here only our main findings.

We first study the impact of 𝜆 on the optimal distributions. To this end, we fix the dimension and analyze how for each possible mutation strength 𝑘 ∈ [1..𝑛] its probability of being chosen depends on 𝜆. that is often used as baseline for comparisons, both in empirical [START_REF] Doerr | Benchmarking discrete optimization heuristics with IOHprofiler[END_REF] and in theoretical [START_REF] Doerr | Theory of Evolutionary Computation: Recent Developments in Discrete Optimization[END_REF] research. Our data shows that for 𝑛 ∈ {8, 11} the generalized (1 + 𝜆) RLS is optimal for 𝜆 ≤ 3 and suboptimal for 𝜆 = 4. Similarly, for 𝑛 ∈ {16, 23, 32} it is optimal for 𝜆 ≤ 4, and for 𝑛 ∈ {64, 91, 100} the threshold is 𝜆 = 5. We are confident that this describes a general trend of a positive correlation between the dimension 𝑛 and the maximal 𝜆 for which the 1-point distribution Pr

[ℓ = 1] = 1 is optimal.
The conditional binomial distribution Bin >0 (𝒏, 1/2) is arbitrarily close to an optimal one for large 𝝀. We plot in Fig. 6 the optimal probability 𝑃 * (𝑘 | 𝑛 = 11, 𝜆) for large 𝜆. The curves for 𝜆 = 512 and 𝜆 = 1 024 approximate the conditional binomial distribution Bin >0 (11, 1/2), which we plot as dotted red line. This can be explained as follows: the (1 + 𝜆) UUSD-EA with the static mutation strength distribution Bin >0 (𝑛, 1/2) is simply the random search algorithm, which samples all search points, except the parent, uniformly at random. For 𝜆 = Ω(2 𝑛 ), this algorithm has a very good chance of sampling every point 𝑦 ∈ {0, 1} 𝑛 , so that it also has a decent chance of hitting the optimum. When 𝜆 is not much larger than 2 𝑛 , the quality of a distribution is significantly influenced by each of its parameters, so for these cases our computations provide distributions that are almost identical to random search in each independent run. When 𝜆 is much bigger than 2 𝑛 , however, the situation changes: while the common sense suggests that the truly optimal distribution gets even closer to pure random search, the quality of a distribution quickly loses sensitivity to its parameters, and we obtain different distributions, which all yield practically indistinguishable expected runtimes. We provide an example for 𝑛 = 3 in Tab. 1, where we observe that for 𝜆 ≤ 128 the average of the computed static mutation strength distributions converge against Bin >0 (𝑛, 1/2). The standard deviation of the independent runs of our optimizer are negligible in this regime. For 𝜆 = 256, 512, and 1 024, however, the maximum standard deviation among the three mutation strengths are 0.07, 0.12, and 0.17, respectively. Fig. 7 plots the results of all 50 independent runs for 𝜆 = 1 024.

𝑷 * (1 | 𝝀, 𝒏) decreases monotonically with increasing 𝝀. For fixed dimension 𝑛, the importance of 1-bit flips significantly Table 1: Average values of the recommended distributions 𝑃 * (𝑘 | 𝑛 = 3, 𝜆) for different 𝜆 and all mutation strengths 𝑘.

For 𝜆 ≤ 128 the standard deviations are negligible, then grow quickly (see Figure 3). Unreliable values are indicated with gray font. We add Bin >0 (3, 1/2) for comparison. decreases as 𝜆 grows. We have seen this for 𝑛 = 16 in Fig. 5. This trend generalizes to other problem dimensions, which can be seen in Fig. 8, where we plot the probabilities 𝑃 * (1 | 𝜆, 𝑛) and 𝑃 * (2 | 𝜆, 𝑛) for 𝑛 ∈ {8, 11, 16, 23, 32, 45, 64}.

𝜆
𝑷 * (2 | 𝝀, 𝒏) is non-monotonic in 𝝀. We clearly see from Fig. 8 that, for fixed dimension 𝑛, the optimal probability to flip two bits is non-monotonic in 𝜆. The 𝜆 at which it becomes non-zero appears to be monotonic in 𝑛, however, the set of 𝜆 we used is not enough to determine the exact threshold. It is 𝜆 = 8 for 𝑛 = 8, 8 < 𝜆 ≤ 16 for 𝑛 ∈ {11, 16, 23}, and it is 16 < 𝜆 ≤ 32 for 𝑛 ∈ {32, 45, 64, 91, 100}. Surprisingly enough, this threshold is always larger than the value at which 𝑃 * (1 | 𝜆, 𝑛) becomes less than one. We do not see any pattern in the points at which Pr[𝑘 = 2] starts to decrease again.

Flipping all bits can be optimal. Intuitively, flipping more than 𝑛/2 bits can be optimal in an elitist algorithm only when OM(𝑥) < 𝑛/2. It may therefore be surprising that even for mutation strength ℓ = 𝑛 (i.e., flipping all bits) the optimal static probability can be non-zero already for comparatively small 𝜆. In Tab. 2 we show for which combinations of 𝑛 and 𝜆 the optimal probability of flipping all bits is non-zero. Note that for some 𝑛 we have not seen any 𝜆 where this probability is nonzero, and this may be related to the parity of 𝑛: for 𝑛 ∈ {23, 45, 91} the small 𝜆 feature rather a nonzero 𝑃 * (𝑛 -1 | 𝑛, 𝜆) instead. So far we do not have any explanation for the observed patterns. We take a set Λ = {6, 7, 8, 16, 32, 64} and show an example in Fig. 9, where we display 𝑃 * (𝑘 | 𝑛 = 64, 𝜆 ∈ Λ) for all 𝑘 ∈ [50, 64] for which there exists at least one 𝜆 ∈ Λ with 𝑃 * (𝑘 | 𝑛 = 64, 𝜆) > 0. Interestingly, 𝑃 * (𝑘 | 𝑛 = 64, 𝜆) = 0 for all tested 𝜆 > 64. We also see that only nine different 𝑘 appear out of [50, 64], of which at most three have a non-zero optimal probability for any of the tested 𝜆.

The number of mutation strengths 𝒌 with 𝑷 * (𝒌 | 𝝀, 𝒏) > 0 increases with 𝝀, but not monotonically. We summarize in Tab. 3 the number of different mutation strengths 𝑘 for which 𝑃 * (𝑘 | 𝜆, 𝑛) > 0. While there is a general trend for an increasing number of such 𝑘 with increasing 𝜆, these trends are non-monotonic. From the previous insights, however, it is clear that for every 𝑛 there exists a threshold 𝜆(𝑛) such that for all 𝜆 > 𝜆(𝑛) the number of 𝑘 with 𝑃 * (𝑘 | 𝜆, 𝑛) > 0 is equal to 𝑛. A closer look for small 𝝀. It is well known that small values of 𝜆 are preferable for the optimization of simple function such as OneMax [START_REF] Jansen | On the Choice of the Offspring Population Size in Evolutionary Algorithms[END_REF]. We therefore take a more detailed look at small values in Tab. [START_REF] Badkobeh | Unbiased Black-Box Complexity of Parallel Search[END_REF], where we list for 𝑛 = 11, 𝜆 ∈ [4..8]∪{16, 32, 64}, and for all possible mutation strengths 𝑘 ∈ [START_REF] Doerr | Theory of Parameter Control Mechanisms for Discrete Black-Box Optimization: Provable Performance Gains Through Dynamic Parameter Choices[END_REF] the optimal probability 𝑃 * (𝑘 | 𝑛, 𝜆). We recall from above that for 𝜆 ≤ 3 the one-point distribution assigning all weight to 𝑘 = 1 is optimal. For 𝜆 = 4 and 𝜆 = 7, the optimal probability 𝑃 * (1 | 𝑛 = 11, 𝜆) is slightly less than 1, and the remaining probability mass is assigned to 𝑘 = 10, i.e., to the operator flipping all bits but one. For 𝜆 = 8 two more values, 𝑘 = 4 and 𝑘 = 11, are also chosen with positive probability. As already discussed in the context of Tab. 3 and Fig. 6, the number of mutation strengths 𝑘 with 𝑃 * (𝑘 | 𝑛 = 11, 𝜆) > 0 increases for increasing 𝜆, and the distribution converges to the conditional binomial distribution Bin >0 (11, 1/2), which we include in the table for reference.

RUNTIME COMPARISON

After having focused on the distributions in the previous section, we now study the runtime of the optimal (1 + 𝜆) UUSD-EA in comparison to other common (1 + 𝜆) UUSD-EAs. We include in our comparison the (1 + 𝜆) EA variants with SBM, SBM >0 , and SBM 0→1 standard bit mutation operators, the (1 + 𝜆) RLS, and the (1 + 𝜆) fastGA with different 𝛽 ∈ {1.3, 1.5, 1.7}. We also considered the variant of the latter which directly samples mutation strengths ℓ proportional to ℓ -𝛽 for the same parameter values.

It is not difficult to see that, for any fixed 𝑛, the expected runtime of the optimal (1 + 𝜆) UUSD-EAs converge to 1 -1/2 𝑛 as 𝜆 → ∞. This is also the case for all (1 + 𝜆) UUSD-EA variants that assign a positive probability to each positive mutation strength 𝑘 ∈ [𝑛], and this for all problems 𝑓 : {0, 1} 𝑛 → R. Since SBM and fast mutation satisfy these requirements, the expected runtime of all (1 + 𝜆) EA variants as well as that of the fastGA variants also converges to 1 -1/2 𝑛 , but at a possibly much different speed. The expected runtime of the generalized (1 + 𝜆) RLS, in contrast, converges to 𝑛/2 on OneMax, since, hand-waivingly, this is the expected distance to the optimum after initialization and the probability to make a progress of 1 in each iteration converged to 1 as 𝜆 → R.

In Table 5 we present the expected runtimes of the optimal

(1 + 𝜆) UUSD-EA(s) on OneMax for all different combinations of 𝑛 and 𝜆 we have. Note that these numbers are the lower bounds for all

(1 + 𝜆) UUSD-EAs, including the algorithms mentioned above. Note also that algorithms obtaining a better expected runtime require adaptive parameter choices.

In Figures 10 we illustrate the regret in the expected runtime of common (1 + 𝜆) UUSD-EAs compared to the optimal one, for 𝑛 ∈ {3, 100}, respectively. Corresponding to our discussion above, we observe that for 𝑛 = 3 all algorithms, with the exception of RLS, converge to the same optimal expected runtime of 0.875 = 1 -1/8 when 𝜆 → ∞, whereas the generalized (1 + 𝜆) RLS converges to 1.5, which corresponds to a factor 12/7 ≈ 1.71... relative disadvantage.

Note that the regrets are not monotone for those algorithms which never flip zero bits, except for the (1 + 𝜆 RLS. As 𝜆 grows from the small values, their regret decreases, most probably as flipping more than one bit becomes a better choice. However, with further increase of 𝜆, the fact that these algorithms flip many bits with a small probability turns to a disadvantage. Note how the heavy-tailed algorithm that samples the mutation strength directly from a power-law distribution with the smallest tested 𝛽 = 1.3 becomes a clear winner at 𝜆 ≥ 7. The similar behavior is seen also for 𝑛 = 100 with the exception that the values of 𝜆 are not large enough to observer the complete convergence picture.

CONCLUSION

We have analyzed in this paper the dependence of 𝑃 * (𝑘 | 𝑛, 𝜆), the optimal probability of flipping 𝑘 bits in the (1 + 𝜆) EA-UUSD, in dependence of 𝑛 and 𝜆. Among other insights, we have shown that the (1 + 𝜆) RLS is optimal when 𝜆 is small, and that the value for which it ceases to be optimal increases with increasing 𝑛. We have also seen that, for fixed 𝑛, the distribution 𝑃 * (𝑘 | 𝑛, 𝜆) converges to the conditional binomial distribution Bin >0 (𝑛, 1/2) when 𝜆 → ∞.

For future work, we consider the following particularly exciting.

1) Formalizing the observations into rigorous results: We believe that some of the observations made in this work could be formalized with reasonable effort.

2) Analyzing benefits of generalized mutation for more complex problems: For practitioners, our work is perhaps most interesting in that it invites to consider mutation operators through the lens of probability distributions over the set of possible radii. This idea should show its full potential on problems that are more complex than OneMax. The fastGA proposed in [START_REF] Doerr | Fast genetic algorithms[END_REF] and fast AIS from [START_REF] Corus | Fast Artificial Immune Systems[END_REF] are compelling examples that shows that such generalization can indeed be very beneficial [START_REF] Meunier | Black-Box Optimization Revisited: Improving Algorithm Selection Wizards through Massive Benchmarking[END_REF][START_REF] Mironovich | Evaluation of heavy-tailed mutation operator on maximum flow test generation problem[END_REF].

3) Transferring the generalizations to variation operators of higher arity: The quest for analyzing more general variation operators is not restricted to mutation alone, but also generalizes to variation operators of higher arity, called "crossover" or "recombination" operators in evolutionary computation. In a simplest extension, one could study effects of changing the binomial distribution associated with the number of bits that is taken from each parent in uniform crossover. We did not yet investigate this idea further, but we hope that a de-coupling of mean and variance, similarly as proposed for variation in [START_REF] Ye | Interpolating Local and Global Search by Controlling the Variance of Standard Bit Mutation[END_REF], may be beneficial.

4) Interplay of generalized mutation with other variation operators:

The benefits of generalized mutation operators are very likely not restricted to mutation-only algorithms, but could also improve algorithms that use variation operators of different arities. First examples demonstrating clear advantages of heavy-tailed mutation in the (1 + (𝜆, 𝜆)) GA [START_REF] Doerr | From black-box complexity to designing new genetic algorithms[END_REF] were recently shown in [START_REF] Antipov | Fast mutation in crossover-based algorithms[END_REF][START_REF] Antipov | Runtime Analysis of a Heavy-Tailed (1 + (𝜆, 𝜆)) Genetic Algorithm on Jump Functions[END_REF].

5) Extensions to the dynamic case:

We studied in this work the case of static distributions 𝑃 * (𝑘 | 𝑛, 𝜆), whereas it is well known that a dynamic choice of the mutation rates, or algorithms' parameters in general, can lead to significant performance gains [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]. Combining the analyses made in [START_REF] Buzdalov | Optimal Mutation Rates for the (1 + 𝜆) EA on OneMax[END_REF] for the optimal dynamic 1pt and the optimal dynamic SBM operators with the approach taken in this work (the generalization to arbitrary distributions) would provide an exact quantification of the disadvantage of static against dynamic mutation operator choices.

Table 5: Expected runtimes (in generations) of the optimal (1 + 𝜆) UUSD-EA(s) on OneMax for different combinations of 𝑛 (rows) and 𝜆 (columns), rounded to two digits. For every fixed 𝑛, the optimal expected runtime converges to 1 -1/2 𝑛 as 𝜆 → ∞.

𝜆

Algorithm 1 : 5 𝑦

 15 The (1+𝜆) unary unbiased static distribution EA (UUSD-EA) maximizing a function 𝑓 : {0, 1} 𝑛 → R. 1 Initialization: Sample 𝑥 ∈ {0, 1} 𝑛 uniformly at random and evaluate 𝑓 (𝑥); 2 Optimization: for 𝑡 = 1, 2, 3, . . . do 3 for 𝑖 = 1, . . . , 𝜆 do 4 Sample ℓ (𝑖) ∼ 𝐷 (𝑛, 𝜆); (𝑖) ← flip ℓ (𝑖) (𝑥);6

( 1 )

 1 𝑓 = 1 -𝑔=𝑓 +1,𝑓 +2,...,𝑛 𝑄(1)𝑔 .

Figure 1 :

 1 Figure 1: Fitness as a function of the number of evaluations for 𝑛 = 16 and 𝜆 = 8. First four runs are shown.

Fig. 1

 1 shows example runs for 𝑛 = 16 and 𝜆 = 8. Runs for all 𝑛 and 𝜆 demonstrate similar behavior. The rest of the paper is based on the data collected for 𝑛 ∈ {3, 5, 8, 11, 16, 23, 32, 45, 64, 91, 100} and 𝜆 ∈ [1..7] ∪ {2 𝑖 | 3 ≤ 𝑖 ≤ 10}. For each (𝑛, 𝜆) pair we performed 50 independent runs of the CMA-ES.

Figure 2 :Figure 3 :

 23 Figure 2: Relative loss of the expected runtime induced by the result of the 50 CMA-ES runs, when compared against the best expected runtime. Plots are for 𝑛 = 3 (left) and 𝑛 = 100 (right), respectively, values are sorted by increasing loss.

4 4

 4 • 10 -6 • 4.2

Figure 4 :

 4 Figure 4: Average wall-clock times for various 𝜆 together with guesses for asymptotic bounds.

9 Figure 5 :

 95 Figure 5: Probabilities of flipping 𝑘 bits in optimal static mutation strength distributions for different population sizes 𝜆 for optimizing OneMax in dimension 𝑛 = 16.

Fig. 5 Figure 6 :

 56 Figure 6: Values 𝑃 * (𝑘 | 𝑛 = 11, 𝜆) for large 𝜆 plotted against 𝑘. The larger 𝜆, the more the curves resemble the binomial distribution Bin >0 (11, 1/2), which is also plotted.

Figure 7 := 1 = 8 , = 2 = 11 , = 1 = 11 , = 2 = 16 , = 1 = 16 , = 2 = 23 , = 1 = 23 , = 2 = 32, = 1 = 32, = 2 = 45, = 1 = 45, = 2 = 64, = 1 = 64, = 2 Figure 8 :

 71821111121611622312321212128 Figure 7: 𝑃 * (𝑘 | 𝑛 = 𝜆 = 024) for each of the 50 independent runs of our optimizer. The expected optimization time is identical for all of them (0.875).

Figure 9 :

 9 Figure 9: 𝑃 * (𝑘 | 𝑛 = 64, 𝜆), in dependence of 𝜆, and for all values 𝑘 ∈ [50, 64] for which there is at least one 𝜆 such that 𝑃 * (𝑘 | 𝑛, 𝜆) > 0. Note that, for all these 𝑘, 𝑃 * (𝑘 | 𝑛, 𝜆) = 0 for 𝜆 ∈ {128, 256, 512, 1024}.

5 ( 1 + 7 ( 1 Figure 10 :

 517110 Figure 10: Regrets of expected runtime of different (1 + 𝜆) UUSD-EAs on OneMax, in dependency of 𝜆, compared to the runtime of the optimal (1 + 𝜆) UUSD-EA computed in Sec. 4, for 𝑛 = 3 (top) and 𝑛 = 100 (bottom). The algorithms that sample the mutation strength ℓ directly from the powerlaw distribution with parameter 𝛽 are denoted as ℓ ∼ pow(𝛽). The notation for other algorithms is standard.

•

  Our implementation of CMA-ES is based on the one from Apache Commons Math, version 3.6.1. The choice of the Java

	10 -2
	10 -6	0	2,000	4,000	6,000	8,000

Table 2 :

 2 Combinations of 𝑛 and 𝜆 for which 𝑃 * (𝑛 | 𝑛, 𝜆) > 0 , i.e., the optimal probability of flipping all bits is non-zero.

	𝑛	3 4	5	6	7	8	𝜆 16 32 64 128 256 512 1024
	8 11 16 23 32 45 64 91 100	+ + + + + + + + + + + + + + + + + + +	+ +	+	+ +	+ + +	+ +	+ + +
	Probability	0.01 0.02		= 52 = 55 = 59 = 61 = 63		= 54 = 57 = 60 = 62 = 64
		0.00		2 3			2 4		2 5	2 6

Table 3 :

 3 Number of different mutation strengths 𝑘 for which 𝑃 * (𝑘 | 𝜆, 𝑛) > 0, in dependence of 𝑛 and 𝜆.

Table 4 :

 4 Values of 𝑃 * (𝑘 | 𝑛 = 11, 𝜆) for selected 𝜆.

	𝜆 7 1 0.9919 0.9459 0.9156 0.8937 0.8555 0.6337 0.4797 0.3483 𝑘 4 5 6 8 16 32 64 Bin >0 (11, 1 2 ) 0.0054 2 0.2483 0.3241 0.3714 0.0269 3 0.0806 4 0.0530 0.1612 5 0.2257 6 0.2257 7 0.1068 0.1865 0.1612 8 0.0162 0.0677 0.0938 0.0806 9 0.0748 0.0217 0.0269 10 0.0081 0.0541 0.0844 0.1063 0.0875 0.0270 0.0054 11 0.0040 0.0005