Evolution strategies are evolutionary algorithms that date back to the 1960s and that are most commonly applied to black-box optimization problems in continuous search spaces. Inspired by biological evolution, their original formulation is based on the application of mutation, recombination and selection in populations of candidate solutions. From the algorithmic viewpoint, evolution strategies are optimization methods that sample new candidate solutions stochastically, most commonly from a multivariate normal probability distribution. Their two most prominent design principles are unbiasedness and adaptive control of parameters of the sample distribution. In this overview the important concepts of success based step-size control, selfadaptation and derandomization are covered, as well as more recent developments like covariance matrix adaptation and natural evolution strategies. The latter give new insights into the fundamental mathematical rationale behind evolution strategies. A broad discussion of theoretical results includes progress rate results on various function classes and convergence proofs for evolution strategies.

Overview

Evolution Strategies [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF][START_REF] Rechenberg | Evolutionsstrategie '94[END_REF][START_REF] Schwefel | Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie[END_REF][START_REF] Schwefel | Evolution and Optimum Seeking[END_REF], sometimes also referred to as Evolutionary Strategies, and Evolutionary Programming [START_REF] Fogel | Artificial Intelligence through Simulated Evolution[END_REF] are search paradigms inspired by the principles of biological evolution. They belong to the family of evolutionary algorithms that address optimization problems by implementing a repeated process of (small) stochastic variations followed by selection: in each generation (or iteration), new offspring (or candidate solutions) are generated from their parents (candidate solutions already visited), their fitness is evaluated, and the better offspring are selected to become the parents for the next generation.

Evolution strategies most commonly address the problem of continuous black-box optimization. The search space is the continuous domain, R n , and solutions in search space are n-dimensional vectors, denoted as x. We consider an objective or fitness function f : R n → R, x → f (x) to be minimized. We make no specific assumptions on f , other than that f can be evaluated for each x, and refer to this search problem as black-box optimization. The objective is, loosely speaking, to generate solutions (x-vectors) with small f -values while using a small number of f -evaluations. 1 In this context, we present an overview of methods that sample new offspring, or candidate solutions, from normal distributions. Naturally, such an overview is biased by the authors' viewpoints, and our emphasis will be on important design principles and on contemporary evolution strategies that we consider as most relevant in practice or future research. More comprehensive historical overviews can be found elsewhere [START_REF] Beyer | Evolution strategies -A comprehensive introduction[END_REF][START_REF] Fogel | Evolutionary Computation: The Fossil Record[END_REF].

In the next section the main principles are introduced and two algorithm templates for an evolution strategy are presented. Section 3 presents six evolution strategies that mark important conceptual and algorithmic developments. Section 4 summarizes im-1 Formally, we like to "converge" to an essential global optimum of f , in the sense that the best f (x) value gets arbitrarily close to the essential infimum of f (i.e., the smallest f -value for which all larger, i.e. worse f -values have sublevel sets with positive volume). c µ/µ,λ the progress coefficient for the (µ/µ, λ)-ES [START_REF] Beyer | The Theory of Evolution Strategies[END_REF] equals the expected value of the average of the largest µ order statistics of λ independent standard normally distributed random numbers and is in the order of 2 log(λ/µ).

C ∈ R n×n a (symmetric and positive definite) covariance matrix

C 1 /2 ∈ R n×n a matrix that satisfies C 1 /2 C 1 /2 = C and is symmetric if not stated otherwise. If C 1 /2
is symmetric, the eigendecomposition C

1 /2 = BΛB with BB = I and diagonal matrix Λ exists and we find C = C

1 /2 C 1 /2 = BΛ 2 B as eigendecomposition of C. e i the ith canonical basis vector f : R n → R fitness or objective function to be minimized I ∈ R n×n the identity matrix (identity transformation) i.i.d. independent and identically distributed N (x, C) a multivariate normal distribution with expectation and modal value x and covariance matrix C, see Section 2.4.

n ∈ N search space dimension P a multiset of individuals, a population s, s σ , s c ∈ R n a search path or evolution path s, s k endogenous strategy parameters (also known as control parameters) of a single parent or the kth offspring; they typically parametrize the mutation, for example with a step-size σ or a covariance matrix C t ∈ N time or iteration index w k ∈ R recombination weights

x, x (t) , x k ∈ R n solution or object parameter vector of a single parent (at iteration t) or of the kth offspring; an element of the search space R n that serves as argument to the fitness function f : R n → R.

diag : R n → R n×n the diagonal matrix from a vector exp α : R n×n → R n×n , A → ∞ i=0 (αA) i / i! is the matrix exponential for n > 1, otherwise the exponential function. If A is symmetric and BΛB = A is the eigendecomposition of A with BB = I and Λ diagonal, we have exp(A) = B exp(Λ)B = B ∞ i=0 Λ i /i! B = I + BΛB + BΛ 2 B /2 + . . . . Furthermore we have exp α (A) = exp(A) α = exp(αA) and exp α (x) = ( e α ) x = e αx .

Main Principles

Evolution strategies derive inspiration from principles of biological evolution. We assume a population, P, of so-called individuals. Each individual consists of a solution or object parameter vector x ∈ R n (the visible traits) and further endogenous parameters, s (the hidden traits), and an associated fitness value, f (x). In some cases the population contains only one individual. Individuals are also denoted as parents or offspring, depending on the context. In a generational procedure, 1. one or several parents are picked from the population (mating selection) and new offspring are generated by duplication and recombination of these parents;

2. the new offspring undergo mutation and become new members of the population;

3. environmental selection reduces the population to its original size.

Within this procedure, evolution strategies employ the following main principles that are specified and applied in the operators and algorithms further below.

Environmental Selection is applied as so-called truncation selection. Based on the individuals' fitnesses, f (x), only the µ best individuals from the population survive. In contrast to roulette wheel selection in genetic algorithms [START_REF] Goldberg | Genetic Algorithms in Search, Optimization and Machine Learning[END_REF], only fitness ranks are used. In evolution strategies, environmental selection is deterministic. In evolutionary programming, like in many other evolutionary algorithms, environmental selection has a stochastic component. Environmental selection can also remove "overaged" individuals first.

Mating Selection and Recombination. Mating selection picks individuals from the population to become new parents. Recombination generates a single new offspring from these parents. Specifically, we differentiate two common scenarios for mating selection and recombination.

fitness-independent mating selection and recombination do not depend on the fitness values of the individuals and can be either deterministic or stochastic. Environmental selection is then essential to drive the evolution toward better solutions.

fitness-based mating selection and recombination, where the recombination operator utilizes the fitness ranking of the parents (in a deterministic way). Environmental selection can potentially be omitted in this case.

Mutation and Parameter Control. Mutation introduces small, random and unbiased changes to an individual. These changes typically affect all variables. The average size of these changes depends on endogenous parameters that change over time. These parameters are also called control parameters, or endogenous strategy parameters, and define the notion of "small", for example via the step-size σ. In contrast, exogenous strategy parameters are fixed once and for all, for example parent number µ. Parameter control is not always directly inspired by biological evolution, but is an indispensable and central feature of evolution strategies.

Unbiasedness is a generic design principle of evolution strategies. Variation resulting from mutation or recombination is designed to introduce new, unbiased "information". Selection on the other hand biases this information towards solutions with better fitness. Under neutral selection (i.e., fitness independent mating and environmental selection), all variation operators are desired to be unbiased. Maximum exploration and unbiasedness are in accord. Evolution strategies are unbiased in the following respects.

• The type of mutation distribution, the Gaussian or normal distribution, is chosen in order to have rotational symmetry and maximum entropy (maximum exploration) under the given variances. Decreasing the entropy would introduce prior information and therefore a bias.

• Object parameters and endogenous strategy parameters are unbiased under recombination and unbiased under mutation. Typically, mutation has expectation zero.

• Invariance properties avoid a bias towards a specific representation of the fitness function, e.g. representation in a specific coordinate system or using specific fitness values (invariance to strictly monotonic transformations of the fitness values can be achieved). Parameter control in evolution strategies strives for invariance properties [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF].

(µ/ρ + , λ) Notation for Selection and Recombination

An evolution strategy is an iterative (generational) procedure. In each generation new individuals (offspring) are created from existing individuals (parents). A mnemonic notation is commonly used to describe some aspects of this iteration. The (µ/ρ + , λ)-ES, where µ, ρ and λ are positive integers, also frequently denoted as (µ + , λ)-ES (where ρ remains unspecified) describes the following.

• The parent population contains µ individuals.

• For recombination, ρ (out of µ) parent individuals are used. We have therefore ρ ≤ µ.

• λ denotes the number of offspring generated in each iteration. In a (µ, λ)-ES, λ ≥ µ must hold and the case λ = µ requires fitness-based mating selection or recombination. In a (µ + λ)-ES, λ = 1 is possible and known as steady-state scenario.

Occasionally, a subscript to ρ is used in order to denote the type of recombination, e.g. ρ I or ρ W for intermediate or weighted recombination, respectively. Without a subscript we tacitly assume intermediate recombination, if not stated otherwise. The notation has also been expanded to include the maximum age, κ, of individuals as (µ, κ, λ)-ES [START_REF] Schwefel | Contemporary evolution strategies[END_REF], where 'plus'-selection corresponds to κ = ∞ and 'comma'selection corresponds to κ = 1.

Two Algorithm Templates

Template 1 gives pseudocode for the evolution strategy. Given is a population, P, of at least µ individuals (x k , s k , f (x k )), k = 1, . . . , µ. Vector x k ∈ R n is a solution vector and s k contains the control or endogenous strategy parameters, for example a success counter or a step-size that primarily serves to control the mutation of x (in Line 6). The values of s k may be identical for all k. In each generation, first λ offspring are generated (Lines 3-6), each by recombination of ρ ≤ µ individuals from P (Line 4), followed by mutation of s (Line 5) and of x (line 6). The new offspring are added to P (Line 7). Overaged individuals are removed from P (Line 8), where individuals from the same generation have, by defi- nition, the same age. Finally, the best µ individuals are retained in P (Line 9). The mutation of the x-vector in Line 6 always involves a stochastic component. Lines 4 and 5 may have stochastic components as well.

