Thorsten Suttorp

Nikolaus Hansen
email: nikolaus.hansen@inria.fr

Christian Igel

Efficient Covariance Matrix Update for Variable Metric Evolution Strategies

Keywords: stochastic optimization, variable metric algorithm, multi-objective optimization, evolutionary algorithm, evolution strategy, covariance matrix adaptation, reinforcement learning

Randomized direct search algorithms for continuous domains, such as Evolution Strategies, are basic tools in machine learning. They are especially needed when the gradient of an objective function (e.g., loss, energy, or reward function) cannot be computed or estimated efficiently. Application areas include supervised and reinforcement learning as well as model selection. These randomized search strategies often rely on normally distributed additive variations of candidate solutions. In order to efficiently search in non-separable and ill-conditioned landscapes the covariance matrix of the normal distribution must be adapted, amounting to a variable metric method. Consequently, Covariance Matrix Adaptation (CMA) is considered state-of-the-art in Evolution Strategies. In order to sample the normal distribution, the adapted covariance matrix needs to be decomposed, requiring in general Θ(n 3 ) operations, where n is the search space dimension. We propose a new update mechanism which can replace a rank-one covariance matrix update and the computationally expensive decomposition of the covariance matrix. The newly developed update rule reduces the computational complexity of the rank-one covariance matrix adaptation to Θ(n 2 ) without resorting to outdated distributions. We derive new versions of the elitist Covariance Matrix Adaptation Evolution Strategy (CMA-ES) and the multi-objective CMA-ES. These algorithms are equivalent to the original procedures except that the update step for the

Introduction