Template 1 The (µ/ρ + , λ)-ES 0 given n, ρ, µ, λ ∈ N + 1 initialize P = {(x k , s k , f (x k )) | 1 ≤ k ≤ µ} 2 while not happy 3 for k ∈ {1, . . . , λ} 4 (x k , s k ) = recombine(select mates(ρ, P)) 5 s k ← mutate s(s k ) 6 x k ← mutate x(s k , x k ) ∈ R n 7 P ← P ∪ {(x k , s k , f (x k )) | 1 ≤ k ≤
When select mates in Line 4 selects ρ = µ individuals from P, it reduces to the identity. If ρ = µ and recombination is deterministic, as is commonly the case, the result of recombine is the same parental centroid for all offspring. The computation of the parental centroid can be done once before the for loop or as the last step of the while loop, simplifying the initialization of the algorithm. Template 2 shows the pseudocode in this case.

In Template 2, only a single parental centroid (x, s) is initialized. Mutation takes this parental centroid as input (notice that s k and x k in Line 4 and 5 are now assigned rather than updated ) and "recombination" is postponed to the end of the loop, computing in Line 8 the new parental centroid. While (x k , s k ) can contain all necessary information for this computation, it is often more transparent to use x and s as additional arguments in Line 8. Selection based on f -values is now limited to mating selection in procedure recombine (that is, procedure select µ best is omitted and µ is the number of individuals in P that are actually used by recombine).

Using a single parental centroid has become the most popular approach, because such algorithms are simpler to formalize, easier to analyze and even perform better in various circumstances as they allow for maximum genetic repair (see below). All instances of evolution strategies given in Section 3 are based on Template 2.

Recombination Operators

In evolution strategies, recombination combines information from several parents to generate a single new offspring. Often, multi-recombination is used, where more than two parents are recombined (ρ > 2). In contrast, in genetic algorithms often two offspring are generated from the recombination of two parents. In evolutionary programming, recombination is generally not used. The most important recombination operators used in evolution strategies are the following.

Discrete or dominant recombination, denoted by (µ/ρ D + , λ), is also known as uniform crossover in genetic algorithms. For each variable (component of the x-vector), a single parent is drawn uniformly from all ρ parents to inherit the variable value. For ρ parents that all differ in each variable value, the result is uniformly distributed across ρ n different x-values. The result of discrete recombination depends on the given coordinate system.

Intermediate recombination, denoted by (µ/ρ I + , λ), takes the average value of all ρ parents (computes the center of mass, the centroid).

Weighted multi-recombination [START_REF] Rudolph | Convergence Properties of Evolutionary Algorithms[END_REF][START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Arnold | Weighted multirecombination evolution strategies[END_REF], denoted by (µ/ρ W + , λ), is a generalization of intermediate recombination, usually with ρ = µ. It takes a weighted average of all ρ parents. The weight values depend on the fitness ranking, in that better parents never get smaller weights than inferior ones. With equal weights, intermediate recombination is recovered. By using comma selection and ρ = µ = λ, where some of the weights may be zero, weighted recombination can take over the role of fitness-based environmental selection and negative weights become a feasible option [START_REF] Rudolph | Convergence Properties of Evolutionary Algorithms[END_REF][START_REF] Arnold | Weighted multirecombination evolution strategies[END_REF]. 2In principle, recombination operators from genetic algorithms, like one-point and two-point crossover or line recombination [START_REF] Mühlenbein | Schlierkamp-Voosen: Predictive models for the breeder genetic algorithm I. Continuous parameter optimization[END_REF] can alternatively be used. However, they have been rarely applied in evolution strategies.

In evolution strategies, the result of selection and recombination is often deterministic (namely, if ρ = µ and recombination is intermediate or weighted). This means that eventually all offspring are generated by mutation from the same single solution vector (the parental centroid) as in Template 2. This leads, for given variances, to maximum entropy because all offspring are independently drawn from the same normal distribution. 3The role of recombination in general is to keep the variation in a population high. Discrete recombination directly introduces variation by generating different solutions. Their distance resembles the distance between the parents. However, discrete recombination, as it depends on the given coordinate system, relies on separability: it can introduce variation successfully only if values of disrupted variables do not strongly depend on each other. Solutions resulting from discrete recombination lie on the vertices of an axis-parallel box.

Intermediate and weighted multi-recombination do not lead to variation within the new population as they result in the same single point for all offspring. However, they do allow the mutation operator to introduce additional variation by means of genetic repair [START_REF] Beyer | Toward a theory of evolution strategies: On the benefits of sex -The (µ/µ, λ) theory[END_REF]: recombinative averaging reduces the effective step length taken in unfavorable directions by a factor of √ µ (or √ µ w in case of weighted recombination), but leaves the step length in favorable directions essentially unchanged, see also Section 4.2. This may allow increased variation by enlarging mutations by a factor of about µ (or µ w ) as revealed in Eq. ( 16), to achieve maximal progress.

Mutation Operators

The mutation operator introduces ("small") variations by adding a point symmetric perturbation to the result of recombination, say a solution vector x ∈ R n . This perturbation is drawn from a multivariate normal distribution 4 , N (0, C), with zero mean (expected value) and covariance matrix C ∈ R n×n .

We have x+N (0, C) ∼ N (x, C), meaning that x determines the expected value of the new offspring individual. We also have x + N (0, C) ∼ x + C 1 /2 N (0, I), meaning that the linear transformation C

1 /2 generates the desired distribution from the vector N (0, I) that has i.i.d. N (0, 1) components. 5Figure 1 shows different normal distributions in dimension n = 2. Their lines of equal density are ellipsoids. Any straight section through the 2dimensional density recovers a 1-dimensional Gaussian bell. Based on multivariate normal distributions, three different mutation operators can be distinguished.

Spherical/isotropic (Figure 1, left) where the covariance matrix is proportional to the identity, i.e., the mutation distribution follows σN (0, I) with step-size σ > 0. The distribution is spherical and invariant under rotations about its mean. Below, Algorithm 1 uses this kind of mutation.

Axis-parallel (Figure 1, middle) where the covariance matrix is a diagonal matrix, i.e., the mutation distribution follows N (0, diag(σ) 2 ), where σ is a vector of coordinate-wise standard deviations and the diagonal matrix diag(σ) 2 has eigenvalues σ 2 i with eigenvectors e i . The principal axes of the ellipsoid are parallel to the coordinate axes. This case includes the previous isotropic case. Below, Algorithms 2, 3, and 4 implement this kind of mutation distribution.

General (Figure 1, right) where the covariance matrix is symmetric and positive definite (i.e.

x Cx > 0 for all x = 0), generally non-diagonal and has (n 2 + n)/2 degrees of freedom (control parameters). The general case includes the previous axis-parallel and spherical cases. Below, Algorithms 5 and 6 implement general multivariate normally distributed mutations.

In the first and the second cases, the variations of variables are independent of each other, they are uncorrelated. This limits the usefulness of the operator in practice. The third case is "incompatible" with discrete recombination: for a narrow, diagonally oriented ellipsoid (not to be confused with a diagonal covariance matrix), a point resulting from selection and discrete recombination lies within this ellipsoid only if each coordinate is taken from the same parent (which happens with probability 1/ρ n-1 ) or from a parent with a very similar value in this coordinate. The narrower the ellipsoid the more similar (i.e. correlated) the value needs to be. As another illustration consider sampling, neutral selection and discrete recombination based on Figure 1, right: after discrete recombination the points (-2, 2) and (2, -2) outside the ellipsoid have the same probability as the points (2, 2) and (-2, -2) inside the ellipsoid. The mutation operators introduced are unbiased in several ways. They are all point-symmetrical and have expectation zero. Therefore, mutation alone will 1.875 2.125 ) with the same eigenvalues (1/4, 4) as the diagonal matrix. Shown are in each subfigure the mean at 0 as small black dot (a different mean solely changes the axis annotations), two eigenvectors of C along the principal axes of the ellipsoids (thin black lines), two ellipsoids reflecting the set of points {x : (x -0) C -1 (x -0) ∈ {1, 4}} that represent the 1-σ and 2-σ lines of equal density, and 100 sampled points (however, a few of them are likely to be outside of the area shown). almost certainly not lead to better fitness values in expectation. The isotropic mutation operator features the same distribution along any direction. The general mutation operator is, as long as C remains unspecified, unbiased towards the choice of a Cartesian coordinate system, i.e. unbiased towards the representation of solutions x, which has also been referred to as invariance to affine coordinate system transformations [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. This however depends on the way how C is adapted (see below).

Parameter Control

Controlling the parameters of the mutation operator is key to the design of evolution strategies. Consider the isotropic operator (Figure 1, left), where the stepsize σ is a scaling factor for the random vector pertur-bation. The step-size controls to a large extent the convergence speed. In situations where larger stepsizes lead to larger expected improvements, a stepsize control technique should aim at increasing the step-size (and decreasing it in the opposite scenario).

The importance of step-size control is illustrated with a simple experiment. Consider a spherical function f (x) = x α , α > 0, and a (1+1)-ES with constant step-size equal to σ = 10 -2 , i.e. with mutations drawn from 10 -2 N (0, I). The convergence of the algorithm is depicted in Fig 2 (constant σ graphs).

We observe, roughly speaking, three stages: up to 600 function evaluations, progress towards the optimum is slow. In this stage the fixed step-size is too small. Between 700 and 800 evaluations, fast progress towards the optimum is observed. In this stage the step-size is close to optimal. Afterwards, the progress ), and of the (1 + 1)-ES with 1/5th success rule (Algorithm 1) on a spherical function f (x) = x α , α > 0 (because of invariance to monotonic f -transformation the same graph is observed for any α > 0). For each algorithm there are three runs in the left plot and three runs in the right plot. The x-axis is linear in the left and in log-scale in the right hand plot. For the (1+1)-ES with constant step-size, σ equals 10 -2 . For the (1+1)-ES with 1/5th success rule, the initial step-size is chosen very small to 10 -9 and the parameter d equals 1 + 10/3. On the left, also the evolution of the step-size of one of the runs of the (1+1)-ES with 1/5th success rule is shown. All algorithms are initialized at 1. Eventually, the (1+1)-ES with 1/5th success rule reveals linear behavior on the left, while the other two algorithms reveal eventually linear behavior in the right hand plot.

decreases and approaches the rate of the pure random search algorithm, well illustrated in the right subfigure. In this stage the fixed step-size is too large and the probability to sample better offspring becomes very small.

The figure also shows runs of the (1 + 1)-ES with 1/5th success rule step-size control (as described in Section 3.1) and the step-size evolution associated to one of these runs. The initial step-size is far too small and we observe that the adaptation technique increases the step-size in the first iterations. Afterwards, step-size is kept roughly proportional to the distance to the optimum, which is in fact optimal and leads to linear convergence in the left subfigure.

Generally, the goal of parameter control is to drive the endogenous strategy parameters close to their optimal values. These optimal values, as we have seen for the step-size in Figure 2, can significantly change over time or depending on the position in search space. In the most general case, the mutation operator has (n 2 + n)/2 degrees of freedom (see Section 2.4). The conjecture is that in the desired scenario lines of equal density of the mutation operator resemble locally the lines of equal fitness [4, p242f]. In case of convex-quadratic fitness functions this resemblance can be perfect and, apart from the step-size, optimal parameters do not change over time (as illustrated in Fig. 3

below).

Control parameters like the step-size can be stored on different "levels". Each individual can have its own step-size value (like in Algorithms 2 and 3), or a single step-size is stored and applied to all individuals in the population. In the latter case, sometimes different populations with different parameter values are run in parallel [START_REF] Herdy | The number of offspring as strategy parameter in hierarchically organized evolution strategies[END_REF].

In the following, six specific evolution strategies are outlined, each of them representing an important achievement in parameter control.

The 1/5th Success Rule

The 1/5th success rule for step-size control is based on an important discovery made very early in the research of evolution strategies [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF]. A similar rule had also been found independently before in [START_REF] Schumer | Adaptive step size random search[END_REF]. As a control mechanism in practice, the 1/5th success rule has been mostly superseded by more sophisticated methods. However, its conceptual insight remains remarkably valuable.

Consider a linear fitness function, for example f : x → x 1 or f : x → i x i . In this case, any point symmetrical mutation operator has a success probability of 1/2: in one half of the cases, the perturbation will improve the original solution, in one half of the cases the solution will deteriorate. Following from Taylors formula, smooth functions are locally linear, that is, they appear to be more and more linear with decreasing neighborhood size. Therefore, the success probability becomes 1/2 for step-size σ → 0. On most non-linear functions, the success rate is indeed a monotonously decreasing function in σ and goes to zero for σ → ∞. This suggests to control the step-size by increasing it for large success rates and decreasing it for small ones. This mechanism can drive the step-size close to the optimal value.

Rechenberg [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF] investigated two simple but quite different functions, the corridor function

f : x → x 1 if |x i | ≤ 1 for i = 2, . . . , n ∞ otherwise
and the sphere function f : x → x 2 i . He found optimal success rates for the (1+1)-ES with isotropic mutation to be ≈ 0.184 > 1/6 and ≈ 0.270 < 1/3, respectively (for n → ∞) [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF]. 6 This leads to approximately 1/5 as being the success value where to switch between decreasing and increasing the step-size.

Algorithm 1 implements the (1+1)-ES with 1/5th success rule in a simple and effective way [START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF]. Lines 4-6 implement Line 8 from Template 2, including selection in Line 7. Line 4 in Algorithm 1 updates the step-size σ of the single parent. The step-size does not change if and only if the argument of exp is zero.

Algorithm 1 The (1+1)-ES with 1/5th Rule 0 given n ∈ N + , d ≈ √ n + 1 1 initialize x ∈ R n , σ > 0 2 while not happy 3 x 1 = x + σ × N (0, I) // mutation 4 σ ← σ × exp 1/d (1 f (x1)≤f (x) -1/5) 5 if f (x 1 ) ≤ f (x) // select if better 6 x = x 1 // x-value of new parent
While this cannot happen in a single generation, we still can find a stationary point for σ: log σ is unbiased if and only if the expected value of the argument of exp is zero. This is the case if E1 f (x1)≤f (x) = 1/5, in other words, if the probability of an improvement with f (x 1 ) ≤ f (x) is 20%. Otherwise, log σ increases in expectation if the success probability is larger than 1/5 and decreases if the success probability is smaller than 1/5. Hence, Algorithm 1 indeed implements the 1/5th success rule.

Self-Adaptation

A seminal idea in the domain of evolution strategies is parameter control via self-adaptation [START_REF] Schwefel | Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie[END_REF]. In self-adaptation, new control parameter settings are generated similar to new x-vectors by recombination and mutation. Algorithm 2 presents an example with adaptation of n coordinate-wise standard deviations (individual step-sizes). First, for conducting the mutation, random events are drawn in Lines 4-6. In Line 7, the step-size vector for each individual undergoes (i) a mutation common for all components, exp(ξ k ), and (ii) a componentwise mutation with exp(ξ k ). These mutations are unbiased, in that E log σ k = log σ. The mutation of x in Line 8 uses the mutated vector σ k . After selection in Line 9, intermediate recombination is applied to compute x and σ for the next generation. By taking the average over σ k we have Eσ = Eσ k in Line 10. However, the application of mutation and recombination on σ introduces a moderate bias such that σ tends to increase under neutral selection [START_REF] Hansen | An analysis of mutative σ-selfadaptation on linear fitness functions[END_REF].

In order to achieve stable behavior of σ, the num-

Algorithm 2 The (µ/µ, λ)-σSA-ES 0 given n ∈ N + , λ ≥ 5n, µ ≈ λ/4 ∈ N, τ ≈ 1/ √ n, τ i ≈ 1/n 1/4 1 initialize x ∈ R n , σ ∈ R n + 2 while not happy 3 for k ∈ {1, . . . , λ} // random numbers i.i.d. for all k 4 ξ k = τ N (0, 1) // global step-size 5 ξ k = τ i N (0, I) // coordinate-wise σ 6 z k = N (0, I) //x-vector change // mutation 7 σ k = σ • exp(ξ k ) × exp(ξ k ) 8 x k = x + σ k • z k 9 P = sel µ best({(x k , σ k , f (x k )) | 1 ≤ k ≤ λ}) // recombination 10 σ = 1 µ σ k ∈P σ k 11 x = 1 µ x k ∈P
x k ber of parents µ must be large enough, which is reflected in the setting of λ. A setting of τ ≈ 1/4 has been proposed in combination with ξ k being uniformly distributed across the two values in {-1, 1}

.

Derandomized Self-Adaptation

Derandomized self-adaptation [START_REF] Ostermeier | A derandomized approach to self-adaptation of evolution strategies[END_REF] addresses the problem of selection noise that occurs with selfadaptation of σ as outlined in Algorithm 2. Selection noise refers to the possibility that very good offspring may be generated with poor strategy parameter settings and vice versa. The problem occurs frequently and has two origins. For µ n the stochastic fluctuations become prohibitive and therefore µ ≈ λ/4 ≥ 1.25n is chosen to make σ-self-adaptation reliable.

• A small/large component in |σ k • z k | (
Derandomization addresses the problem of selection noise on σ directly without resorting to a large parent number. The derandomized (1, λ)-σSA-ES is outlined in Algorithm 3 and addresses selection noise twofold. Instead of introducing new variations in σ by means of exp(ξ k ), the variations from z k are directly used for the mutation of σ in Line 7. The variations are dampened compared to their use in the mutation of x (Line 6) via d and d i , thereby mimicking the effect of intermediate recombination on σ [START_REF] Ostermeier | A derandomized approach to self-adaptation of evolution strategies[END_REF][START_REF] Runarsson | Reducing random fluctuations in mutative self-adaptation, Parallel Problem Solving from Nature[END_REF]. The order of the two mutation equations becomes irrelevant.

For Algorithm 3 also a (µ/µ, λ) variant with recombination is feasible. However, in particular in the (µ/µ I , λ)-ES, σ-self-adaptation tends to generate too small step-sizes. A remedy for this problem is to use non-local information for step-size control.

Non-Local Derandomized Step-Size Control (CSA)

When using self-adaptation, step-sizes are associated with individuals and selected based on the fitness of each individual. However, step-sizes that serve individuals well by giving them a high likelihood to be selected are generally not step-sizes that maximize the progress of the entire population or the parental

Algorithm 3 Derandomized (1, λ)-σSA-ES 0 given n ∈ N + , λ ≈ 10, τ ≈ 1/3, d ≈ √ n, d i ≈ n 1 initialize x ∈ R n , σ ∈ R n + 2 while not happy 3 for k ∈ {1, . . . , λ} // random numbers i.i.d. for all k 4 ξ k = τ N (0, 1) 5 z k = N (0, I)
// mutation, re-using random events 6

x k = x + exp(ξ k ) × σ • z k 7 σ k = σ • exp 1/di |z k | E|N (0, 1)| -1 7b × exp 1/d (ξ k ) 8 (x 1 , σ 1 , f (x 1 )) ← select single best( {(x k , σ k , f (x k )) | 1 ≤ k ≤ λ}) // assign new parent 9 σ = σ 1 10 x = x 1 centroid x.
We will see later that, for example, the optimal step-size may increase linearly with µ (Section 4.2 and Eq. ( 16)). With self-adaptation on the other hand, the step-size of the µth best offspring is typically even smaller than the step-size of the best offspring. Consequently, Algorithm 3 assumes often too small step-sizes and can be considerably improved by using non-local information about the evolution of the population. Instead of single (local) mutation steps z, an exponentially fading record, s σ , of mutation steps is taken. This record, referred to as search path or evolution path, can be pictured as a sequence or sum of consecutive successful z-steps that is non-local in time and space. A search path carries information about the interrelation between single steps. This information can improve the adaptation and search procedure remarkably. Algorithm 4 outlines the (µ/µ I , λ)-ES with cumulative path length control, also denoted as cumulative step-size adaptation (CSA), and additionally with non-local individual step-size adaptation [START_REF] Ostermeier | Stepsize adaptation based on non-local use of selection information, Parallel Problem Solving from Nature[END_REF][START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF].

In the (µ/µ, λ)-ES with search path, Algorithm 4, the factor ξ k for changing the overall step-size has Algorithm 4 The (µ/µ, λ)-ES with Search Path

0 given n ∈ N + , λ ∈ N, µ ≈ λ/4 ∈ N, c σ ≈ µ/(n + µ), d ≈ 1 + µ/n, d i ≈ 3n 1 initialize x ∈ R n , σ ∈ R n + , s σ = 0 2 while not happy 3 for k ∈ {1, . . . , λ} 4 z k = N (0, I) // i.i.d. for each k 5 x k = x + σ • z k 6 P ← sel µ best({(x k , z k , f (x k )) | 1 ≤ k ≤ λ})
// recombination and parent update 7 Using a search path is justified in two ways. First, it implements a low-pass filter for selected z-steps, removing high frequency (most likely noisy) information. Second, and more importantly, it utilizes information that is otherwise lost: even if all single steps have the same length, the length of s σ can vary, because it depends on the correlation between the directions of z-steps. If single steps point into similar directions, the path will be up to almost 2/c σ times longer than a single step and the step-size will increase. If they oppose each other the path will be up to almost c σ /2 times shorter and the step-size will decrease. The same is true for single components of s σ .

s σ ← (1 -c σ ) s σ + 7b c σ (2 -c σ ) √ µ µ z k ∈P z k 8 σ ← σ • exp 1/di |s σ | E|N (0, 1)| -1 8b × exp cσ/d s σ E N (0, I) -1 9 x = 1 µ x k ∈P
The factors c σ (2 -c σ ) and √ µ in Line 7b guaranty unbiasedness of s σ under neutral selection, as usual.

All evolution strategies described so far are of somewhat limited value, because they feature only isotropic or axis-parallel mutation operators. In the remainder we consider methods that entertain not only an n-dimensional step-size vector σ, but also correlations between variables for the mutation of x.

Addressing Dependencies Between Variables

The evolution strategies presented so far sample the mutation distribution independently in each component of the given coordinate system. The lines of equal density are either spherical or axis-parallel ellipsoids (compare Figure 1). This is a major drawback, because it allows to solve problems with a long or elongated valley efficiently only if the valley is aligned with the coordinate system. In this section we discuss evolution strategies that allow to traverse non-axis-parallel valleys efficiently by sampling distributions with correlations.

Full Covariance Matrix Algorithms that adapt the complete covariance matrix of the mutation distribution (compare Section 2.4) are correlated mutations [START_REF] Schwefel | Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie[END_REF], the generating set adaptation [START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF], the covariance matrix adaptation (CMA) [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF], a mutative invariant adaptation [START_REF] Ostermeier | An evolution strategy with coordinate system invariant adaptation of arbitrary normal mutation distributions within the concept of mutative strategy parameter control[END_REF], and some instances of natural evolution strategies [START_REF] Wierstra | Natural evolution strategies[END_REF][START_REF] Sun | Schmidhuber: Efficient natural evolution strategies[END_REF][START_REF] Glasmachers | Exponential natural evolution strategies[END_REF]. Correlated mutations and some natural evolution strategies are however not invariant under changes of the coordinate system [START_REF] Hansen | Invariance, self-adaptation and correlated mutations and evolution strategies, Parallel Problem Solving from Nature[END_REF][START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Glasmachers | Exponential natural evolution strategies[END_REF]. In the next sections we outline two evolution strategies that adapt the full covariance matrix reliably and are invariant under coordinate system changes: the covariance matrix adaptation evolution strategy (CMA-ES) and the exponential natural evolution strategy (xNES).

Restricted Covariance Matrix Algorithms that adapt non-diagonal covariance matrices, but are restricted to certain matrices, are the momentum adaptation [START_REF] Ostermeier | An evolution strategy with momentum adaptation of the random number distribution, Parallel Problem Solving from Nature (PPSN II)[END_REF], direction adaptation [START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF], main vector adaptation [START_REF] Poland | Main vector adaptation: A CMA variant with linear time and space complexity[END_REF], and limited memory CMA-ES [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF]. These variants are limited in their capability to shape the mutation distribution, but they might be advantageous for larger dimensional problems, say larger than a hundred.

Covariance Matrix Adaptation (CMA)

The covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF][START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Hansen | Evaluating the CMA evolution strategy on multimodal test functions, Parallel Problem Solving from Nature[END_REF] is a de facto standard in continuous domain evolutionary computation. The CMA-ES is a natural generalization of Algorithm 4 in that the mutation ellipsoids are not constrained to be axis-parallel, but can take on a general orientation.

The CMA-ES is also a direct successor of the generating set adaptation [START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF], replacing self-adaptation to control the overall step-size with cumulative step-size adaptation.

The (µ/µ W , λ)-CMA-ES is outlined in Algorithm 5. Two search paths are maintained, s σ and s c . The first path, s σ , accumulates steps in the coordinate system where the mutation distribution is isotropic and which can be derived by scaling in the principal axes of the mutation ellipsoid only. The path generalizes s σ from Algorithm 4 to non-diagonal covariance matrices and is used to implement cumulative step-size adaptation, CSA, in Line 10 (resembling Line 8b in Algorithm 4). Under neutral selection, s σ ∼ N (0, I) and log σ is unbiased.

The second path, s c , accumulates steps, disregarding σ, in the given coordinate system. 7 The covariance matrix update consists of a rank-one update, based on the search path s c , and a rank-µ update with µ nonzero recombination weights w k . Under neutral selection the expected covariance matrix equals the covariance matrix before the update.

The updates of x and C follow a common principle. The mean x is updated such that the likelihood of successful offspring to be sampled again is maximized (or increased if c m < 1). The covariance matrix C is updated such that the likelihood of successful steps (x k -x)/σ to appear again, or the likelihood Algorithm 5 The (µ/µ W , λ)-CMA-ES

0 given n ∈ N + , λ ≥ 5, µ ≈ λ/2, w k = w (k)/ µ k w (k), w (k) = log(λ/2 + 1/2) - log rank(f (x k )), µ w = 1/ µ k w 2 k , c σ ≈ µ w /(n + µ w ), d ≈ 1 + µ w /n, c c ≈ (4 + µ w /n)/(n + 4 + 2µ w /n), c 1 ≈ 2/(n 2 + µ w ), c µ ≈ µ w /(n 2 + µ w ), c m = 1 1 initialize s σ = 0, s c = 0, C = I, σ ∈ R n + , x ∈ R n 2 while not happy 3 for k ∈ {1, . . . , λ} 4 z k = N (0, I) // i.i.d. for all k 5 x k = x + σC 1 /2 × z k 6 P = sel µ best({(z k , f (x k )) | 1 ≤ k ≤ λ}) 7 x ← x + c m σ C 1 /2 z k ∈P w k z k 8 s σ ← (1 -c σ ) s σ + // search path for σ c σ (2 -c σ ) √ µ w z k ∈P w k z k 9 s c ← (1 -c c ) s c + // search path for C h σ c c (2 -c c ) √ µ w z k ∈P w k C 1 /2 z k 10 σ ← σ exp cσ/d s σ E N (0, I) -1 11 C ← (1 -c 1 + c h -c µ ) C + c 1 s c s c + c µ z k ∈P w k C 1 /2 z k (C 1 /2 z k )
where

h σ = 1 sσ 2 /n<2+4/(n+1) , c h = c 1 (1 - h 2 σ )c c (2 -c c ), and C 1 /2 is the unique symmetric positive definite matrix obeying C 1 /2 × C 1 /2 = C. All c-coefficients are ≤ 1.
to sample (in direction of) the path s c , is increased. A more fundamental principle for the equations is given in the next section.

Using not only the µ best but all λ offspring can be particularly useful for the "rank-µ" update of C in line 11 where negative weights w k for inferior offspring are advisable. Such an update has been introduced as active CMA [START_REF] Jastrebski | Improving evolution strategies through active covariance matrix adaptation[END_REF].

The factor c m in Line 7 can be equally written as a mutation scaling factor κ = 1/c m in Line 5, compare [START_REF] Beyer | Mutate large, but inherit small! On the analysis of rescaled mutations in ( 1, λ)-ES with noisy fitness data, Parallel Problem Solving from Nature[END_REF]. This means that the actual mutation steps are larger than the inherited ones, resembling the derandomization technique of damping step-size changes to address selection noise as described in Section 3.3.

An elegant way to replace Line 10 is

σ ← σ exp (cσ/d)/2 s σ 2 n -1 (1)
and often used in theoretical investigations of this update as those presented in Section 4.2.

A single run of the (5/5 W , 10)-CMA-ES on a convex-quadratic function is shown in Fig. 3. For sake of demonstration the initial step-size is chosen far too small (a situation that should be avoided in practice) and increases quickly for the first 400 fevaluations. After no more than 5500 f -evaluations the adaptation of C is accomplished. Then the eigenvalues of C (square roots of which are shown in the lower left) reflect the underlying convex-quadratic function and the convergence speed is the same as on the sphere function and about 60% of the speed of the (1+1)-ES as observed in Fig. 2. The resulting convergence speed is about ten thousand times faster than without adaptation of C and at least one thousand times faster compared to any of the algorithms from the previous sections.

Natural Evolution Strategies

The idea of using natural gradient learning [START_REF] Amari | Natural gradient works efficiently in learning[END_REF] in evolution strategies has been proposed in [START_REF] Wierstra | Natural evolution strategies[END_REF] and further pursued in [START_REF] Sun | Schmidhuber: Stochastic search using the natural gradient[END_REF][START_REF] Glasmachers | Exponential natural evolution strategies[END_REF]. Natural evolution strategies (NES) put forward the idea that the update of all distribution parameters can be based on the same fundamental principle. NES have been proposed as a more principled alternative to CMA-ES and characterized by operating on Cholesky factors of a covariance matrix. Only later was it discovered that also CMA-ES implements the underlying NES principle of natural gradient learning [START_REF] Akimoto | Bidirectional relation between CMA evolution strategies and natural evolution strategies, Parallel Problem Solving from Nature[END_REF][START_REF] Glasmachers | Exponential natural evolution strategies[END_REF].

For simplicity, let the vector θ represent all parameters of the distribution to sample new offspring. In case of a multivariate normal distribution as above, we have a bijective transformation between θ and mean and covariance matrix of the distribution, θ ↔ (x, σ 2 C). Standard Deviations in All Coordinates i with α i = 10 3(i-1)/(n-1) , y = Rx, where R is a random matrix with R R = I, for n = 10. Shown is the evolution of various parameters against the number of function evaluations. Upper left: best (thick blue line), median and worst fitness value that reveal the final convergence phase after about 5500 function evaluations where the ellipsoid function has been reduced to the simple sphere; minimal and maximal coordinate-wise standard deviation of the mutation distribution and in between (mostly hidden) the step-size σ that is initialized far too small and increases quickly in the beginning, that increases afterwards several times again by up to one order of magnitude and decreases with maximal rate during the last 1000 f -evaluations; axis ratio of the mutation ellipsoid (square root of the condition number of C) that increases from 1 to 1000 where the latter corresponds to α n /α 1 . Lower left: sorted principal axis lengths of the mutation ellipsoid disregarding σ (square roots of the sorted eigenvalues of C, see also Fig. 1) that adapt to the (local) structure of the underlying optimization problem; they finally reflect almost perfectly the factors α -1 i up to a constant factor. Upper right: x (distribution mean) that is initialized with all ones and converges to the global optimum in zero while correlated movements of the variables can be observed. Lower right: standard deviations in the coordinates disregarding σ (square roots of diagonal elements of C) showing the R-dependent projections of the principal axis lengths into the given coordinate system. The straight lines to the right of the vertical line at about 6300 only annotate the coordinates and do not reflect measured data.

We consider a probability density p(.|θ) over R n parametrized by θ, a non-increasing function W f θ : R → R, 8 and the expected W f θ -transformed fitness [START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]

J(θ) = E(W f θ (f (x))) x ∼ p(.|θ) = R n W f θ (f (x)) p(x|θ)dx , (2)
where the expectation is taken under the given sample distribution. The maximizer of J w.r.t. p(.|θ) is, for any fixed W f θ , a Dirac distribution concentrated on the minimizer of f . A natural way to update θ is therefore a gradient ascent step in ∇ θ J direction. However, the "vanilla" gradient ∇ θ J depends on the specific parametrization chosen in θ. In contrast, the natural gradient, denoted by ∇ θ , is associated to the Fisher metric that is intrinsic to p and independent of the chosen θ-parametrization. Developing ∇ θ J(θ) under mild assumptions on f and p(.|θ) by exchanging differentiation and integration, recognizing that the gradient ∇ θ does not act on W f θ , using the loglikelihood trick ∇ θ p(.|θ) = p(.|θ) ∇ θ ln p(.|θ) and finally setting θ = θ yields 9

∇ θ J(θ) = E W f θ (f (x)) ∇ θ ln p(x|θ) . (3)
A Monte-Carlo approximation of the expected value by the average finally yields the comparatively simple expression

∇ θ J(θ) ≈ 1 λ λ k=1 preference weight W f θ (f (x k )) ∇ θ ln p(x k |θ) intrinsic candidate direction (4)
8 More specifically, W f θ : y → w(Pr z∼p(.|θ) (f (z) ≤ y)) computes the p θ -quantile, or cumulative distribution function, of f (z) with z ∼ p(.|θ) at point y, composed with a nonincreasing predefined weight function w : [0, 1] → R (where w(0) > w(1/2) = 0 is advisable). The value of W f θ (f (x)) is invariant under strictly monotonous transformations of f . For x ∼ p(.|θ) the distribution of W f θ (f (x)) ∼ w(U [0, 1]) depends only on the predefined w; it is independent of θ and f and therefore also (time-)invariant under θ-updates. Given λ samples x k , we have the rank-based consistent estimator

W f θ (f (x k )) ≈ w rank(f (x k ))-1/2 λ
. 9 We set θ = θ because we will estimate W θ using the current samples that are distributed according to p(.|θ) for a natural gradient update of θ, where x k ∼ p(.|θ) is sampled from the current distribution. The natural gradient can be computed as ∇ θ = F -1 θ ∇ θ , where F θ is the Fisher information matrix expressed in θcoordinates. For the multivariate Gaussian distribution, ∇ θ ln p(x k |θ) can indeed be easily expressed and computed efficiently. We find that in CMA-ES (Algorithm 5), the rank-µ update (Line 11 with c 1 = 0) and the update in Line 7 are natural gradient updates of C and x, respectively [START_REF] Akimoto | Bidirectional relation between CMA evolution strategies and natural evolution strategies, Parallel Problem Solving from Nature[END_REF][START_REF] Glasmachers | Exponential natural evolution strategies[END_REF], where the kth largest w k is a consistent estimator for the kth largest W f θ (f (x k )) [START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF].

While the natural gradient does not depend on the parametrization of the distribution, a finite step taken in the natural gradient direction does. This becomes relevant for the covariance matrix update, where natural evolution strategies take a different parametrization than CMA-ES. Starting from Line 11 in Algorithm 5, we find for

c 1 = c h = 0 C ← (1 -c µ ) C + c µ z k ∈P w k C 1 /2 z k (C 1 /2 z k ) = C 1 /2 (1 -c µ ) I + c µ z k ∈P w k z k z k C 1 /2 w k =1 = C 1 /2 I + c µ z k ∈P w k z k z k -I C 1 /2 cµ 1 ≈ C 1 /2 exp cµ z k ∈P w k z k z k -I C 1 /2 . ( 5
)
The term bracketed between the matrices C 1 /2 in the lower three lines is a multiplicative covariance matrix update expressed in the natural coordinates, where the covariance matrix is the identity and C

1 /2 serves as coordinate system transformation into the given coordinate system. Only the lower two lines of Eq. ( 5) do not rely on the constraint k w k = 1 in order to satisfy a stationarity condition on C. 10 The last line of Eq. ( 5) is used in the exponential natural evolution Algorithm 6 The Exponential NES (xNES) [START_REF] Glasmachers | Exponential natural evolution strategies[END_REF] and guarantees positive definiteness of C even with negative weights, independent of c µ and of the data z k . The xNES is depicted in Algorithm 6.

0 given n ∈ N + , λ ≥ 5, w k = w (k)/ k |w (k)|, w (k) ≈ log(λ/2 + 1/2) -log rank(f (x k )), η c ≈ (5 + λ)/(5 n 1.5 ) ≤ 1, η σ ≈ η c , η x ≈ 1 1 initialize C 1 /2 = I, σ ∈ R + , x ∈ R n 2 while not happy 3 for k ∈ {1, . . . , λ} 4 z k = N (0, I) // i.i.d. for all k 5 x k = x + σC 1 /2 × z k 6 P = {(z k , f (x k )) | 1 ≤ k ≤ λ} 7 x ← x + η x σ C 1 /2 z k ∈P w k z k 8 σ ← σ exp ησ/2 z k ∈P w k z k 2 n -1 9 C 1 /2 ← C 1 /2 × exp ηc/2 z k ∈P w k z k z k - z k 2 n I strategy, xNES
In xNES, sampling is identical to CMA-ES and environmental selection is omitted entirely. Line 8 resembles the step-size update in [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF]. Comparing the updates more closely, with c σ = 1 Eq. ( 1) uses

µ w k w k z k 2 n -1 whereas xNES uses k w k z k 2 n - 1
for updating σ. For µ = 1 the updates are the same. For µ > 1, the latter only depends on the lengths of the z k , while the former depends on their lengths and directions. Finally, xNES expresses the update Eq. ( 5) in Line 9 on the Cholesky factor C 1 /2 , which does not remain symmetric in this case

(C = C 1 /2 × C 1 /2
still holds). The term -z k 2 /n keeps the determinant of C 1 /2 (and thus the trace of log C

1 /2 ) constant and is of rather cosmetic nature. Omitting the term is equivalent to using η σ + η c instead of η σ in line 8.

The exponential natural evolution strategy is a very elegant algorithm. Like CMA-ES it can be interpreted as an incremental Estimation of Distribution Algorithm [START_REF] Pelikan | Introduction to estimation of distribution algorithms[END_REF]. However, it performs generally inferior compared to CMA-ES because it does not use search paths for updating σ and C.

Further Aspects

Internal Parameters Adaptation and selfadaptation address the control of the most important internal parameters in evolution strategies. Yet, all algorithms presented have hidden and exposed parameters in their implementation. Many of them can be set to reasonable and robust default values. The population size parameters µ and λ however change the search characteristics of an evolution strategy significantly. Larger values, in particular for parent number µ, often help address highly multimodal or noisy problems more successfully.

In practice, several experiments or restarts are advisable, where different initial conditions for x and σ can be employed. For exploring different population sizes, a schedule with increasing population size, IPOP, is advantageous [START_REF] Harik | A parameter-less genetic algorithm[END_REF][START_REF] Lobo | The parameterless genetic algorithm in practice[END_REF][START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF], because runs with larger populations take typically more function evaluations. Preceding long runs (large µ and λ) with short runs (small µ and λ) leads to a smaller (relative) impairment of the later runs than vice versa.

Internal computational complexity Algorithms presented in Sections 3.1-3.4 that sample isotropic or axis-parallel mutation distributions have an internal computational complexity linear in the dimension. The internal computational complexity of CMA-ES and xNES is, for constant population size, cubic in the dimension due to the update of C 1 /2 . Typical implementations of the CMA-ES however have quadratic complexity, as they implement a lazy update scheme for C 1 /2 , where C is decomposed into C 1 /2 C 1 /2 only after about n/λ iterations. An exact quadratic update for CMA-ES has also been proposed [START_REF] Suttorp | Efficient covariance matrix update for variable metric evolution strategies[END_REF]. While never considered in the literature, a lazy update for xNES to achieve quadratic complexity seems feasible as well.

Invariance Selection and recombination in evolution strategies are based solely on the ranks of offspring and parent individuals. As a consequence, the behavior of evolution strategies is invariant under order-preserving (strictly monotonous) transformations of the fitness function value. In particular, all spherical unimodal functions belong to the same function class, which the convex-quadratic sphere function is the most pronounced member of. This function is more thoroughly investigated in Section 4.

All algorithms presented are invariant under translations and Algorithms 1, 5 and 6 are invariant under rotations of the coordinate system, provided that the initial x is translated and rotated accordingly.

Parameter control can introduce yet further invariances. All algorithms presented are scale invariant due to step-size adaptation. Furthermore, ellipsoidal functions that are in the reach of the mutation operator of the evolution strategies presented in Sections 3.2 to 3.7 are eventually transformed, effectively, into spherical functions. These evolution strategies are invariant under the respective affine transformations of the search space, given the initial conditions are chosen respectively.

Variants Evolution strategies have been extended and combined with other approaches in various ways. We mention here constraint handling [START_REF] Michalewicz | Evolutionary algorithms for constrained parameter optimizationproblems[END_REF][START_REF] Mezura-Montes | Constraint-handling in nature-inspired numerical optimization: Past, present, and future[END_REF], fitness surrogates [START_REF] Emmerich | Metamodel-assisted evolution strategies, Parallel Problem Solving from Nature[END_REF], multi-objective variants [START_REF] Igel | Covariance matrix adaptation for multi-objective optimization[END_REF][START_REF] Hansen | Improved step size adaptation for the MO-CMA-ES[END_REF], and exploitation of fitness values [START_REF] Salomon | Evolutionary algorithms and gradient search: Similarities and differences[END_REF].

Theory

There is ample empirical evidence, that on many unimodal functions evolution strategies with step-size control, as those outlined in the previous section, converge fast and with probability one to the global optimum. Convergence proofs supporting this evidence are discussed in Section 4.3. On multimodal functions on the other hand, the probability to converge to the global optimum (in a single run of the same strategy) is generally smaller than one (but larger than zero), as suggested by observations and theoretical results [START_REF] Rudolph | Self-adaptive mutations may lead to premature convergence[END_REF]. Without parameter control on the other hand, elitist strategies always converge to the essential global optimum, 11 however at a much slower rate (compare random search in Figure 2).

In this section we use a time index t to denote iteration and assume, for notational convenience and without loss of generality (due to translation invariance), that the optimum of f is in x * = 0. This simplifies writing x (t) -x * to simply x (t) and then x (t) measures the distance to the optimum of the parental centroid in time step t.

Linear convergence plays a central role for evolution strategies. For a deterministic sequence x (t) linear convergence (towards zero) takes place if there exists a c > 0 such that

lim t→∞ x (t+1) x (t) = exp(-c) (6)
which means, loosely speaking, that for t large enough, the distance to the optimum decreases in every step by the constant factor exp(-c). Taking the logarithm of Eq. ( 6), then exchanging the logarithm and the limit and taking the Cesàro mean yields lim

T →∞ 1 T T -1 t=0 log x (t+1) x (t) = 1 T log x (T ) / x (0) = -c . ( 7
)
For a sequence of random vectors we define linear convergence based on Eq. ( 7) as follows.

Definition 1 (linear convergence). The sequence of random vectors x (t) converges almost surely linearly to 0 if there exists a c > 0 such that

-c = lim T →∞ 1 T log x (T )
x (0) a.s.

= lim

T →∞ 1 T T -1 t=0 log x (t+1) x (t) a.s. ( 8
)
The sequence converges in expectation linearly to 0 if there exists a c > 0 such that

-c = lim t→∞ E log x (t+1) x (t) . ( 9
)
The constant c is the convergence rate of the algorithm.

Linear convergence hence means that asymptotically in t, the logarithm of the distance to the optimum decreases linearly in t like -ct. This behavior has been observed in Figure 2 for the (1+1)-ES with 1/5th success rule on a unimodal spherical function.

Note that λ function evaluations are performed per iteration and it is then often useful to consider a convergence rate per function evaluation, i.e. to normalize the convergence rate by λ.

The progress rate measures the reduction of the distance to optimum within a single generation [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF].

Definition 2 (progress rate). The normalized progress rate is defined as the expected relative reduction of x (t)

ϕ * = n E x (t) -x (t+1) x (t) x (t) , s (t) = n 1 -E x (t+1) x (t) x (t) , s (t) , ( 10
)
where the expectation is taken over x (t+1) , given (x (t) , s (t) ). In situations commonly considered in theoretical analyses, ϕ * does not depend on x (t) and is expressed as a function of strategy parameters s (t) .

Definitions 1 and 2 are related, in that for a given

x (t) ϕ * ≤ -n log E x (t+1) x (t) (11) ≤ -n E log x (t+1) x (t) = nc . ( 12
)
Therefore, progress rate ϕ * and convergence rate nc do not agree and we might observe convergence (c > 0) while ϕ * < 0. However for n → ∞, we typically have ϕ * = nc [START_REF] Auger | Reconsidering the progress rate theory for evolution strategies in finite dimensions[END_REF].

The normalized progress rate ϕ * for evolution strategies has been extensively studied in various situations, see Section 4.2. Scale-invariance and (sometimes artificial) assumptions on the step-size typically ensure that the progress rates do not depend on t.

Another way to describe how fast an algorithm approaches the optimum is to count the number of function evaluations needed to reduce the distance to the optimum by a given factor 1/ or, similarly, the runtime to hit a ball of radius around the optimum, starting, e.g., from distance one. Definition 3 (runtime). The runtime is the first hitting time of a ball around the optimum. Specifically, the runtime in number of function evaluations as a function of reads

λ × min t : x (t) ≤ × x (0) = λ × min t : x (t) x (0) ≤ . (13)
Linear convergence with rate c as given in Eq. ( 9) implies that, for → 0, the expected runtime divided by log(1/ ) goes to the constant λ/c.

Lower Runtime Bounds

Evolution strategies with a fixed number of parent and offspring individuals cannot converge faster than linearly and with a convergence rate of O(1/n). This means that their runtime is lower bounded by a constant times log(1/ n ) = n log(1/ ) [START_REF] Teytaud | General lower bounds for evolutionary algorithms, Parallel Problem Solving from Nature[END_REF][START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF][START_REF] Teytaud | Lower bounds for evolution strategies[END_REF][START_REF] Jägersküpper | Lower bounds for hit-andrun direct search[END_REF][START_REF] Jägersküpper | Lower bounds for randomized direct search with isotropic sampling[END_REF]. This result can be obtained by analyzing the branching factor of the tree of possible paths the algorithm can take. It therefore holds for any optimization algorithm taking decisions based solely on a bounded number of comparisons between fitness values [START_REF] Teytaud | General lower bounds for evolutionary algorithms, Parallel Problem Solving from Nature[END_REF][START_REF] Fournier | Lower bounds for comparison based evolution strategies using VC-dimension and sign patterns[END_REF][START_REF] Teytaud | Lower bounds for evolution strategies[END_REF]. More specifically, the runtime of any (1 + , λ)-ES with isotropic mutations cannot be asymptotically faster than ∝ n log(1/ ) λ/ log(λ) [START_REF]Jägersküpper: Probabilistic runtime analysis of (1 + ,λ) evolution strategies using isotropic mutations[END_REF]. Considering more restrictive classes of algorithms can provide more precise non-asymptotic bounds [START_REF] Jägersküpper | Lower bounds for hit-andrun direct search[END_REF][START_REF] Jägersküpper | Lower bounds for randomized direct search with isotropic sampling[END_REF]. Different approaches address in particular the (1+1)and (1, λ)-ES and precisely characterize the fastest convergence rate that can be obtained with isotropic normal distributions on any objective function with any step-size adaptation mechanism [START_REF] Jebalia | Log-linear convergence and optimal bounds for the (1+1)-ES, Evolution Artificielle[END_REF][START_REF] Auger | Reconsidering the progress rate theory for evolution strategies in finite dimensions[END_REF][START_REF] Auger | Theory of evolution strategies: A new perspective[END_REF][START_REF] Auger | Analyzing the impact of mirrored sampling and sequential selection in elitist evolution strategies[END_REF].

Considering the sphere function, the optimal convergence rate is attained with distance proportional step-size, that is, a step-size proportional to the distance of the parental centroid to the optimum, σ = const × x = σ * x /n. Optimal step-size and optimal convergence rate according to Eqs. ( 8) and ( 9) can be expressed in terms of expectation of some random variables that are easily simulated numerically. The convergence rate of the (1+1)-ES with distance proportional step-size is shown in Figure 4 as a function of the normalized step-size σ * = nσ/ x . The peak of each curve is the upper bound for the convergence rate that can be achieved on any function with any form of step-size adaptation. As for the general bound, the evolution strategy converges linearly and the convergence rate c decreases to zero like 1/n for n → ∞ [55, [START_REF] Auger | Mirrored sampling in evolution strategies with weighted recombination[END_REF][START_REF] Auger | Analyzing the impact of mirrored sampling and sequential selection in elitist evolution strategies[END_REF], which is equivalent to linear scaling of the runtime in the dimension. The asymptotic limit for the convergence rate of the (1+1)-ES, as shown in the lowest curve in Figure 4, coincides with the progress rate expression given in the next section.

Progress Rates

This section presents analytical approximations to progress rates of evolution strategies for sphere, ridge, and cigar functions in the limit n → ∞. Both onegeneration results and those that consider multiple time steps and cumulative step-size adaptation are considered.

The first analytical progress rate results date back to the early work of Rechenberg [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF] and Schwefel [START_REF] Schwefel | Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie[END_REF], who considered the sphere and corridor models and very simple strategy variants. Further results have since been derived for various ridge functions, several classes of convex quadratic functions, and more general constrained linear problems. The strategies that results are available for have increased in complexity as well and today include multi-parent strategies employing recombination as well as several step-size adaptation mechanisms. Only strategy variants with isotropic mutation distributions have been considered up to this point. However, parameter control strate- Figure 4: Normalized convergence rate nc versus normalized step-size nσ/ x of the (1+1)-ES with distance proportional step-size for n = 2, 3, 5, 10, 20, ∞ (top to bottom). The peaks of the graphs represent the upper bound for the convergence rate of the (1+1)-ES with isotropic mutation (corresponding to the lower runtime bound). The limit curve for n to infinity (lower black curve) reveals the optimal normalized progress rate of ϕ * ≈ 0.202 of the (1+1)-ES on sphere functions for n → ∞.

gies that successfully adapt the shape of the mutation distribution (such as CMA-ES) effectively transform ellipsoidal functions into (almost) spherical ones, thus lending extra relevance to the analysis of sphere and sphere-like functions.

The simplest convex quadratic functions to be optimized are variants of the sphere function (see also the discussion of invariance in Section 3.8)

f (x) = x 2 = n i=1 x 2 i = R 2 ,
where R denotes the distance from the optimal solution. Expressions for the progress rate of evolution strategies on sphere functions can be computed by decomposing mutation vectors into two components z and z as illustrated in Fig. 5. Component z is the projection of z onto the negative of the gradient vec-

(1+1)-ES on Sphere Functions

The normalized progress rate of the (1 + 1)-ES on sphere functions is

ϕ * = σ * √ 2π e -1 8 σ * 2 - σ * 2 4 1 -erf σ * √ 8 ( 14
)
in the limit of n → ∞ [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF]. The expression in square brackets is the success probability (i.e., the probability that the offspring candidate solution is superior to its parent and thus replaces it). The first term in Eq. ( 14) is the contribution to the normalized progress rate from the component z of the mutation vector that is parallel to the gradient vector.

The second term results from the component z that is perpendicular to the gradient direction. The black curve in Figure 4 illustrates how the normalized progress rate of the (1+1)-ES on sphere functions in the limit n → ∞ depends on the normalized mutation strength. For small normalized mutation strengths, the normalized progress rate is small as the short steps that are made do not yield significant progress. The success probability is nearly one half. For large normalized mutation strengths, progress is near zero as the overwhelming majority of steps result in poor offspring that are rejected. The normalized progress rate assumes a maximum value of ϕ * = 0.202 at normalized mutation strength σ * = 1.224. The range of step-sizes for which close to optimal progress is achieved is referred to as the evolution window [START_REF] Rechenberg | Evolutionstrategie: Optimierung technischer Systeme nach Prinzipien der biologischen Evolution[END_REF]. In the runs of the (1+1)-ES with constant step-size shown in Fig. 2, the normalized step-size initially is to the left of the evolution window (large relative distance to the optimal solution) and in the end to its right (small relative distance to the optimal solution), achieving maximal progress at a point in between.

(µ/µ, λ)-ES on Sphere Functions

The normalized progress rate of the (µ/µ, λ)-ES on sphere functions is described by

ϕ * = σ * c µ/µ,λ - σ * 2 2µ ( 15
)
in the limit n → ∞ [START_REF] Rechenberg | Evolutionsstrategie '94[END_REF]. The term c µ/µ,λ is the expected value of the average of the µ largest order statistics of λ independent standard normally distributed random numbers. For λ fixed, c µ/µ,λ decreases with increasing µ. For fixed truncation ratio µ/λ, c µ/µ,λ approaches a finite limit value as λ and µ increase [START_REF] Beyer | Toward a theory of evolution strategies: On the benefits of sex -The (µ/µ, λ) theory[END_REF][START_REF] Beyer | The Theory of Evolution Strategies[END_REF].

It is easily seen from Eq. ( 15) that the normalized progress rate of the (µ/µ, λ)-ES is maximized by normalized mutation strength

σ * = µc µ/µ,λ . (16)
The normalized progress rate achieved with that setting is

ϕ * = µc 2 µ/µ,λ 2 . ( 17
)
The progress rate is negative if σ * > 2µc µ/µ,λ . Figure 6 illustrates how the optimal normalized progress rate per offspring depends on the population size parameters µ and λ. Two interesting observations can be made from the figure.

• For all but the smallest values of λ, the (µ/µ, λ)-ES with µ > 1 is capable of significantly more rapid progress per offspring than the (1, λ)-ES. This contrasts with findings for the (µ/1, λ)-ES, the performance of which on sphere functions for n → ∞ monotonically deteriorates with increasing µ [START_REF] Beyer | The Theory of Evolution Strategies[END_REF].

• For large λ, the optimal truncation ratio is µ/λ = 0.27, and the corresponding progress per offspring is 0.202. Those values are identical to the optimal success probability and resulting normalized progress rate of the (1 + 1)-ES.

Beyer [START_REF] Beyer | The Theory of Evolution Strategies[END_REF] shows that the correspondence is no coincidence and indeed exact. The step-sizes that the two strategies employ differ widely, however. The optimal step-size of the (1+1)-ES is 1.224; that of the (µ/µ, λ)-ES is µc µ/µ,λ and for fixed truncation ratio µ/λ increases (slightly superlinearly) with the population size. For example, optimal step-sizes of (µ/µ, 4µ)-ES for µ ∈ {1, 2, 3} are 1.029, 2.276, and 3.538, respectively. If offspring candidate solutions can be evaluated in parallel, the (µ/µ, λ)-ES is preferable to the (1+1)-ES, which does not benefit from the availability of parallel computational resources.

Equation (15) holds in the limit n → ∞ for any finite value of λ. In finite but high dimensional search spaces, it can serve as an approximation to the normalized progress rate of the (µ/µ, λ)-ES on sphere functions in the vicinity of the optimal step-size provided that λ is not too large. A better approximation for finite n is derived in [START_REF] Beyer | Toward a theory of evolution strategies: On the benefits of sex -The (µ/µ, λ) theory[END_REF][START_REF] Beyer | The Theory of Evolution Strategies[END_REF] (however compare also [START_REF] Auger | Reconsidering the progress rate theory for evolution strategies in finite dimensions[END_REF]).

The improved performance of the (µ/µ, λ)-ES for µ > 1 compared to the strategy that uses µ = 1 is a consequence of the factor µ in the denominator of the term in Eq. ( 15) that contributes negatively to the normalized progress rate. The components z of mutation vectors selected for survival are correlated and likely to point in the direction opposite to the gradient vector. The perpendicular components z in the limit n → ∞ have no influence on whether a candidate solution is selected for survival and are thus uncorrelated. The recombinative averaging of mutation vectors results in a length of the z -component similar to those of individual mutation vectors. However, the squared length of the components perpendicular to the gradient direction is reduced by a factor of µ, resulting in the reduction of the negative term in Eq. ( 15) by a factor of µ. Beyer [START_REF] Beyer | Toward a theory of evolution strategies: On the benefits of sex -The (µ/µ, λ) theory[END_REF] has coined the term genetic repair for this phenomenon.

Weighted recombination (compare Algorithms 5 and 6) can significantly increase the progress rate of (µ/µ, λ)-ES on sphere functions. If n is large, the kth best candidate solution is optimally associated with a weight proportional to the expected value of the kth largest order statistic of a sample of λ independent standard normally distributed random numbers. The resulting optimal normalized progress rate per offspring candidate solution for large values of λ then approaches a value of 0.5, exceeding that of optimal unweighted recombination by a factor of almost two and a half [START_REF] Arnold | Weighted multirecombination evolution strategies[END_REF]. The weights are symmetric about zero. If only positive weights are employed and µ = λ/2 , the optimal normalized progress rate per offspring with increasing λ approaches a value of 0.25. The weights in Algorithms 5 and 6 closely resemble those positive weights.

(µ/µ, λ)-ES on Noisy Sphere Functions

Noise in the objective function is most commonly modeled as being Gaussian. If evaluation of a candidate solution x yields a noisy objective function value f (x) + σ N (0, 1), then inferior candidate solutions will sometimes be selected for survival and superior ones discarded. As a result, progress rates decrease with increasing noise strength σ . Introducing normalized noise strength σ * = σ n/(2R 2 ), in the limit n → ∞, the normalized progress rate of the (µ/µ, λ)-ES on noisy sphere functions is

ϕ * = σ * c µ/µ,λ √ 1 + ϑ 2 - σ * 2 2µ ( 18
)
where ϑ = σ * /σ * is the noise-to-signal ratio that the strategy operates under [START_REF] Arnold | Local performance of the (µ/µ I , λ)-ES in a noisy environment, Foundations of Genetic Algorithms[END_REF]. Noise does not impact the term that contributes negatively to the strategy's progress. However, it acts to reduce the magnitude of the positive term stemming from the contributions of mutation vectors parallel to the gradient direction. Notice that unless the noise scales such that σ * is independent of the location in search space (i.e., the standard deviation of the noise term increases in direct proportion to f (x), such as in a multiplicative noise model with constant noise strength), Eq. [START_REF] Yao | Evolutionary programming made faster[END_REF] describes progress in single time steps only rather than a rate of convergence. Figure 7 illustrates for different offspring population sizes λ how the optimal progress rate per offspring depends on the noise strength. The curves have been obtained from Eq. ( 18) for optimal values of σ * and µ. As the averaging of mutation vectors results in a vector of reduced length, increasing λ (and µ along with it) allows the strategy to operate using larger and larger step-sizes. Increasing the step-size reduces the noise-to-signal ratio ϑ that the strategy operates under and thereby reduces the impact of noise on selection for survival. Through genetic repair, the (µ/µ, λ)-ES thus implicitly implements the rescaling of mutation vectors proposed in [START_REF] Rechenberg | Evolutionsstrategie '94[END_REF] for the (1, λ)-ES in the presence of noise. (Compare c m and η x in Algorithms 5 and 6 that, for values smaller than one, implement the explicit rescaling). It needs to be emphasized though that in finite-dimensional search spaces, the ability to increase λ without violating the assumptions made in the derivation of Eq. ( 18) is severely limited. Nonetheless, the benefits resulting from genetic repair are significant, and the performance of the (µ/µ, λ)-ES is much more robust in the presence of noise than that of the (1+1)-ES.

Cumulative Step-Size Adaptation

All progress rate results discussed up to this point consider single time steps of the respective evolution strategies only. Analyses of the behavior of evolution strategies that include some form of step-size adaptation are considerably more difficult. Even for objective functions as simple as sphere functions, the state of the strategy is described by several variables with nonlinear, stochastic dynamics, and simplifying assumptions need to be made in order to arrive at quantitative results.

In the following we consider the (µ/µ, λ)-ES with cumulative step-size adaptation (Algorithm 4 with Eq. ( 1) in place of Lines 8 and 8b for mathematical convenience) and parameters set such that c σ → 0 as n → ∞ and d = Θ(1). The state of the strategy on noisy sphere functions with σ * = const (i.e., noise that decreases in strength as the optimal solution is approached) is described by the distance R of the parental centroid from the optimal solution, normalized step-size σ * , the length of the search path s parallel to the direction of the gradient vector of the objective function, and that path's overall squared length. After initialization effects have faded, the distribution of the latter three quantities is time invariant. Mean values of the time invariant distribution can be approximated by computing expected values of the variables after a single iteration of the strategy in the limit n → ∞ and imposing the condition that those be equal to the respective values before that iteration. Solving the resulting system of equations for σ * ≤ √ 2µc µ/µ,λ yields

σ * = µc µ/µ,λ 2 - σ * µc µ/µ,λ 2 (19)
for the average normalized mutation strength assumed by the strategy [START_REF] Arnold | Noisy Optimization with Evolution Strategies[END_REF][START_REF] Arnold | Performance analysis of evolutionary optimization with cumulative step length adaptation[END_REF]. The corresponding normalized progress rate

ϕ * = √ 2 -1 2 µc 2 µ/µ,λ 2 - σ * µc µ/µ,λ 2 (20)
is obtained from Eq. ( 18). Both the average mutation strength and the resulting progress rate are plotted against the noise strength in Fig. 8. For small noise strengths cumulative step-size adaptation generates mutation strengths that are larger than optimal. The evolution window continually shifts toward smaller values of the step-size, and adaptation remains behind its target. However, the resulting mutation strengths achieve progress rates within 20 percent of optimal ones. For large noise strengths the situation is reversed and the mutation strengths generated by cumulative step-size adaptation are smaller than optimal. However, increasing the population size parameters µ and λ allows shifting the operating regime of the strategy toward the left hand side of the graphs in Fig. 8, where step-sizes are near optimal. As above, it is important to keep in mind the limitations of the results derived in the limit n → ∞. In finite dimensional search spaces the ability to compensate for large amounts of noise by increasing the population size is more limited than Eqs. ( 19) and ( 20) suggest.

Parabolic Ridge Functions

A class of test functions that poses difficulties very different from those encountered in connection with sphere functions are ridge functions,

f (x) = x 1 + ξ n i=2 x 2 i α/2 = x 1 + ξR α ,
which include the parabolic ridge for α = 2. The x 1axis is referred to as the ridge axis, and R denotes the distance from that axis. Progress can be made by minimizing the distance from the ridge axis or by proceeding along it. The former requires decreasing step-sizes and is limited in its effect as R ≥ 0. The latter allows indefinite progress and requires that the step-size does not decrease to zero. Short and long term goals may thus be conflicting, and inappropriate step-size adaptation may lead to stagnation. As an optimal solution to the ridge problem does not exist, the progress rate ϕ of the (µ/µ, λ)-ES on ridge functions is defined as the expectation of the step made in the direction of the negative ridge axis. For constant step-size, the distance R of the parental centroid from the ridge axis assumes a time-invariant limit distribution. An approximation to the mean value of that distribution can be obtained by identifying that value of R for which the expected change is zero. Using this value yields

ϕ = 2µc 2 µ/µ,λ nξ(1 + 1 + (2µc µ/µ,λ /(nξσ)) 2 ) ( 21
)
for the progress rate of the (µ/µ, λ)-ES on parabolic ridge functions [START_REF] Oyman | Analysis of the (µ/µ, λ)-ES on the parabolic ridge[END_REF]. The strictly monotonic behavior of the progress rate, increasing from a value of zero for σ = 0 to ϕ = µc 2 µ/µ,λ /(nξ) for σ → ∞, is fundamentally different from that observed on sphere functions. However, the derivative of the progress rate with regard to the step-size for large values of σ tends to zero. The limited time horizon of any search as well as the intent of using ridge functions as local rather than global models of practically relevant objective functions both suggest that it may be unwise to increase the step-size without bounds.

The performance of cumulative step-size adaptation on parabolic ridge functions can be studied using the same approach as described above for sphere functions, yielding

σ = µc µ/µ,λ √ 2nξ (22)
for the (finite) average mutation strength [START_REF] Arnold | Evolution strategies with cumulative step length adaptation on the noisy parabolic ridge[END_REF]. From Eq. ( 21), the corresponding progress rate

ϕ = µc 2 µ/µ,λ 2nξ (23)
is greater than half of the progress rate attained with any finite step size.

Cigar Functions

While parabolic ridge functions provide an environment for evaluating whether step-size adaptation mechanisms are able to avoid stagnation, the ability to make continual meaningful positive progress with some constant nonzero step-size is of course atypical for practical optimization problems. A class of ridgelike functions that requires continual adaptation of the mutation strength and is thus a more realistic model of problems requiring ridge following are cigar functions

f (x) = x 2 1 + ξ n i=2 x 2 i = x 2 1 + ξR 2
with parameter ξ ≥ 1 being the condition number of the Hessian matrix. Small values of ξ result in sphere-like characteristics, large values in ridge-like ones. As above, R measures the distance from the x 1 -axis. Assuming successful adaptation of the step-size, evolution strategies exhibit linear convergence on cigar functions. The expected relative per iteration change in objective function value of the population centroid is referred to as the quality gain ∆ and determines the rate of convergence. In the limit n → ∞ it is described by

∆ * =          σ * 2 2µ(ξ -1) if σ * < 2µc µ/µ,λ ξ -1 ξ c µ/µ,λ σ * - σ * 2 2µ otherwise
where σ * = σn/R and ∆ * = ∆n/2 [START_REF] Arnold | On the behaviour of evolution strategies optimising cigar functions[END_REF]. That relationship is illustrated in Fig. 9 for several values of the conditioning parameter. The parabola for ξ = 1 reflects the simple quadratic relationship for sphere functions seen in Eq. ( 15). (For the case of sphere functions, normalized progress rate and normalized quality gain are the same.) For cigar functions with large values of ξ, two separate regimes can be identified. For small step-sizes, the quality gain of the strategy is limited by the size of the steps that can be made in the direction of the x 1 -axis. The x 1component of the population centroid virtually never changes sign. The search process resembles one of ridge following, and we refer to the regime as the ridge regime. In the other regime, the step-size is such that the quality gain of the strategy is effectively limited by the ability to approach the optimal solution in the subspace spanned by the x 2 , . . . , x naxes. The x 1 -component of the population centroid changes sign much more frequently than in the ridge regime, as is the case on sphere functions. We thus refer to the regime as the sphere regime.

The approach to the analysis of the behavior of cumulative step-size adaptation explained above for sphere and parabolic ridge functions can be applied to cigar functions as well, yielding

σ * = √ 2µc µ/µ,λ
for the average normalized mutation strength generated by cumulative step-size adaptation [START_REF] Arnold | On the behaviour of evolution strategies optimising cigar functions[END_REF]. The corresponding normalized quality gain is

∆ * =          ( √ 2 -1)µc 2 µ/µ,λ if ξ < √ 2 √ 2 -1 µc 2 µ/µ,λ ξ -1 otherwise .
Both are compared with optimal values in Fig. 10. For small condition numbers, (µ/µ, λ)-ES operate in the sphere regime and are within 20 percent of the optimal quality gain as seen above. For large condition numbers, the strategy operates in the ridge regime and achieves a quality gain within a factor of two of the optimal one, in accordance with the findings for the parabolic ridge above.

Further Work

Further research regarding the progress rate of evolution strategies in different test environments in- cludes work analyzing the behavior of mutative selfadaptation for linear [START_REF] Hansen | An analysis of mutative σ-selfadaptation on linear fitness functions[END_REF], spherical [START_REF] Beyer | Towards a theory of evolution strategies: Self-adaptation[END_REF], and ridge functions [START_REF] Meyer-Nieberg | Mutative selfadaptation on the sharp and parabolic ridge[END_REF]. Hierarchically organized evolution strategies have been studied when applied to both parabolic ridge and sphere functions [START_REF] Arnold | Hierarchically organised evolution strategies on the parabolic ridge[END_REF][START_REF] Beyer | On strategy parameter control by meta-ES[END_REF]. Several step-size adaptation techniques have been compared for ridge functions, including, but not limited to, parabolic ones [START_REF] Arnold | Step length adaptation on ridge functions[END_REF]. A further class of convex quadratic functions for which quality gain results have been derived is characterized by the occurrence of only two distinct eigenvalues of the Hessian, both of which occur with high multiplicity [START_REF] Arnold | On the use of evolution strategies for optimising certain positive definite quadratic forms[END_REF][START_REF] Beyer | Performance of the (µ/µ I , λ)-σSA-ES on a class of PDQFs[END_REF]. An analytical investigation of the behavior of the (1+1)-ES on noisy sphere functions finds that failure to reevaluate the parental candidate solution results in the systematic overvaluation of the parent and thus in potentially long periods of stagnation [START_REF] Arnold | Local performance of the (1 + 1)-ES in a noisy environment[END_REF]. Contrary to what might be expected, the increased difficulty of replacing parental candidate solutions can have a positive effect on progress rates as it tends to prevent the selection for survival of offspring candidate solutions solely due to favorable noise values. The convergence behavior of the (1+1)-ES on finite dimensional sphere functions is studied by Jebalia et al. [START_REF] Jebalia | Loglinear convergence and divergence of the scaleinvariant (1+1)-ES in noisy environments[END_REF] who show that the additive noise model is inappropriate in finite dimensions unless the parental candidate solution is reevaluated, and who suggest a multiplicative noise model instead. An analysis of the behavior of (µ, λ)-ES (without recombination) for noisy sphere functions finds that in contrast to the situation in the absence of noise, strategies with µ > 1 can outperform (1, λ)-ES if there is noise present [START_REF] Arnold | On the benefits of populations for noisy optimization[END_REF]. The use of non-singleton populations increases the signal-to-noise ratio and thus allows for more effective selection of good candidate solutions. The effects of non-Gaussian forms of noise on the performance of (µ/µ, λ)-ES applied to the optimization of sphere functions have also been investigated [START_REF] Arnold | A general noise model and its effects on evolution strategy performance[END_REF]. Finally, there are some results regarding the optimization of time-varying objectives [START_REF] Arnold | Optimum tracking with evolution strategies[END_REF] as well as analyses of simple constraint handling techniques [START_REF] Arnold | On the behaviour of the (1 + 1)-ES for a simple constrained problem, Parallel Problem Solving from Nature[END_REF][START_REF] Arnold | On the behaviour of the (1, λ)-ES for a simple constrained problem[END_REF][START_REF] Arnold | Analysis of a repair mechanism for the (1, λ)-ES applied to a simple constrained problem[END_REF].

Convergence Proofs

In the previous section we have described theoretical results that involve approximations in their derivation and consider the limit for n → ∞. In this section, exact results are discussed.

Convergence proofs with only mild assumptions on the objective function are easy to obtain for evolution strategies with a step-size that is effectively bounded from below and above (and, for non-elitist strategies, when additionally the search space is bounded) [START_REF] Auger | Theory of evolution strategies: A new perspective[END_REF][START_REF] Rudolph | Convergence Properties of Evolutionary Algorithms[END_REF]. In this case, the expected runtime to reach an -ball around the global optimum (see also Definition 3) cannot be faster than ∝ 1/ n , as obtained with pure random search for → 0 or n → ∞. 12 Similarly, convergence proofs can be obtained for adaptive strategies that include provisions for using a fixed step-size and covariance matrix with some constant probability.

Convergence proofs for strategy variants that do not explicitly ensure that long steps are sampled for a sufficiently long time typically require much stronger restrictions on the set of objective functions that they hold for. Such proofs however have the potential to reveal much faster, namely linear convergence. Evolution strategies with the artificial distance proportional step-size, σ = const× x , exhibit, as shown above, linear convergence on the sphere function with an associated runtime proportional to log(1/ ) [START_REF] Bienvenüe | Global convergence for evolution strategies in spherical problems: Some simple proofs and difficulties[END_REF][START_REF] Jebalia | Log-linear convergence and optimal bounds for the (1+1)-ES, Evolution Artificielle[END_REF][START_REF] Jebalia | Loglinear convergence and divergence of the scaleinvariant (1+1)-ES in noisy environments[END_REF][START_REF] Auger | Analyzing the impact of mirrored sampling and sequential selection in elitist evolution strategies[END_REF]. This result can be easily proved by using a law of large numbers, because x (t+1) / x (t) are independent and identically distributed for all t.

Without the artificial choice of step-size, σ/ x becomes a random variable. If this random variable is a homogeneous Markov chain and stable enough to satisfy the law of large numbers, linear convergence is maintained [START_REF] Bienvenüe | Global convergence for evolution strategies in spherical problems: Some simple proofs and difficulties[END_REF][START_REF] Auger | Theory of evolution strategies: A new perspective[END_REF]. The stability of the Markov chain associated to the self-adaptive (1, λ)-ES on the sphere function has been shown in dimension n = 1 [START_REF] Auger | Convergence results for (1,λ)-SA-ES using the theory of ϕ-irreducible Markov chains[END_REF] providing thus a proof of linear convergence of this algorithm. The extension of this proof to higher dimensions is straightforward.

Proofs that are formalized by upper bounds on the time to reduce the distance to the optimum by a given factor can also associate the linear dependency of the convergence rate in the dimension n. The (1 + λ)and the (1, λ)-ES with common variants of the 1/5th success rule converge linearly on the sphere function with a runtime of O(n log(1/ ) λ/ √ log λ) [START_REF] Jägersküpper | Algorithmic analysis of a basic evolutionary algorithm for continuous optimization[END_REF][START_REF]Jägersküpper: Probabilistic runtime analysis of (1 + ,λ) evolution strategies using isotropic mutations[END_REF]. When λ is smaller than O(n) the (1 + λ)-ES with a modified success rule is even √ log λ times faster and therefore matches the general lower runtime bound Ω(n log(1/ ) λ/ log(λ)) [START_REF]Jägersküpper: Probabilistic runtime analysis of (1 + ,λ) evolution strategies using isotropic mutations[END_REF]Theorem 5]. On convexquadratic functions, the asymptotic runtime of the (1+1)-ES is the same as on the sphere function and, at least in some cases, proportional to the condition number of the problem [START_REF]Jägersküpper: How the (1+1) ES using isotropic mutations minimizes positive definite quadratic forms[END_REF].

Convergence proofs of modern evolution strategies with recombination, of CSA-ES, CMA-ES or xNES are not yet available, however we believe that some of them are likely to be achieved in the coming decade.

= i z 2 i 1 .= µ k=1 |w k | 2 / µ k=1 w 2 k 1 ρ

 2121 portant theoretical results. Symbols and Abbreviations Throughout this chapter, vectors like z ∈ R n are column vectors, their transpose is denoted as z , and transformations like exp(z), z 2 , or |z| are applied component-wise. Further symbols are |z|= (|z 1 |, |z 2 |, . . . ) absolute value taken component wise z Euclidean length of a vector ∼ equality in distribution ∝ in the limit proportional to • binary operator giving the component-wise product of two vectors or matrices (Hadamard product), such that for a, b ∈ R n we have a • b ∈ R n and (a • b) i = a i b i . the indicator function, 1 α = 0 if α is false or 0 or empty, and 1 α = 1 otherwise. λ ∈ N number of offspring, offspring population size µ ∈ N number of parents, parental population size µ w , the variance effective selection mass or effective number of parents, where always µ w ≤ µ and µ w = µ if all recombination weights w k are equal in absolute value (1+1) elitist selection scheme with one parent and one offspring, see Section 2.1 (µ + , λ), e.g. (1+1) or (1, λ), selection schemes, see Section 2.1 (µ/ρ, λ) selection scheme with recombination (if ρ > 1), see Section 2.∈ N number of parents for recombination σ > 0 a step-size and/or standard deviation σ ∈ R n + a vector of step-sizes and/or standard deviations ϕ ∈ R a progress measure, see Definition 2 and Section 4.2

Figure 1 :

 1 Figure 1: Three 2-dimensional multivariate normal distributions N (0, C) ∼ C1 /2 N (0, I). The covariance matrix C of the distribution is, from left to right, the identity I (isotropic distribution), the diagonal matrix 1/4 0 0 4 (axis-parallel distribution) and ( 2.125 1.875 1.875 2.125 ) with the same eigenvalues (1/4, 4) as the diagonal matrix. Shown are in each subfigure the mean at 0 as small black dot (a different mean solely changes the axis annotations), two eigenvectors of C along the principal axes of the ellipsoids (thin black lines), two ellipsoids reflecting the set of points {x : (x -0) C -1 (x -0) ∈ {1, 4}} that represent the 1-σ and 2-σ lines of equal density, and 100 sampled points (however, a few of them are likely to be outside of the area shown).

Figure 2 :

 2 Figure2: Runs of the (1+1)-ES with constant step-size, of pure random search (uniform in [-0.2, 1] 10 ), and of the (1 + 1)-ES with 1/5th success rule (Algorithm 1) on a spherical function f (x) = x α , α > 0 (because of invariance to monotonic f -transformation the same graph is observed for any α > 0). For each algorithm there are three runs in the left plot and three runs in the right plot. The x-axis is linear in the left and in log-scale in the right hand plot. For the (1+1)-ES with constant step-size, σ equals 10 -2 . For the (1+1)-ES with 1/5th success rule, the initial step-size is chosen very small to 10 -9 and the parameter d equals 1 + 10/3. On the left, also the evolution of the step-size of one of the runs of the (1+1)-ES with 1/5th success rule is shown. All algorithms are initialized at 1. Eventually, the (1+1)-ES with 1/5th success rule reveals linear behavior on the left, while the other two algorithms reveal eventually linear behavior in the right hand plot.

  x k disappeared (compared to Algorithm 3) and the update of σ is postponed until after the for loop. Instead of the additional random variate ξ k , the length of the search path s σ determines the global stepsize change in Line 8b. For the individual step-size change, |z k | is replaced by |s σ |.

  5884720664009635e-08 blue:abs(f), cyan:f-min(f), green:sigma, red:axis ratio 0

Figure 3 :α 2 i y 2

 32 Figure 3: A single run of the (5/5 W , 10)-CMA-ES on the rotated ellipsoid function x → n i=1 α 2i y 2 i with α i = 10 3(i-1)/(n-1) , y = Rx, where R is a random matrix with R R = I, for n = 10. Shown is the evolution of various parameters against the number of function evaluations. Upper left: best (thick blue line), median and worst fitness value that reveal the final convergence phase after about 5500 function evaluations where the ellipsoid function has been reduced to the simple sphere; minimal and maximal coordinate-wise standard deviation of the mutation distribution and in between (mostly hidden) the step-size σ that is initialized far too small and increases quickly in the beginning, that increases afterwards several times again by up to one order of magnitude and decreases with maximal rate during the last 1000 f -evaluations; axis ratio of the mutation ellipsoid (square root of the condition number of C) that increases from 1 to 1000 where the latter corresponds to α n /α 1 . Lower left: sorted principal axis lengths of the mutation ellipsoid disregarding σ (square roots of the sorted eigenvalues of C, see also Fig.1) that adapt to the (local) structure of the underlying optimization problem; they finally reflect almost perfectly the factors α -1

Figure 5 :

 5 Figure 5: Decomposition of mutation vector z into a component z in the direction of the negative of the gradient vector of the objective function and a perpendicular component z .

Figure 6 :

 6 Figure 6: Maximal normalized progress per offspring of the (µ/µ, λ)-ES on sphere functions for n → ∞ plotted against the truncation ratio. The curves correspond to, from bottom to top, λ = 4, 10, 40, 100, ∞. The dotted line represents the maximal progress rate of the (1+1)-ES.

Figure 7 :

 7 Figure 7: Optimal normalized progress rate per offspring of the (µ/µ, λ)-ES on noisy sphere functions for n → ∞ plotted against the normalized noise strength. The solid lines depict results for, from bottom to top, λ = 4, 10, 40, 100, ∞ and optimally chosen µ. The dashed line represents the optimal progress rate of the (1+1)-ES [67].

Figure 8 :

 8 Figure 8: Normalized mutation strength and normalized progress rate of the (µ/µ, λ)-ES with cumulative step size adaptation on noisy sphere functions for n → ∞ plotted against the normalized noise strength. The dashed lines depict optimal values.

Figure 9 :

 9 Figure 9: Normalized quality gain of (µ/µ, on cigar functions for n → ∞ plotted against the normalized mutation strength for ξ ∈ {1, 4, 100}. The vertical line represents the average normalized mutation strength generated by cumulative step-size adaptation.

Figure 10 :

 10 Figure 10: Normalized mutation strength and normalized quality gain of the (µ/µ, λ)-ES with cumulative step-size adaptation on cigar functions for n → ∞ plotted against the condition number of the The dashed curves represent optimal values.

The sum of weights must be either one or zero, or recombination must be applied to the vectors x k -x and the result added to x.

[START_REF] Schwefel | Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie[END_REF] With discrete recombination, the offspring distribution is generated from a mixture of normal distributions with different mean values. The resulting distribution has lower entropy unless it has a larger overall variance.

Besides normally distributed mutations, Cauchy mutations[START_REF] Kappler | Are evolutionary algorithms improved by large mutations?, Parallel Problem Solving from Nature[END_REF][START_REF] Rudolph | Local convergence rates of simple evolutionary algorithms with Cauchy mutations[END_REF][START_REF] Yao | Evolutionary programming made faster[END_REF] have also been proposed in the context of evolution strategies and evolutionary programming.

[START_REF] Fogel | Artificial Intelligence through Simulated Evolution[END_REF] Using the normal distribution has several advantages. The N (0, I) distribution is the most convenient way to implement an isotropic perturbation. The normal distribution is stable: sums of independent normally distributed random variables are again normally distributed. This facilitates the design and analysis of algorithms remarkably. Furthermore, the normal distribution has maximum entropy under the given variances.

Optimality here means to achieve the largest expected approach of the optimum in a single generation.

Whenever sσ is large and therefore σ is increasing fast, the coefficient hσ prevents sc from getting large and quickly changing the distribution shape via C. Given hσ ≡ 1, under neutral selection sc ∼ N (0, C). The coefficient c h in line 11 corrects for the bias on sc introduced by events hσ = 0.

For a given C on the right hand side of Eq. (5), we have under neutral selection the stationarity condition E(Cnew) = C for the first three lines and E(log(Cnew)) = log(C) for the last line, where log is the inverse of the matrix exponential exp.

On a bounded domain and with mutation variances bounded away from zero, non-elitist strategies generate a subsequence of x-values converging to the essential global optimum.

If the mutation distribution is not normal and exhibits a singularity in zero, convergence can be much faster than with random search even when the step-size is bounded away from zero[START_REF] Zhigljavsky | Theory of Global Random Search[END_REF].