Evolution strategies (ESs, [START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF][START_REF] Schwefel | Evolution and Optimum Seeking. Sixth-Generation Computer Technology Series[END_REF][START_REF] Beyer | Evolution strategies: A comprehensive introduction[END_REF] are randomized direct search algorithms. They are one of the major branches of evolutionary algorithms, a class of algorithms which draws inspiration from principles of neo-Darwinian evolution theory. Although ESs can be applied to various kinds of machine learning problems, the most elaborate variants are specialized for real-valued parameter spaces. Exemplary applications include supervised learning of feed-forward and recurrent neural networks, direct policy search in reinforcement learning, and model selection for kernel machines (e.g., [START_REF] Mandischer | A comparison of evolution strategies and backpropagation for neural network training[END_REF][START_REF] Igel | Optimization of dynamic neural fields[END_REF][START_REF] Schneider | Evolutionary adaptation of nonlinear dynamical systems in computational neuroscience[END_REF][START_REF] Igel | Neuroevolution for reinforcement learning using evolution strategies[END_REF][START_REF] Friedrichs | Evolutionary tuning of multiple SVM parameters[END_REF][START_REF] Kassahun | Efficient reinforcement learning through evolutionary acquisition of neural topologies[END_REF][START_REF] Pellecchia | Making driver modeling attractive[END_REF][START_REF] Mersch | Evolutionary optimization of sequence kernels for detection of bacterial gene starts[END_REF]Siebel and Sommer, 2007;Heidrich-Meisner and Igel, 2008a,b,c;Glasmachers and Igel, 2008, see below).

Evolution strategies and many other evolutionary algorithms for real-valued optimization, such as variants of evolutionary programming [START_REF] Fogel | Evolutionary Computation: Toward a New Philosophy of Machine Intelligence[END_REF], rely on Gaussian random mutations. Candidate solutions are altered by adding random vectors drawn according to multi-variate zero-mean normal distributions. Adaptation of the covariance matrices of these random variations during optimization means learning and employing an appropriate metric for the search process. It is well-known that an appropriate adaptation of the mutation strength (step size adaptation) is virtually indispensable and assists the approach of an optimum with increasing precision [START_REF] Schumer | Adaptive step size random search[END_REF][START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF][START_REF] Schwefel | Evolution and Optimum Seeking. Sixth-Generation Computer Technology Series[END_REF]. Additionally, the adaptation of the shape of the mutation distribution is appropriate. Adapting an arbitrary ellipsoidal shape [START_REF] Schwefel | Evolution and Optimum Seeking. Sixth-Generation Computer Technology Series[END_REF][START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF][START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] amounts to a variable-metric approach in which a Mahalanobis metric is learned. Adapting the distribution shape can improve the search performance by orders of magnitude, especially for non-separable or badly scaled objective functions [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Beyer | Evolution strategies: A comprehensive introduction[END_REF][START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF]. The covariance matrix adaptation ES (CMA-ES) implements both, adaptation of mutation strength and adaptation of distribution shape, explicitly, based on three major concepts: (a) the likelihood of successful steps is increased; (b) the correlation between successive steps in the past is exploited such that future successive steps tend to become uncorrelated (perpendicular); (c) invariance properties are sustained. The CMA-ES is regarded as one of the most powerful evolutionary algorithms for continuous optimization [START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF][START_REF] Beyer | Evolution strategies[END_REF].

In order to sample a multi-variate normal distribution, the covariance matrix has to be factorized. The covariance matrix update and the sampling of the normal dis-tribution require Θ(n 2 ) computations. Numerical routines for covariance matrix factorization require in general Θ(n 3 ) computations in n-dimensional search spaces. Consequently, on high dimensional objective functions that can be evaluated quickly, the factorization dominates the overall computational time. In order to reduce the computational burden on such objective functions to quadratic time, the factorization needs to be postponed until after Ω(n) iterations and consequently outdated distributions must be sampled.

In the following, we present a new update scheme for the covariances of Gaussian mutations in evolutionary algorithms that solves this problem. It operates directly on the factorization of the covariance matrix and reduces the computational complexity to Θ(n 2 ), the theoretical minimum for a matrix update. The update scheme also proves that it is always possible to update the covariance matrix and the distribution in each iteration with an asymptotic computational complexity of Θ(n 2 ). In addition, it simplifies the algebraic operations in the CMA-ES.

In the remainder of the introduction, we review performance evaluations and machine learning applications involving ESs with an emphasis on the CMA-ES. Then, we summarize previous work on reducing the time complexity of the covariance matrix update in ESs. In Section 2, we discuss Gaussian mutations and the adaptation of the mutation distribution in ESs. After that, we introduce our new update scheme. We incorporate the update scheme into the elitist (1+1)-CMA-ES and into the steadystate multi-objective CMA-ES (MO-CMA-ES, Igel et al, 2007a,b). Then, we introduce a new non-elitist strategy, denoted as (µ/µw, λ)-Cholesky-CMA-ES, which is a simplified version of the original non-elitist CMA-ES with rank-one updating. We empirically investigate the performance of the new algorithm in comparison with the original (µ/µw, λ)-CMA-ES with rank-one updating.

Performance Assessment and Machine Learning Applications of the CMA-ES

The CMA-ES is recognized as one of the most competitive evolutionary algorithms for real-valued optimization. Numerous performance evaluations and performance comparisons have been conducted across different suites of benchmark problems [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Hansen | Evaluating the CMA evolution strategy on multimodal test functions[END_REF][START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF]Auger and Hansen, 2005a,b;[START_REF] Whitley | Comparing the niches of CMA-ES, CHC and Pattern Search using diverse benchmarks[END_REF]. The number of objective function evaluations in the optimization process scales, with few exceptions, between linearly and quadratically with the search space dimension n [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF][START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF]. Linear scaling must be regarded as a general lower bound [START_REF] Jägersküpper | Lower Bounds for Hit-and-Run Direct Search[END_REF][START_REF] Jägersküpper | Lower bounds for randomized direct search with isotropic sampling[END_REF].

Arguably the most comprehensive performance comparison of CMA-ES with other bio-inspired, evolutionary and hybrid search methods for real parameter spaces took place in the Session on Real-Parameter Optimization at the 2005 IEEE Congress on Evolutionary Computation. 1 Performance results on 25 benchmark functions in 10, 30 and 50 dimensions were published for 11 algorithms (from 17 submissions), among them differential evolution, real-coded genetic algorithms, estimation of distribution algorithms, and hybrid particle swarm optimization. The (µ/µw, λ)-CMA-ES was applied in a restart setting, where the population size was increased by a factor of two before each restart (Auger and Hansen, 2005b). The restarted CMA-ES was the clear winner of the competition. It outperformed the best competitors by solving the most functions. In addition, it solved them for the most part considerably faster (often by an order of magnitude) whenever a quantitative assessment was available. On both the subset of all unimodal functions and on the subset of all multi-modal functions, the CMA-ES achieved the best overall results.

Summarizing the performance assessments with reference to benchmark function properties, we find that only on essentially separable2 objective functions CMA-ES can be significantly outperformed by other bio-inspired or hybrid algorithms. On essentially non-separable, ill-conditioned functions, CMA-ES generally outperforms them by orders of magnitude.

The CMA-ES has been successfully applied to a wide range of real-world problems (see the list by [START_REF] Hansen | [END_REF]. Most recently, [START_REF] Winter | Registration of CT and intraoperative 3D ultrasound images of the spine using evolutionary and gradientbased methods[END_REF] conducted a comparison of the CMA-ES with gradient-based approaches -including a conjugate gradient algorithm and the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method -in the domain of medical image registration, where the ES showed superior performance. In sum, [START_REF] Beyer | Evolution strategies[END_REF] observed that "CMA-ESs represent the state-of-the-art in evolutionary optimization in real-valued R n search spaces".

Using ESs for parameter optimization is particularly advisable in the presence of noise [START_REF] Arnold | Noisy Optimization With Evolution Strategies[END_REF][START_REF] Arnold | A comparison of evolution strategies with other direct search methods in the presence of noise[END_REF], where standard gradient-based approaches fail. [START_REF] Arnold | A comparison of evolution strategies with other direct search methods in the presence of noise[END_REF] compared ESs to several other direct search methods in a noisy environment and the results support the evolutionary approach. In machine learning, noisy objective functions are common. For example, they arise when the reward to be optimized is determined by a finite number of episodes (or rollouts), the standard scenario in reinforcement learning. They also arise in supervised learning and model selection when the learning system performance is assessed through randomly generated samples of data.

Evolution strategies and in particular the CMA-ES have been successfully applied in various fields of machine learning. The following examples highlight this success, emphasizing our own work.

Supervised Learning

The CMA-ES is well-suited for supervised learning in scenarios where the objective function is highly noisy and/or multi-modal, in particular if gradient information cannot be obtained efficiently. Although [START_REF] Mandischer | A comparison of evolution strategies and backpropagation for neural network training[END_REF] considered ESs for feedforward neural network training, the strength of this approach becomes more apparent when training recurrent systems. For example, [START_REF] Igel | Optimization of dynamic neural fields[END_REF] and [START_REF] Schneider | Evolutionary adaptation of nonlinear dynamical systems in computational neuroscience[END_REF] compared the CMA-ES with the BFGS based on analytically derived gradients for learning the parameters of neural fields, and the CMA-ES clearly outperformed the quasi-Newton method. [START_REF] Mierswa | Evolutionary learning with kernels: a generic solution for large margin problems[END_REF][START_REF] Mierswa | Controlling overfitting with multi-objective support vector machines[END_REF] investigated ESs for support vector machine (SVM) training, in particular in the case of non positive semi-definite kernel functions. In his recent work, he applied a multi-objective ES to SVM training, directly addressing the trade-off between model complexity and empirical risk [START_REF] Mierswa | Controlling overfitting with multi-objective support vector machines[END_REF]. This approach does not require an a priori chosen regularization parameter and can be viewed as a combination of training and model selection (see below).

Reinforcement Learning

Evolutionary algorithms are reinforcement learning (RL) methods in the definition of [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. They are typically used for direct policy search. [START_REF] Igel | Neuroevolution for reinforcement learning using evolution strategies[END_REF] brought forward the CMA-ES in the framework of RL. He found that the CMA-ES with rank-one updating outperforms alternative evolutionary RL approaches on variants of the pole balancing benchmark in fully and partially observable environments. In a more recent study by [START_REF] Gomez | Efficient non-linear control through neuroevolution[END_REF], these results were compared to 8-12 (depending on the task) other RL algorithms including value-function and policy gradient approaches. On the four test problems where the CMA-ES was considered, it ranked first, second (twice), and third. [START_REF] Pellecchia | Making driver modeling attractive[END_REF] applied the CMA-ES to model human car driving behavior using neural attractor dynamics to represent the policies. Siebel and Sommer (2007) and [START_REF] Kassahun | Efficient reinforcement learning through evolutionary acquisition of neural topologies[END_REF] employed the CMA-ES for RL in robotics. They combined the CMA-ES with evolutionary topology optimization to evolve artificial neutral networks.

Recently, Heidrich-Meisner and Igel (2008a,b,c) performed a systematic comparison between the CMA-ES and policy gradient methods with variable metrics. They discuss similarities and differences between these related approaches. Preliminary experiments -where the methods operated on the same class of policies and the same benchmark problems -indicate that the CMA-ES learns at least as quickly while being much more robust regarding the choice of hyperparameters and initial policies.

Model Selection

Real-valued evolutionary algorithms are often used for model selection of SVMs (see e.g., [START_REF] Runarsson | Asynchronous parallel evolutionary model selection for support vector machines[END_REF][START_REF] Friedrichs | Evolutionary tuning of multiple SVM parameters[END_REF][START_REF] Igel | Multi-objective model selection for support vector machines[END_REF][START_REF] Suttorp | Multi-objective optimization of support vector machines[END_REF][START_REF] Mersch | Evolutionary optimization of sequence kernels for detection of bacterial gene starts[END_REF][START_REF] Glasmachers | Uncertainty handling in model selection for support vector machines[END_REF]. Compared to gradient-based optimization, they do not require the class of kernel functions to have a differentiable structure, they are applicable even if the score function for assessing parameter performance is not differentiable, and they are less prone to converge to suboptimal solutions. Unlike nested grid-search, they also allow for adapting multiple hyperparameters. The CMA-ES was successfully applied to adapt general Gaussian kernels [START_REF] Friedrichs | Evolutionary tuning of multiple SVM parameters[END_REF][START_REF] Glasmachers | Uncertainty handling in model selection for support vector machines[END_REF] as well as string kernels [START_REF] Mersch | Evolutionary optimization of sequence kernels for detection of bacterial gene starts[END_REF]. Successful multi-objective model selection for SVMs using real-valued evolutionary algorithms is described by [START_REF] Igel | Multi-objective model selection for support vector machines[END_REF] and [START_REF] Suttorp | Multi-objective optimization of support vector machines[END_REF].

Previous Work in Complexity Reduction of the Covariance Matrix Update in ESs

Several heuristics have been proposed to reduce the time complexity of the covariance matrix update. Already in the first article on CMA-ES, [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF] suggested to postpone the update of the covariance matrix for n iterations. This, however, leads to sampling outdated distributions. Nevertheless, in practice, the update is only performed every τ generations. For τ = o(n) the computational complexity is still ω(n 2 ) if the algorithm relies on an eigenvalue decomposition, while τ = ω(n) is not advisable. Only when τ = Θ(n) is the computational complexity Θ(n 2 ), the same as for the update scheme developed in this article. This approach has three drawbacks. First, sampling outdated search distributions reduces the search performance in terms of progress achieved per objective function evaluation and therefore increases the overall running time. This is shown in Figure 1. Second, the decomposition leads to peak computations every τ generations. Third, an internal parameter, the update frequency, needs to be chosen properly. Fig. 1 Additionally needed function evaluations (in percent) to reach 10 -30 on f elli (see Table 2) when outdated distributions are sampled. The results refer to the default (µ/µ w, λ)-CMA-ES without rank-µ update of the covariance matrix [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF].

Curves for the standard variant with rank-µ update look virtually the same. The abscissa shows the number of steps the update is postponed divided by problem dimensionality. The lower curve depicts the median, the error bars the 10 th and 90 th percentile, for n = 100 dimensions and 11 trials for each case. The upper line shows the results for n = 200 dimensions from 3 trials. Waiting for n iteration steps leads to a loss between 12 and 14 percent. [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF] derived a distribution update in quadratic time under the constraint that only the variations from the previous iteration (i.e., from a single generation) are considered. This update rule was applied within simplified variants of the elitist CMA-ES for single-and multi-objective optimization [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF](Igel et al, , 2007b)). The simplification was necessary, because the covariance update in the CMA-ES is originally based on a weighted average of previously taken steps (the evolution path, see below) and not just on the last iteration. Depending on the optimization problem, this simplification can lead to a significant performance decrease both in the singleand multi-objective case as described by [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF]Igel et al ( , 2007b) ) and demonstrated in Section 4.1 below (Figure 4). [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] proposed to adapt the covariance only in an l-dimensional subspace, where l < n is chosen a priori. This restriction allows for an update step in O(l 2 n). [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] showed that the performance decreases significantly if the necessary dimensionality is underestimated. On a generic test function (f elli from below, where n = 30) the variant slows down by more than a factor of ten, even with l = n 2 (Knight and Lunacek, 2007, erratum). Algorithms that can learn covariances in linear time were proposed previously [START_REF] Ostermeier | An evolution strategy with momentum adaptation of the random number distribution[END_REF][START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF][START_REF] Poland | Main vector adaptation: A CMA variant with linear time and space complexity[END_REF]. Not surprisingly, these variants also suffer from a limited model complexity because the number of updated parameters is linear in n. They perform comparatively poorly whenever a full covariance matrix is required from the underlying search space landscape.

The algorithms presented in the following do not suffer from any of these limitations. They learn arbitrary covariance matrices, they do not sample outdated distributions, and they utilize an evolution path.

Limitations

Any search and optimization algorithm must exploit properties of the objective function, and the question arises on which classes of problems the CMA-ES will perform badly [START_REF] Igel | On classes of functions for which No Free Lunch results hold[END_REF]. First, the CMA-ES relies on neighborhood, assuming a correlation between the distance of solutions and the dissimilarity in their objective function values [START_REF] Jones | Fitness distance correlation as a measure of problem difficulty for genetic algorithms[END_REF]. The distance measure is based on the variable metric. When the neighborhood information is useless or misleading, CMA-ES will perform poorly (like most search algorithms). Second, invariance properties are a major design criterion for the CMA-ES (see below). Invariance, however, always coincides with lack of specialization. For example, a class of optimization problems might share a common coordinate system that is optimal for solving all instances. On such problems, algorithms that a priori utilize the distinguished coordinate system are potentially better than coordinate system invariant algorithms. Indeed, CMA-ES can be outperformed on separable problems, where the given coordinate system is highly distinguished. Similarly, it can be beneficial to sacrifice invariance under order preserving transformation of the function value (which is due to the rank-based selection in ESs). For example, on noise-free quadratic functions, quasi-Newton methods are typically an order of magnitude faster than CMA-ES, because they exploit the function values explicitly. Finally, the internal computations of CMA-ES scale at least with n 2 . On high dimensional quickly to evaluate objective functions, linearly scaling algorithms can be advantageous, even though they lack the ability to learn all mutual dependencies between variables (Ros and [START_REF] Hansen | [END_REF]. It is the objective of this work to reach the n 2 lower bound without any restrictions on the dependencies that can be learned and without resorting to outdated distributions.

Gaussian Mutations in Evolution Strategies

We consider ESs for real-valued optimization [START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF][START_REF] Schwefel | Evolution and Optimum Seeking. Sixth-Generation Computer Technology Series[END_REF][START_REF] Beyer | Evolution strategies: A comprehensive introduction[END_REF][START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF], which are quite well understood theoretically (cf. [START_REF] Beyer | The Theory of Evolution Strategies[END_REF][START_REF] Auger | Convergence results for the (1, λ)-SA-ES using the theory of ϕirreducible Markov chains[END_REF][START_REF] Jägersküpper | How the (1+1) ES using isotropic mutations minimizes positive definite quadratic forms[END_REF][START_REF] Jägersküpper | Lower bounds for randomized direct search with isotropic sampling[END_REF]. Let the optimization problem be defined by an objective function f : R n → Y to be minimized, where n denotes the dimensionality of the search space (space of candidate solutions, decision space) and Y the space of cost values. Evolution strategies are random search methods, which iteratively sample a set of candidate solutions from a probability distribution over the search space, evaluate these points using f , and construct a new probability distribution over the search space based on the gathered information. In typical ESs, this search distribution is parameterized by a set of candidate solutions, the parent population, and by parameters of the variation operators that are used to create new candidate solutions (the offspring) from the parent population.

In a canonical ES the objective vector x (g+1) i ∈ R n of the ith offspring at generation g is created by

x (g+1) i ← c (g) i recombination + v (g) i mutation .
The vector c is a realization of an ndimensional random vector distributed according to a zero-mean Gaussian distribution with covariance matrix

C (g) i , that is, v (g) i ∼ N 0, C (g) i .
Sampling the mutation distribution to generate v (g) i is usually conducted in two steps. First, the standard normal distribution is sampled to generate a realization of an ndimensional normally distributed random vector z (g) i ∼ N (0, I) with unit covariance matrix and zero mean. Second, this random vector is rotated and scaled by a linear transformation

A (g) i ∈ R n×n such that A (g) i z (g) i ∼ N (0, C (g) i ) for z (g) i ∼ N (0, I) .
To obtain A (g) i and to sample a general multivariate normal distribution, the covariance matrix

C (g) i is decomposed into Cholesky factors C (g) i = A (g) i A (g) i T .
Every symmetric nonnegative definite matrix, such as covariance matrices, has realvalued Cholesky factors. We do not require the matrices A (g) i to be of a special form [START_REF] Grewal | Kalman Filtering: Theory and Practice[END_REF], for example triangular. The factorization is in general not unique (by imposing certain constraints, e.g., A (g) i lower triangular with nonnegative diagonal elements, a unique decomposition can be defined). Numerical routines for computing Cholesky factors by a triangular Cholesky, eigenvalue, or singular value decomposition algorithm for a general covariance matrix require Θ(n 3 ) steps.

One of the decisive features of ESs is that the covariance matrices are subject to adaptation. The general strategy is to alter the covariance matrices such that steps promising a large fitness gain are sampled more often. There are typically two ways how the adaptation of the matrices is realized. First, the covariance matrix or its Cholesky factors can be parameterized, and these parameters can then be adapted. Either the parameterization or the adaptation rule has to ensure that the resulting matrices remain positive definite. In self-adaptive ESs the parameterization guarantees positive definiteness (e.g., see [START_REF] Rudolph | On correlated mutations in evolution strategies[END_REF], and the parameters are changed by mutation. Here we consider a second, more efficient way, where the covariance matrix is directly altered by additive updates of the form

C (g+1) = αC (g) + βV (g) .
The matrix V (g) ∈ R n×n is assumed to be positive definite and α, β ∈ R + are weighting factors. Let v (g) ∈ R be a step in the search space promising large fitness gain. To increase the probability that v (g) is sampled in the next iteration, the rank-one update

C (g+1) = αC (g) + βv (g) v (g) T (1)
is appropriate. This update rule shifts the mutation distribution towards the line distribution N 0, v (g) v (g) T , which is the Gaussian distribution with the highest probability to generate v (g) among all normal distributions with zero mean. After the update, the new covariance matrix has to be decomposed into Cholesky factors in order to sample the distribution. If the covariance matrix updates occur frequently in the ES the time needed for the factorization can dominate the computation time of the ES even for a moderate number of dimensions n.

Efficient Covariance Matrix Update

In general, each factorization of a covariance matrix requires ω(n 2 ) operations. Thus, in ESs with additive covariance matrix adaptation (e.g., using ( 1)) the Cholesky factorization of the covariance matrix is the computationally dominating factor apart from the fitness function evaluations. Therefore, we propose not to factorize the covariance matrix, but to use an incremental rank-one update rule for the Cholesky factorization. The idea is never to compute the covariance matrix explicitly, but to operate on Cholesky factors only. We consider Cholesky factors that are general n × n-matrices, and give a general update rule for ESs. The proposed technique makes use of efficient rank-one matrix updates (Hager, 1989), which are for example frequently used in the domain of Kalman filtering [START_REF] Grewal | Kalman Filtering: Theory and Practice[END_REF]. First, we derive an intermediate update rule for the special case where v (g) = A (g) z (g) and z (g) is given 3 [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF]. Under this assumption we can rewrite the rank-one update of the covariance matrix (1) as

C (g+1) = αC (g) + βA (g) z (g) A (g) z (g) T .
Our goal is to turn this update for C (g) into an update for A (g) . To achieve this, we derive the following lemma.

Lemma 1 [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF] Let u ∈ R n be an arbitrary vector. Then, the matrix I + uu T can be decomposed into

I + uu T = (I + ςuu T )(I + ςuu T ) with ς = 1 2 for u = 0 and ς = 1 u 2 1 + u 2 -1 otherwise.
3 The vector z (g) is, for example, the realization of a standard normally distributed random variable and an individual has been mutated by adding a vector proportional to v (g) = A (g) z (g) .

Proof

We have

(I + ςuu T )(I + ςuu T ) = I + 2ςuu T + ς 2 u 2 uu T = I + 2ς + ς 2 u 2 uu T = I + uu T .
The last step holds because 2ς + ς 2 u 2 = 1 by definition of ς. ⊓ ⊔ This result allows us to prove the following theorem, which shows how to realize an incremental update of the Cholesky factors given already appropriately factorized update vectors.

Theorem 1 [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF] Let C t ∈ R n×n be a symmetric positive definite matrix with Cholesky factorization C t = A t A T t . Assume further that C t is updated using

C t+1 = αC t + βv t v T t ,
where v t ∈ R n is given in the decomposed form v t = A t z t , and α, β ∈ R + . For z t = 0 a Cholesky factor of the matrix C t+1 can be computed by

A t+1 = √ αA t + √ α z t 2 1 + β α z t 2 -1 [A t z t ] z T t , (2)
for z t = 0 we have A t+1 = √ αA t .
Proof For z t = 0 it holds

C t+1 = αC t + βv t v T t = αA t A T t + βA t z t z T t A T t = αA t I + z t β α z T t A T t = √ αA t I + ςz t β α z T t I + ςz t β α z T t A T t √ α , where ς = α β 1 zt 2 1 + β α z t 2 -1 .
For the factorization in the last equality we used Lemma 1 with u = β α z t . It directly follows that

A t+1 = √ αA t I + ςz t β α z T t = √ αA t + β √ α ςA t z t z T t = √ αA t + √ α z t 2 1 + β α z t 2 -1 [A t z t ] z T t
computes the update of the Cholesky factor. The proof for the case z t = 0 is obvious.

⊓ ⊔

The square brackets in the last equation indicate the order of computation, showing how to achieve a time complexity of Θ(n 2 ). The numerical stability of direct updates of Cholesky factors is likely to be better than updates requiring subsequent decompositions (e.g., see the discussion in Grewal and Andrews, 1993, Chapter 6). Now we derive, as a new extension to the previously derived results, a corresponding update for the inverse of the Cholesky factors, which can then be used to bring arbitrary updates v t v T t in the form A t z t [A t z t ] T suitable for Theorem 1.

Theorem 2 Let A t , A t+1 ∈ R n×n , z t ∈ R n , and α, β ∈ R + be given such that Equation (2) of Theorem 1 holds. Let A -1 t be the inverse of A t . Then the inverse of A t+1 can be computed by

A -1 t+1 = 1 √ α A -1 t - 1 √ α z t 2   1 - 1 1 + β α z t 2   z t z T t A -1 t (3)
for z t = 0 and by

A -1 t+1 = 1 √ α A -1 t for z t = 0.
Proof Let z t = 0 and define a := √ α and b :=

√ α zt 2 1 + β α z t 2 -1 .
With these definitions the inverse of A t+1 can be calculated as follows:

A -1 t+1 = aA t + bA t z t z T t -1 = 1 a A -1 t - 1 a A -1 t bA t z t z T t 1 a A -1 t 1 + z T t 1 a A -1 t (bA t z t ) = 1 a A -1 t - b/a 2 1 + b/a • z t 2 z t z T t A -1 t = 1 a A -1 t - √ α 1 + β α z t 2 -1 α z t 2 • z t z T t A -1 t 1 + 1 + β α z t 2 -1 = 1 √ α A -1 t - 1 √ α z t 2   1 - 1 1 + β α z t 2   z t z T t A -1 t
The first equality follows from the definition of the update step for the Cholesky factor given by Theorem 1. For the second equality the Sherman-Morrison formula (B + 1989). Note, that if A 0 is invertible so is A t for all t ∈ N 0 . The proof for the case

xy T ) -1 = B -1 -B -1 xy T B -1 1+y T B -1 x is used with B = aA t , x = bA t z t and y = z t (Hager,
z t = 0 is obvious.
⊓ ⊔ Now we can combine the theorems to derive our main result:

Corollary 1 For covariance matrix updates of the form

C (g+1) = αC (g) + βv (g) v (g) T
the Cholesky factors of C (g+1) and their inverses can be computed in Θ(n 2 ) time given the decomposition C (g) = A (g) A (g) T and A (g) -1

. Let z (g) = A (g) -1 v (g) . For v (g) = 0 we have

A (g+1) -1 = 1 √ α A (g) -1 - 1 √ α z (g) 2   1 - 1 1 + β α z (g) 2   • z (g) z (g) T A (g) -1 (4)
and

A (g+1) = √ αA (g) + √ α z (g) 2 • 1 + β α z (g) 2 -1 v (g) z (g) T . (5)
For v (g) = 0, we have

A (g+1) -1 = 1 √ α A (g) -1 and A (g+1) = √ αA (g) .
Proof Theorem 1 can be applied for arbitrary vectors v (g) because the inverse of the Cholesky factor is known. In this case, we can compute g) with computational complexity Θ(n 2 ). The update of the Cholesky factor using this decomposition and Theorem 1 scales also with Θ(n 2 ). The required inverse of the Cholesky factor A (g) -1 can be recursively updated using Theorem 2 in Θ(n 2 ) steps.

z (g) = A (g) -1 v (g) and decom- pose v (g) into v (g) = A (g) A (g) -1 v (g) = A (g) z (
⊓ ⊔ Note that in applications of Corollary 1, many of the terms in the update equations are already available and need not to be re-computed. For an analogous rank-µ update 4 (e.g., see [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF][START_REF] Hansen | Evaluating the CMA evolution strategy on multimodal test functions[END_REF] the time complexity can be reduced accordingly to Θ(µn 2 ), where µ < n.

Applications

In this section, we show how the new covariance matrix update can be incorporated into state-of-the-art ESs. We consider recent variants of the covariance matrix adaptation ES (CMA-ES). The CMA-ES is our method of choice for real-valued evolutionary optimization for several reasons. First, the method follows appealing design principles (e.g., see [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], in particular it puts strong emphasis on invariance properties. Second, the method is quasi-parameter free (Auger and Hansen, 2005b). Third, as explicated above, the method performed well in benchmark scenarios (e.g., Auger and Hansen, 2005b) and, last but not least, also in many real-world applications [START_REF] Hansen | [END_REF]. 4 Rank-µ updates have the general form

C (g+1) = αC (g) + µ i=1 β i v (g) i v (g) i T with v (g) i ∈ R n and β i > 0 for 1 ≤ i ≤ µ. The n × n matrix V (g) = µ i=1 β i v (g) i v (g) i T
has rank µ if the vectors are linearly independent. The update can be realized iteratively,

C (g+1) = αC (g) + β 1 v (g) 1 v (g) 1 T + β 2 v (g) 2 v (g) 2 T
. . . . In practice, the rank-µ update is also used in the CMA-ES, where it is particularly useful with large offspring populations. Then, not only the weighted population mean, but also individual steps are considered in the update of the covariance matrix.

The key idea of all CMA-ES algorithms is to alter the mutation distribution such that the probability is increased, to reproduce steps in search space that led to the actual population (i.e., produced offspring that were selected). Using a covariance matrix update enables the algorithm to detect correlations between the objective variables and to become invariant under transformations of the search space. The algorithms implement several important concepts. The first one is known as derandomization: the adaptation of step size, variance, and covariance parameters depends deterministically on the realized selected steps in search space -in contrast to mutative self-adaptation [START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF][START_REF] Schwefel | Evolution and Optimum Seeking. Sixth-Generation Computer Technology Series[END_REF], where the dependence is stochastic. The second principle is cumulation: the trick of taking the search path of the population over the past generations into account (evolution path), where the influence of previous steps decays exponentially. Further, all CMA-ES variants decouple the update of a global step size and a covariance matrix, that is, the mutation distribution is parameterized by

v (g) i ∼ σ (g) N 0, C (g) i = σ (g) A (g) i N (0, I) . ( 6
)
The Cholesky factor A determines the shape of the distribution. Its update is based on the idea to reinforce successful directions. The update of step size σ, taking place on a considerably shorter time scale, is supposed to ensure nearly maximal progress for a given distribution shape.

The update of the mutation distribution is conducted such that the CMA-ES is invariant under transformations of the optimization problem, in particular invariant under rotation and translation of the search space (of course, apart from the initialization). If the mutation distribution is initialized accordingly, any affine transformation of the search space does not affect the performance of the CMA-ES. The single-objective CMA-ES is invariant under order-preserving transformations of the fitness function value. The multi-objective CMA-ES is invariant under order-preserving affine transformations f i → a i f i +b i with a i ∈ R + and b i ∈ R if the reference point for computing the hypervolume (x ref , see Section 4.2.1) is transformed accordingly or can be neglected, as it is the case in the standard MO-CMA-ES implementation.5

Evolution strategies use rank-based selection. Given µ parents and λ offspring, a (µ, λ) selection strategy chooses the µ best of the offspring to form the next parent population. A (µ+λ) strategy chooses the µ best individuals from the union of parents and offspring. The latter is an elitist selection method, because the best solution found so far is always a member of the population.

In the following, we present the elitist CMA-ES and the multi-objective MO-CMA-ES using the efficient learning of the metric. For completeness, we concisely explain the algorithms based on the original publications. Then, we introduce a new variant of the non-elitist CMA-ES compatible with the incremental covariance matrix update. For notational convenience, we do not consider the case of the mutation vector being zero (i.e., z = 0) in the description of the algorithms, a case which occurs with probability zero.

Elitist CMA-ES

The elitist CMA-ES is a combination of the well-known (1+1) selection scheme [START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF][START_REF] Schwefel | Evolution and Optimum Seeking. Sixth-Generation Computer Technology Series[END_REF][START_REF] Beyer | Evolution strategies: A comprehensive introduction[END_REF] with the covariance matrix adaptation as proposed for the non-elitist CMA-ES by Hansen andOstermeier (1996, 2001). In the elitist CMA-ES the adaptation of the covariance matrix employs the rank-one update from the original, non-elitist CMA-ES. However, the adaptation of the global step size is handled differently; the cumulative step size adaptation of the non-elitist CMA-ES is replaced by a success rule based step size control.

In this section we consider scalar fitness functions f : R n → R, x → f (x) to be minimized. We describe the original (1+1)-CMA-ES as proposed by [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF] and then incorporate the Cholesky update following Corollary 1. In the original algorithm, parent and offspring encode candidate solution vectors x parent , x offspring ∈ R n . We keep track of an averaged success rate p succ ∈ [0, 1], the global step size σ ∈ R + , an evolution path pc ∈ R n , and the covariance matrix C ∈ R n×n .

The algorithm is described in three routines. In the main part, termed (1+1)-CMA-ES and described in pseudo-code in the Algorithm 1 below, a new candidate solution is sampled and the parent solution x parent is replaced depending on whether the new solution x offspring is at least as good as x parent : Here the indicator function [•] is one if its argument is true and zero otherwise. Thus, [f(x offspring ) ≤ f (x parent )] is one if the last mutation has been successful and zero otherwise. The step size is updated depending on the success with a learning rate cp (0 < cp ≤ 1) using a target success rate p target succ :

Algorithm 1: (1+1)-CMA-ES initialize x parent , σ, C, p succ ← p target succ , pc ← 0 1 repeat 2 determine A such that C = AA T 3 z ∼ N (0, I) 4 x offspring ← x parent + σAz 5 updateStepSize σ, [f(x offspring ) ≤ f (x parent )], p succ 6 if f x offspring ≤ f (x parent ) then
Procedure updateStepSize(σ, λsucc, p succ )

p succ ← (1 -cp) p succ + cpλsucc 1 σ ← σ • exp 1 d p succ -p target succ 1 -p target succ 2
This procedure is rooted in the 1/5-success-rule proposed by [START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF]. The implementation follows a version proposed by [START_REF] Kern | Learning probability distributions in continuous evolutionary algorithms -A comparative review[END_REF], smoothed by means of line 1 and generalized in that it can also be applied for λ > 1. It implements the well-known heuristic that the step size should be increased if the success rate (i.e., the fraction of offspring not worse than the parent) is high, and the step size should be decreased if the success rate is low. The damping parameter d controls the extend of the step size changes.

If the new candidate solution is not worse than its parent, the covariance matrix is updated as in the original CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]:

Procedure updateCov(C, y, p succ , pc) if p succ < p thresh then 1 pc ← (1 -cc)pc + cc(2 -cc) y 2 C ← (1 -ccov)C + ccov • pcpc T 3 else 4 pc ← (1 -cc)pc 5 C ← (1 -ccov)C + ccov • pcpc T + cc(2 -cc)C 6
The update of the evolution path pc depends on the value of p succ . If the smoothed success rate p succ is high, that is, above p thresh < 0.5, the update of the evolution path pc is stalled. This prevents a too fast increase of axes of C when the step size is far too small, for example, in a (close to) linear surrounding. If the smoothed success rate p succ is low, the update of pc is accomplished obeying an exponential smoothing. The constants cc and ccov (0 ≤ ccov < cc ≤ 1) are learning rates for the evolution path and the covariance matrix, respectively. The factor cc(2cc) normalizes the variance of pc viewed as a random variable. The evolution path pc is used to update the covariance matrix. The new covariance matrix is a weighted mean of the old covariance matrix and the outer product pcpc T . In the second case (line 5), the second summand in the update of pc is missing and the length of pc shrinks. Although of minor relevance, the term cc(2cc)C (line 6) compensates for this shrinking in C.

Initial values are set to p succ = p target succ , pc = 0, and C = I, where p target succ is given in Table 1. The initial candidate solution x parent ∈ R n and the initial σ ∈ R + Table 1 Default parameters for the (1+1)-CMA Evolution Strategy.

Step size control:

d = 1 + n 2 , p target succ = 2 11
, cp = 1 12 Covariance matrix adaptation:

c c = 2 n + 2 , ccov = 2 n 2 + 6
, p thresh = 0.44 must be chosen problem dependent. The optimum should presumably be within the cube x parent ± 2 σ(1, . . . , 1) T . The (external) strategy parameters of the (1+1)-CMA-ES and their default values, taken from Igel et al (2007b) with λ set to one, are given in Table 1. Most default values are adopted from the precursor algorithms. In particular, the parameters for the covariance matrix adaptation are similar to those for the standard non-elitist CMA-ES. They have been validated by simulations on simple test functions for different dimensions, where the parameters were varied in a large interval. It has to be stressed that the parameters in Table 1 have to be rather viewed as constants and that in practice the user just needs to select the initial search point and the initial global step size (this is also true for the corresponding parameters of the ESs described below). We now discuss the setting for each parameter from Table 1 in turn.

d ≈ n/2 is the damping parameter for the step size controlling the amplitude of step size changes. In order to avoid large stochastic fluctuations the damping parameter should be as large as possible. On the other hand it must be small enough to permit fast convergence. According to the updateStepSize procedure, the step size changing multiplier can range, for the given target success rate, between exp(-1/(4.5 d)) ≈ 0.80 1/d and exp(1/d) ≈ 2.7 1/d . The step-size changing factor is monotonous in the success rate. Given p target succ = 2/11 and d = n/2, we get for success rates of 0, 1/10, 1/7, and 1/5.5, the respective changing factors of about 0.64, 0.82, 0.91, and 1 (constant step size) in n iteration steps. The factors justify the choice for d. A step size reduction of approximately 0.82 in n iterations corresponds to the maximal possible convergence rate that can be achieved with an (1+1)-ES and optimal step size [START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF]. p target succ ≈ 0.2 is the target success rate. A target success rate of 1/5 is well-established for the (1+1)-ES [START_REF] Beyer | Evolution strategies: A comprehensive introduction[END_REF] and was derived from optimal success rates on simple fitness functions [START_REF] Schumer | Adaptive step size random search[END_REF][START_REF] Rechenberg | Evolutionsstrategie: Optimierung Technischer Systeme nach Prinzipien der Biologischen Evolution. Frommann-Holzboog Ros R, Hansen N[END_REF]. We use a slightly smaller target success rate, because it generally leads to a larger adapted step size, which in turn is more likely to escape local optima, though we believe that this effect is small. cp ≈ 0.1 is the success rate averaging parameter, that is, the weight for the most recent value in the exponential smoothing of success events. The weight is chosen such that the success rate is smoothed over a significant period of about fifteen iterations, that is, a period with about two to three expected successes. While a large value for cp leads to larger stochastic fluctuations in σ, a small value impedes the response time of the adaptation process. cc ≈ 2/n is the weight for the most recent entry in the exponential smoothing for the evolution path that is used in the covariance matrix update (it appears in the decay term as prefactor 1cc). The backward time horizon of the cumulation process is approximately c -1 c . Only for c -1 c ∝ n the performance scale-up of the optimization on ridge like functions becomes linear with n [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF] which rules out considerably larger values, for example cc = 1 √ n . Too small values for cc on the other hand would require an undesirable reduction of the learning rate for the covariance matrix. ccov ≈ 2/(n 2 + 6) is the learning rate for the covariance matrix, where c -1 cov reflects the degrees of freedom of the covariance matrix with an additional adjustment for small n. p thresh < 0.5 is the threshold for stalling the update of evolution path and covariance matrix. Its value is chosen close to, but significantly smaller than one half: the update is stalled for large success probabilities. A success probability of one half occurs in a linear environment and will rarely be exceeded. Now we apply Corollary 1 and replace the covariance matrix update by the Cholesky update. The new main routine reads:

Algorithm 2: (1+1)-Cholesky-CMA-ES initialize x parent , σ, A, A inv , p succ ← p target succ , pc ← 0 1 repeat 2 z ∼ N (0, I) 3 x offspring ← x parent + σAz 4 updateStepSize σ, [f(x offspring ) ≤ f (x parent )], p succ 5 if f x offspring ≤ f (x parent ) then 6 x parent ← x offspring 7 updateCholesky (A, A inv , Az, p succ , pc) 8
until stopping criterion is met 9

In the (1+1)-Cholesky-CMA-ES the covariance matrix is never explicitly calculated and the update of the covariance is replaced by the corresponding update of the Cholesky factors and their inverse:

Procedure updateCholesky(A, A inv , z, p succ , pc) if p succ < p thresh then 1 pc ← (1 -cc)pc + cc(2 -cc) z 2 α ← (1 -ccov) 3 else 4 pc ← (1 -cc)pc 5 α ← (1 -ccov) + ccov • cc(2 -cc) 6 β ← ccov 7 w ← A inv • pc 8 A ← √ αA + √ α w 2 1 + β α w 2 -1 pcw T 9 A inv ← 1 √ α A inv - 1 √ α w 2 1 - 1 1+ β α w 2 w w T A inv 10
The resulting algorithm is equivalent to the original version, except that the computational complexity of a single generation is reduced from Θ(n 3 ) to Θ(n 2 ). The memory requirements of both variants of the elitist CMA-ES are the same. A further advantage of the (1+1)-CMA-ES with Cholesky update, termed (1+1)-Cholesky-CMA-ES in the remainder of this article, is its simple implementation. All considerations in this section are not restricted to the case of a single offspring, but also apply to the general (1+λ)-CMA-ES sampling λ offspring in each generation. The (1+λ)-CMA-ES is described by Igel et al (2007a), and it is straightforward to derive the corresponding (1+λ)-Cholesky-CMA-ES.

Experimental Evaluation

In this section, we present empirical results showing the validity and efficiency of our new update scheme. All experiments have been implemented using the Shark machine learning library [START_REF] Igel | Shark[END_REF].

Computational Complexity. We experimentally demonstrate the gain in computational efficiency achieved by using the new incremental update rule compared to the origi-Table 2 Single-objective benchmark functions used in our experiments. We define y = Ox, where O ∈ R n×n is an orthogonal matrix created randomly anew for each trial (each column vector is uniformly distributed on the unit hypersphere). In the noisy fitness function f sphereCauchy , ξ denotes a random variable drawn according to the standard Cauchy distribution.

Noisy Sphere, n = 20,

I init = [0.1, 0.3] n f sphereCauchy (x) = n i=1 x 2 i + 1 2n ξ n i=1 x 2 i Ellipsoid, n = 20, I init = [0.1, 0.3] n , c = 10 6 f elli (x) = n i=1 c i-1 n-1 y 2 i Cigar, n = 20, I init = [0.1, 0.3] n , c = 10 6 f cigar (x) = y 2 1 + n i=2 cy 2 i Generalized Rosenbrock's function, n = 20, I init = [0.1, 0.3] n f Rosenbrock (x) = n-1 i=1 100(y 2 i -y i+1 ) 2 + (y i -1) 2
nal (1+1)-CMA-ES, which decomposes the covariance matrix in every iteration. We measured the average time (on an Intel Xeona CPU with 2.40 GHz running Linux) needed for a covariance matrix update on the objective function f elli (see Table 2) depending on the problem dimension n. The average was computed over the first 100 updates.

The original elitist CMA-ES used a Cholesky factorization algorithm (Grewal and Andrews, 1993, section 6.4, p. 217), which was efficiently implemented.

The results are presented in Figure 2. The double log scale of the plot nicely shows the different exponents of the polynomial scaling of the algorithms. Already for 20 dimensions the new update is almost ten times faster. For n = 1000 the factor is about five hundred.

Stability. We performed several experiments on different objective functions to check the numerical stability of our update rules. Because the inverse of the Cholesky factor is updated independently of the Cholesky factor, one must suspect that they drift apart. Therefore, we monitored the evolution of AA inv -I F during the optimization. Here

• F denotes the Frobenius matrix norm. As a baseline comparison we additionally monitored AA -1 -I F , where the inverse A -1 of A is computed in every step using an accurate standard method (based on singular value decomposition, SVD). We considered the quadratic test problem f elli having a Hessian matrix with condition number of 10 6 for n = 3 and n = 200, and the Rosenbrock function f Rosenbrock that requires adaptation to a changing local metric for n = 20. As done in practice, we started with the unit covariance matrix, that is, perfect initial decomposition.

The results are shown in Figure 3. Both methods lose accuracy with time. Indeed, in the first 20 n generations the iterative updates are even more accurate than the inversion algorithm. After that, the iterative updates accumulate more errors, but even after 20 n 2 generations the error never exceeds 10 -11 .

Thus, we regard our update rules as highly accurate and numerically stable. Of course, it is possible to reset the iterative update every Ω(n) iterations computing the matrix inversion explicitly without affecting the run-time complexity. Our results show that this is not necessary at least before 50 n 2 iterations. Necessity of the Evolution Path. The new (1+1)-Cholesky-CMA-ES is a significant improvement over the method previously proposed by the authors in [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF]. Both methods implement an incremental, computational efficient update of the Cholesky matrix, but only the new one can utilize an evolution path. [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF] give several examples of the decline in optimization performance when the evolution path is omitted. Here, we replicated one of these experiments comparing the (1+1)-Cholesky-CMA-ES to the previously presented work.

Figure 4 shows the fitness evolution on the 20-dimensional f elli benchmark function based on 25 independent trials. The new method can learn the topology of the objective function about four times faster (within about 13 000 function evaluations). Also the second phase, the log-linear convergence to the optimum, is slightly faster. The shown result is, to our experience, representative for different dimensionalities and also for many objective functions.

The evolution path facilitates especially the learning of a single elongated axis in the shape of the distribution, the main task on the cigar benchmark function. Figure 5 shows the number of function evaluations to reach the target function value 10 -15 on f cigar depending on dimension n. With evolution path, the graph resembles virtually 300 n, the scaling is linear. Without evolution path the graph is close to 150 n 1.8 , the scaling is almost quadratic. On some objective functions only small differences can be observed and omitting the evolution path can even be (slightly) advantageous [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF]. Nevertheless in general, by exploiting the information on correlations between the steps in the evolution path, the new (1+1)-Cholesky-CMA-ES considerably outperforms the old one.

Multi-objective CMA-ES

Multi-objective optimization (MOO, or vector optimization) is concerned with the simultaneous optimization of multiple (scalar) objectives (e.g., [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF][START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]Jin, 2006). The goal of MOO is usually to find or to approximate the set of Pareto-optimal solutions. A solution is Pareto-optimal if it cannot be improved in one objective without getting worse in another one. A diverse set of Pareto-optimal solutions provides insights into the trade-offs between the objectives and serves as the basis for (human) decision-making. In recent years, evolutionary multi-objective algorithms have become popular for MOO (e.g., [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]Jin, 2006).

We consider real-valued MOO with m objectives. Each point x of the search space is assigned m objective function values f : R n → R m , x → (f 1 (x), . . . , fm(x)) T . We say that a solution x dominates a solution x ′ and write x ≺ x ′ , if and only if x is at least as good as x ′ in each objective and is strictly better than x ′ in at least one objective. The multi-objective CMA-ES (MO-CMA-ES) considers a population of µ solutions, where every solution reproduces and updates its search strategy as in the (1+1)-CMA-ES. The main additional problem to solve is the question of how to rank a set of solutions in order to determine the best µ candidate solutions forming the next parent population. The general approach to multi-objective ranking as taken in the MO-CMA-ES is introduced next, before we describe the search algorithm using the new update rule.

Ranking vector-valued solutions

The ranking in the MO-CMA-ES relies on non-dominated sorting as in the NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF][START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] and on the order induced by the contributing hypervolume as in the SMS-EMOA [START_REF] Beume | SMS-EMOA: Multiobjective selection based on dominated hypervolume[END_REF].

In a first step, the elements in a population X of candidate solutions are ranked according to their level of non-dominance. We denote the non-dominated solutions in

X by ndom(X) = {x ∈ X | ∄x ′ ∈ X : x ′ ≺ x}, where x ′ ≺ x means that x ′ dominates x. The Pareto front of X is then given by {(f 1 (x), . . . , fm(x)) | x ∈ ndom(X)}.
The elements in ndom(X) get rank 1. The other ranks are defined recursively by considering the set without the solutions with lower ranks (cf. [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]Igel et al, 2007a). Let dom 0 (X) = X, dom l (X) = dom l-1 (X) \ ndom l (X), and ndom l (X) = ndom(dom l-1 (X)) for l ∈ {1, . . . }. For x ∈ X the level of non-dominance r(x, X) is i iff x ∈ ndom i (X). Level of non-dominance is the first sorting criterion. We measured the number of iterations needed to reach an objective function value smaller than 10 -15 . The curves show the medians of 25 trials for n < 100 and of 5 trials for n > 100 on a double log scale. Error bars denote the 10 th and 90 th percentile.

A second sorting criterion is needed to rank the solutions that have the same level of non-dominance. This criterion is very important, as usually (in particular in the optimization of continuous objective functions) after some generations the number of non-dominated solutions in the population exceeds the number of solutions to be selected.

In the MO-CMA-ES, the contributing hypervolume serves as second sorting criterion. The hypervolume measure or S-metric was introduced by Zitzler and Thiele (1998) in the domain of evolutionary MOO. For a reference point x ref , the hypervolume Sx ref (X ′ ) of a set X ′ is defined as volume of the union of hypercuboids

Sx ref (X ′ ) = x∈ndom(X ′ ) {(f 1 (x ′ ), . . . , fm(x ′ )) | x ≺ x ′ ≺ x ref } The contributing hypervolume ∆ S (x, X ′ ) of a point x ∈ ndom(X ′ ) is given accord- ingly by ∆ S (x, X ′ ) := Sx ref (X ′ ) -Sx ref (X ′ \ {x}).
The rank s(x, X ′ ) of an individual a can be defined recursively based on its contribution to the hypervolume, where ties are broken at random. The individual contributing least to the hypervolume of X ′ gets the lowest rank. The individual contributing least to the hypervolume of X ′ without the individual with the lowest rank is assigned the second lowest rank and so on. We call x ∈ X ′ a boundary element if ∆ S (x, X ′ ) depends on the choice of the reference point x ref . The point x ref is chosen such that all elements in X ′ dominate x ref and that ∆ S (x, X ′ ) > ∆ S (x ′ , X ′ ) holds for any bound-ary element x ∈ X ′ and any non boundary element x ′ ∈ X ′ . That is, the individuals at the "boundaries" of the Pareto front of X ′ are preferably selected. Let a lower rank be worse. Formally (assuming that argmin breaks ties randomly), for

x ∈ ndom(X ′ ) it holds s(x, X ′ ) = 1 if x = argmin x ′ ∈X ′ {∆ S (x ′ , X ′ )} and s(x, X ′ ) = k if x = argmin x ′ ∈X ′ {∆ S (x ′ , X ′ \ {x ′′ | s(x ′′ , X ′ ) < k})}.
Based on this ranking and the level of non-dominance the relation

x ≺ X x ′ ⇔ r(x, X) < r(x ′ , X) or (r(x, X) = r(x ′ , X)) ∧ (s(x, ndom r(x,X) (X)) > s(x ′ , ndom r(x ′ ,X) (X))) ,
is defined on X. That is, x is better than x ′ when compared using ≺ X if and only if either x has a better level of non-dominance or x and x ′ are on the same level but x contributes more to the hypervolume when considering the points at that level of non-dominance.

Computing the hypervolume for many objectives is computationally demanding [START_REF] Klee | Can the measure of [a i , b i ] be computed in less than O(n • log n) steps?[END_REF][START_REF] Overmars | New upper bounds in Klee's measure problem[END_REF]. However, for 1 < m ≤ 3 special algorithms exist which solve this problem in O(µ m-2 ln µ), where µ is the number of points in the objective space [START_REF] Fonseca | An improved dimension-sweep algorithm for the hypervolume indicator[END_REF]. For m > 3, the current best algorithm scales with O(µ ln µ + µ m/2 ) [START_REF] Beume | Faster S-metric calculation by considering dominated hypervolume as Klee's measure problem[END_REF]. If computing the contributing hypervolume for selection turns out to be too time consuming in an application, the contributing hypervolume can be replaced by the change in quality measured by the ǫ-indicator [START_REF] Zitzler | Performance assesment of multiobjective optimizers: An analysis and review[END_REF], which can be computed efficiently.

The MO-CMA-ES

In the following, we describe the MO-CMA-ES working on Cholesky decompositions of the covariance matrices, for the original version we refer to Igel et al (2007a,b). The variant described here uses steady-state selection, that is, only one offspring is generated per generation. We consider steady-state selection where all µ members of the population are potential parents (i.e., (µ+1)-selection). This corresponds to the selection scheme used by [START_REF] Emmerich | An EMO algorithm using the hypervolume measure as selection criterion[END_REF] and [START_REF] Beume | SMS-EMOA: Multiobjective selection based on dominated hypervolume[END_REF]. Conceptually, this selection scheme has the advantage that the problem of selecting the µ out of µ + λ possible points such that the hypervolume is maximized can be solved easily for λ = 1.

An individual in the MO

-CMA-ES is a 6-tuple a = [x, p succ , σ, pc, A, A inv ] ,
where, x ∈ R n is the point in the search space, p succ ∈ R + 0 the average success rate, σ ∈ R + the global step size, pc ∈ R n the evolution path, and A and A inv are a Cholesky factor of the covariance matrix and its inverse. The steady-state MO-CMA algorithm can be described as follows.

Algorithm 3: steady-state MO-CMA initialize a k for k = 1, . . . , µ 1 Q ← a k 1 ≤ k ≤ µ 2 repeat 3 i ← U (1, µ) 4 a ← Q ≺:i 5 a ′ ← a 6 z ∼ N (0, I) 7 x ′ ← x + σAz 8 Q ← Q ∪ {a ′ } 9 if a ′ = Q ≺:(µ+1 ) then 10 updateStepSize σ ′ , a ′ ≺ Q a , p ′ succ 11 updateCholesky A ′ , A ′ inv , x ′ -x σ , p ′ succ , p ′ c 12 updateStepSize σ, a ′ ≺ Q a , p succ 13 Q ← Q \ Q ≺:(µ+1 )
14 until stopping criterion is met 15 At the beginning the population Q is initialized with µ parents. In the iteration loop, first a parent is chosen randomly in lines 4-5, where U (1, µ) denotes the discrete uniform distribution over {1, . . . , µ} and

Q ≺:i is the ith best individual in Q according to ≺ Q . One offspring a ′ = x ′ , p ′ succ , σ ′ , p ′ c , A ′ , A ′
inv is generated from the parent and added to the current population Q (lines 6-9). Then, step size and Cholesky factor of the offspring are updated, in case it is not the worst individual (lines 10-12). Finally, the step size of the parent is updated and the worst individual is removed from the population Q.

The described algorithm is equivalent to the MO-CMA-ES from Igel et al (2007a) while reducing its computational complexity from Θ(n 3 ) to Θ(n 2 ) per generation step. Igel et al (2007b) showed that using the CMA-ES covariance update rule without evolution path decreases the performance of the algorithm. However, in the new method an incremental covariance update is realized in combination with an evolution path. Hence, the search behaviors of the original generational and steady-state MO-CMA-ESs do not change when using the incremental update rule. All performance comparisons conducted for the original algorithms carry over to the new variants with computationally efficient covariance matrix update. This was validated empirically on the benchmark scenario described by Igel et al (2007b).

The (µ+1)-selection has proven to provide excellent performance on a great number of benchmark problems (Igel et al, 2007a,b) compared to NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF] and NSDE [START_REF] Iorio | Solving rotated multi-objective optimization problems using differential evolution[END_REF]. Recently, [START_REF] Voß | Scalarization versus indicator-based selection in multi-objective CMA evolution strategies[END_REF] compared the MO-CMA-ES with scalarization approaches for multi-objective optimization. For scalarization, the Tchebycheff method as well as the weighted-sum approach were considered [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF] and the (1+1)-CMA-ES served as single objective optimizer. The MO-CMA-ES clearly outperformed the scalarization approaches. For the detailed performance evaluation of the MO-CMA-ESs we refer to the comparisons by Igel et al (2007a,b) and [START_REF] Voß | Scalarization versus indicator-based selection in multi-objective CMA evolution strategies[END_REF].

Non-Elitist CMA-ES

In this section, we combine the incremental Cholesky update with a slightly simplified non-elitist (µ/µw, λ)-CMA-ES [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Beyer | Evolution strategies[END_REF][START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF], where µ/µw denotes weighted recombination of µ parental individuals. Non-elitism avoids systematic overvaluation6 in the presence of noise, and is thus of particular importance for many machine learning applications, as argued in the introduction.

The resulting new algorithm is a computationally more efficient variant of the original CMA-ES with rank-one updating [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF]. The need to simplify the original algorithm is due to its step size adaptation and is discussed in the following.

Cumulative Step Size Adaptation

In the (µ/µw, λ)-CMA-ES the global step size σ is adjusted using cumulative step size adaptation (CSA, sometimes denoted as path length control ). Let the offspring be generated using weighted global intermediate recombination of the selected parents, followed by Gaussian mutation with covariance matrix σ 2 C = σ 2 AA T . That is, prior to the mutation, the weighted center of mass x w = µ i=1 w i x i:λ is computed, where x i:λ is the ith best parent and w = (w 1 , . . . , wµ) T ∈ R µ are weighting coefficients with w 1 ≥ w 2 ≥ • • • ≥ wµ and w 1 = 1. Accordingly, we define z w = µ i=1 w i z i:λ as the weighted mean of the µ realizations of the standard normally distributed random vector that led to the µ best offspring. Thus, the step taken by the weighted mean of the population can be expressed as σA z w . The matrix σA scales and rotates the random vectors. We can remove all scaling effects by considering the orthogonal matrix B instead of σA, when the Cholesky factor A is given in the form A = BD with B ∈ R n×n orthogonal and D ∈ R n×n diagonal.

Standard CSA [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF][START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF] keeps track of a "conjugate" evolution path

pσ ← (1 -cσ)pσ + cσ(2 -cσ)µ eff B z w (7)
with learning rate cσ ∈ ]0, 1]. The evolution path pσ is a weighted sum of random vectors originally distributed according to N (0, I). Under random selection (when x and f (x) are independent) the normalization in (7) ensures that pσ, viewed as a random variable, is also distributed according to N (0, I) (this is explained in detail in the review by [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF]. The factor √ µ eff compensates for the loss of variance due to computing the weighted center of mass during recombination. The effective parent number, µ eff , also called variance effective selection mass, is given by µ eff =

1 µ i=1 w 2 i
. It is used for both step size control and covariance matrix adaptation.

In pσ, the influence of previous steps decays exponentially fast, controlled by the learning rate cσ. The conjugate evolution path can (and should) be learned on a faster timescale than the covariance matrix. [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF] propose the following learning rate

cσ = µ eff + 2 n + µ eff + 3 . ( 8
)
The update of the global step size σ is realized by

σ ← σ • exp cσ dσ pσ χn -1 . (9)
Here χn = E N (0, I) is the expected length of an n-dimensional, normally distributed random vector with covariance matrix I. The damping parameter dσ decouples the adaptation rate from the strength of the variation. Because of the proper normalization in ( 7), the expected length of pσ would be χn under random selection. Therefore, the global step size is increased if the steps leading to selected individuals are larger and/or more correlated than expected and decreased if they are smaller and/or more anticorrelated than expected in the absence of selection pressure.

Unfortunately, the incremental update rules of Corollary 1 do not provide the matrix B needed in (7). 7 Therefore, we omit B in the update of the global step size σ, similarly as in [START_REF] Beyer | Evolution strategies[END_REF],

pσ ← (1 -cσ)pσ + cσ(2 -cσ)µ eff z w . (10)
In this formulation, pσ viewed as a random variable would be still distributed according to N (0, I) with expected length χn under random selection. Hence, large or small selected steps still lead to increase or decrease of the step size, respectively. However, the identification of the degree of correlation of successive steps becomes an approximation, because B changes over time. Compared to the decay of the conjugate evolution path pσ in (10) the changes in B are small. To further diminish their effect, we increase the learning rate for the evolution path to

cσ = √ µ eff √ n + √ µ eff . (11)
In ( 11), the dependency of cσ on n is taken to its upper limit 1 √ n [START_REF] Hansen | Verallgemeinerte individuelle Schrittweitenregelung in der Evolutionsstrategie. Eine Untersuchung zur entstochastisierten, koordinatensystemunabhängigen Adaptation der Mutationsverteilung[END_REF]. Compared to (8), the interdependency between µ eff and n is preserved, in that for µ eff ≫ n the coefficient gets close to one.

When the Cholesky factor A remains constant, as it can be approximately the case in the final phase of optimization, the rules ( 9) and ( 10) result in the same update behavior for σ. Otherwise, we expect the approximation to be sufficiently accurate and only insignificantly affect the performance. This is confirmed by our experiments in Section 4.3.3.

7 In general, the incrementally updated Cholesky factor A is not symmetric. A symmetric Cholesky factor with A = BDB T would allow for an alternative formulation of (7) by replacing B z w with σ -1 A inv x (g+1) w x (g) w , thus not requiring B.

Non-Elitist CMA-ES with Incremental Cholesky Update

The simplified non-elitist (µ/µw, λ)-CMA-ES with Cholesky update, referred to as (µ/µw, λ)-Cholesky-CMA-ES in the following, is completely described in the following algorithm.

Algorithm 4: (µ/µw, λ)-Cholesky-CMA-ES initialize σ, x k for k = 1, . . . , λ, A ← I, A inv ← I, pσ ← 0, pc ← 0 1 repeat 2 x w ← µ i=1 w i x i:λ 3 for k = 1, . . . , λ do 4 z k ∼ N (0, I) 5 x k ← x w + σAz k 6 z w ← µ i=1 w i z i:λ 7 pc ← (1 -cc)pc + cσ(2 -cσ)µ eff A z w 8 v ← A inv • pc 9 A inv ← 1 √ 1-ccov A inv - 1 √ 1-ccov v 2 1 - 1 1+ ccov 1-ccov v 2 v v T A inv 10 A ← √ 1 -ccovA + √ 1-ccov v 2 1 + ccov 1-ccov v 2 -1 pcv T 11 pσ ← (1 -cσ)pσ + cσ(2 -cσ)µ eff z w 12 σ ← σ • exp cσ dσ pσ χn -1 13 until stopping criterion is met 14
In each iteration, λ offspring are generated by weighted global intermediate recombination (lines 4-7). Then, instead of the covariance matrix update in the original algorithm, the Cholesky factor and its inverse are updated according to Corollary 1 (lines 8-11). Finally, the global step size σ is adapted by the simplified CSA (lines 12 and 13).

The parameters for the (µ/µw, λ)-Cholesky-CMA-ES are given in Table 3. They are taken from [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF] for the case µcov = 1, using the increased learning rate for the conjugate evolution path (11) instead of (8). We describe the reasoning Table 3 Parameters for the (µ/µ w , λ)-Cholesky-CMA-ES. All parameters equal the parameters of the original variant, except for c σ .

Selection and recombination: λ = 4 + ⌊3 ln n⌋ , µ = ⌊λ/2⌋ , w i = ln(µ + 1) -ln i µ ln(µ + 1) -µ j=1 ln j for i = 1, . . . , µ

Step size control:

cσ = √ µ eff √ n + √ µ eff , dσ = 1 + 2 • max 0, µ eff -1 n + 1 -1 + cσ
Covariance matrix adaptation:

cc = 4 n + 4 , ccov = 2 (n + √ 2) 2
behind the choice of selection and recombination parameters, and of dσ in the following.

The remaining parameters were already discussed above (parameter cσ) or in Sect. 4.1 (parameters cc and ccov).

λ is the population size and its default value slowly increases with the dimension. The small setting for λ emphasizes fast convergence and reflects the point of view that an ES is a robust local search procedure. Smaller settings for λ can lead to a speed up, but too small settings endanger the working of the algorithm. With increasing population size the search becomes more and more global (given an adequate value for µ eff ). For considerably larger λ the rank-µ update for the covariance matrix should be used. µ and (w i ) i=1,...,µ denote parent number and recombination weights. The chosen recombination weights approximate the optimal setting on the sphere function [START_REF] Arnold | Weighted multirecombination evolution strategies[END_REF] for positive optimal weights. Negative weights are disregarded as, to our experience, they can be detrimental on non-spherical fitness functions. In general, in order to exploit the selection information effectively, the variance effective selection mass, µ eff , should lie between λ/5 and λ/2. With the given setting, µ eff := 1/ µ i=1 w 2 i equals approximately λ/3. dσ ≈ 1 controls the magnitude of step-size changes similar to d in Table 1. Equation ( 9) is formulated such that a good choice for dσ does not (heavily) depend on cσ. Smaller values for dσ, for example dσ = 0.5, can speed up the convergence. Too small values render the adaptation scheme unstable. Larger values will slow down the convergence without further harm. For large values of µ eff , dσ is increased and a rank-µ covariance matrix update should be used.

The original (µ/µw, λ)-CMA-ES requires an eigendecomposition of the covariance matrix. The eigendecomposition is more expensive and more difficult to implement than the Cholesky decomposition. As our main achievement, the (µ/µw, λ)-Cholesky-CMA-ES as described in this section does not need any of these procedures.

Experimental Evaluation

In order to evaluate the effect of the simplification of the CSA update, we compared the (µ/µw, λ)-Cholesky-CMA-ES to the (µ/µw, λ)-CMA-ES (using rank-one updating only) on 20-dimensional benchmark functions listed in Table 2.

The quadratic benchmark functions f elli and f cigar with rotated coordinate systems and condition numbers 10 6 are considered to investigate the learning of local metrics. They are particularly difficult, because they are non-separable and have high condition numbers. The ellipsoid f elli is our prototypical benchmark function having eigenvalues distributed equidistantly on the log-scale. The function f cigar tests the learning of a single elongated axis, see Section 4.1.1. The generalized Rosenbrock function f Rosenbrock with rotated coordinate system is difficult, because it is non-separable and the appropriate local metric changes in the course of optimization. Continuous adaptation of the covariance matrix is required. The sphere function with multiplicative Cauchy noise f sphereCauchy has been selected to demonstrate the performance in the case of noise. Replacing the Cauchy with a Gaussian noise distribution leads to similar results. As in all experiments, σ was initialized to 0.2/3 (one third of the initialization interval I init as given in Table 2). The (µ/µw, λ)-CMA-ES is run with the original value of cσ for updating the "conjugate" evolution path and with the increased learning rate, equations ( 8) and (11), respectively.

The results, based on 25 trials per test function and algorithm, are presented in Figures 6,7, and 8. The algorithms performed very similar as expected. On f elli , the (µ/µw, λ)-Cholesky-CMA-ES is slightly faster (about 5%, the difference is statistically significant) than the (µ/µw, λ)-CMA-ES with original learning rate cσ. The difference can be entirely attributed to the new choice of parameter cσ. The results on f cigar are qualitatively the same as on f elli and are therefore omitted. On the other two functions no relevant performance differences are observed. The shown results are representative and we have found, up to now, no case with a performance break down due to the simplification of CSA.

The standard (µ/µw, λ)-CMA-ES that includes the rank-µ update, reaches the function value of 10 -10 on f elli and f Rosenbrock around 20 % faster, and it performs on par on f sphereCauchy (results not shown).

Discussion

We presented a new, general rule for covariance matrix updates in variable metric realvalued direct search algorithms. The update rule operates directly on a factorization of the covariance matrix making the usually needed, repeated decompositions of the covariance matrix unnecessary. By simultaneously maintaining the inverse of the factors, the rule can replace arbitrary rank-one updates. This brings the following advantages:

-For n-dimensional problems, the new update requires Θ(n 2 ) computations and therefore reaches the asymptotically computational lower bound for the matrix update.

-We verified experimentally that the proposed matrix update is numerically stable and efficient in practice. For search space dimensions between n = 10 and 1000 the update proves to be about n/2 times faster than the otherwise necessary Cholesky decomposition.

-The proposed update is completely specified without hidden or numerically involved procedures such as an eigenvalue decomposition. Therefore, implementations in low level programming languages and even in hardware become easily possible.

The new learning rule can replace the original rank-one update of the covariance matrix and the subsequent matrix decomposition in the single-and multi-objective CMA-ES: -It can be used with the elitist (1+λ)-CMA-ES (Igel et al, 2007b) and the MO-CMA-ES (Igel et al, 2007b) without changing the algorithms. As opposed to previously introduced incremental update rules [START_REF] Igel | A computational efficient covariance matrix update and a (1+1)-CMA for evolution strategies[END_REF], restricted updates [START_REF] Ostermeier | An evolution strategy with momentum adaptation of the random number distribution[END_REF][START_REF] Hansen | On the adaptation of arbitrary normal mutation distributions in evolution strategies: The generating set adaptation[END_REF][START_REF] Poland | Main vector adaptation: A CMA variant with linear time and space complexity[END_REF][START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF]Ros and Hansen, 2008), or to using outdated distributions (Hansen andOstermeier, 1996, 2001), here an efficient Θ(n 2 ) update is achieved without loss in performance. -We combined the new update procedure with a slightly simplified non-elitist CMA-ES. The new algorithm performs in our experiments on par with the standard (µ/µw, λ)-CMA-ES with rank-one update [START_REF] Hansen | Completely derandomized self-adaptation in evolution strategies[END_REF], but lacks its geometrical interpretation of conjugate perpendicularity.

The matrix update can be straightforwardly extended to rank-µ updating of the covariance matrix [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF][START_REF] Hansen | Evaluating the CMA evolution strategy on multimodal test functions[END_REF]) leading to an algorithm with time complexity Ω(µn 2 ). However, whether the rank-µ update can be conducted more efficiently than outlined above, whether the (µ/µw, λ)-Cholesky-CMA-ES works well with larger values for µ, and detailed investigations of rank-µ algorithms remain subjects for future work. Two drawbacks of our new method can be recognized:

-In combination with the non-elitist CMA-ES the original algorithm needs to be slightly modified. This drawback can be rectified with a symmetrical factorization of the covariance matrix. Therefore, the question remains to be addressed in future whether an efficient matrix update can be found where the Cholesky factor is symmetrical. -The eigenvalues of the search distribution are not directly available for inspection.

However, the eigenspectrum gives insight into the structure of the underlying problem. But because rare additional computations of the eigenspectrum for inspection do not have any effect on the algorithm performance, this drawback is of minor relevance.

Concluding our discussion, we believe that the easier implementation and the improved performance for large scale problems will make the CMA-ES more popular in particular in the machine learning community.

Conclusions

Evolution Strategies (ESs) with variable metric are powerful tools in machine learning for search and optimization in continuous domains. They sample new candidate solutions according to a multi-variate normal distribution and adapt the covariance matrix to allow for efficient search even when the problem is non-separable and illconditioned. We proposed an incremental learning rule for the covariance matrix. We proved that it can equally replace the rank-one update together with the Cholesky decomposition of the covariance matrix, for example implemented in elitist single-and multi-objective covariance matrix adaptation ESs (CMA-ESs). The new learning rule reduces the computational complexity for a rank-one update of the search distribution to Θ(n 2 ), the asymptotically tight bound. The new method is considerably easier to implement, considerably faster in large dimensions and provides a significant improvement for high dimensional optimization and machine learning problems with fast computable performance measures.

  recombination scheme. For example, in case of weighted global intermediate recombination c (g) i is the weighted center of mass of the objective vectors in the current parent population. If no recombination is used, c

7x

  parent ← x offspring 8 updateCov (C, Az, p succ , pc) 9 until stopping criterion is met 10

Fig. 2

 2 Fig.2Comparison of the time needed for covariance updates in the elitist CMA-ESs depending on the problem dimension n. We measured the average time of the covariance matrix updates on f elli , based on the first 100 covariance updates. Shown is the median, the error bars indicate the 10 th and 90 th percentile. The standard elitist CMA-ES uses an efficient Cholesky decomposition for each update, while our new methods performs the updates of the Cholesky factor and its inverse incrementally. Note the double log scale.

Fig. 3

 3 Fig.3Accuracy of the incremental covariance matrix update compared to standard matrix inversion. Shown are typical trials on f elli with n = 3 and n = 200 as well as on f Rosenbrock with n = 20. The trials were stopped after an objective function value smaller than 10 -15 was reached. The divergence was checked every 10 th generation. Note the double log scale.

Fig. 4

 4 Fig.4Comparison of the elitist (1+1)-CMA-ES on the 20-dimensional f elli with and without evolution path. The curves show the medians of 25 trials. Both methods use an efficient incremental rank-one update, but the new algorithm learns the underlying metric about four times faster.

Fig. 5

 5 Fig.5Comparison of the elitist (1+1)-CMA-ES on f cigar with and without evolution path. We measured the number of iterations needed to reach an objective function value smaller than 10 -15 . The curves show the medians of 25 trials for n < 100 and of 5 trials for n > 100 on a double log scale. Error bars denote the 10 th and 90 th percentile.

  Fig. 6 Comparison of the (µ/µ w, λ)-CMA-ES with rank-one updating and the simplified (µ/µ w, λ)-Cholesky-CMA-ES on the benchmark function f elli for n = 20. The (µ/µw, λ)-CMA-ES with original and increased learning rate cσ are denoted by original and fast (µ/µw, λ)-CMA-ES, respectively. The trajectories show the medians of 25 trials, the final error bars the 10 th and 90 th percentile. The final difference between original (µ/µw, λ)-CMA-ES and (µ/µw, λ)-Cholesky-CMA-ES is significant with p < 10 -3 in the Wilcoxon rank-sum test. The curves of fast (µ/µw, λ)-CMA-ES and (µ/µw, λ)-Cholesky-CMA-ES can hardly be distinguished.

  Fig. 8 Comparison of the (µ/µw, λ)-CMA-ES with rank-one updating and the simplified (µ/µ w, λ)-Cholesky-CMA-ES on the noisy benchmark function f sphereCauchy for n = 20. The (µ/µ w , λ)-CMA-ES with original and increased learning rate c σ are denoted by original and fast (µ/µ w, λ)-CMA-ES, respectively. The trajectories show the medians of 25 trials, the final error bars the 10 th and 90 th percentile. The three curves can hardly be distinguished.

For details see http://www.ntu.edu.sg/home/EPNSugan/index_files/CEC-05/CEC05.htm.

Separable objective functions can be optimized coordinate-wise by n independent onedimensional optimizations and therefore are not subject to the curse of dimensionality.

We thank the anonymous reviewer for pointing out the extended invariance properties of the MO-CMA-ES.

Overvaluation[START_REF] Arnold | Noisy Optimization With Evolution Strategies[END_REF] refers to the effect that in the course of evolution individuals with their fitness sampled from the tail of the noise distribution are sustained. In order to compete with the previous parent, an individual needs to experience a similarly extreme noise event. This leads to continuously decreasing success rates and reduced progress.

Acknowledgments

We thank the editor and anonymous reviewers for their highly appreciated feedback.