Ilya Loshchilov
email: ilya.loshchilov@epfl.ch

A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization

Keywords: I.2.8 [Computing Methodologies]: Artificial Intelligence Problem Solving, Control Methods, and Search Evolution strategies, CMA-ES, large scale optimization, Cholesky update

, the LM-CMA-ES samples candidate solutions according to a covariance matrix reproduced from m direction vectors selected during the optimization process. The decomposition of the covariance matrix into Cholesky factors allows to reduce the time and memory complexity of the sampling to O(mn), where n is the number of decision variables. When n is large (e.g., n > 1000), even relatively small values of m (e.g., m = 20, 30) are sufficient to efficiently solve fully nonseparable problems and to reduce the overall run-time.

INTRODUCTION

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) is designed to learn dependencies between decision variables by adapting a covariance matrix which defines the sampling distribution of candidate solutions [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF]. This algorithm constantly demonstrates good performance at various platforms for comparing continuous optimizers such as Author's version .

the Black-Box Optimization Benchmarking (BBOB) workshop [START_REF] Auger | Comparison Tables of All Algorithms on All Noiseless Functions[END_REF][START_REF] Loshchilov | Bi-population CMA-ES agorithms with surrogate models and line searches[END_REF] and the Special Session at Congress on Evolutionary Computation [START_REF] García | A study on the use of non-parametric tests for analyzing the evolutionary algorithms' behaviour: a case study on the CEC'2005 Special Session on Real Parameter Optimization[END_REF][START_REF] Loshchilov | CMA-ES with restarts for solving CEC 2013 benchmark problems[END_REF]. The CMA-ES was also extended to noisy [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF], expensive [START_REF] Auger | Benchmarking the local metamodel CMA-ES on the noiseless BBOB'2013 test bed[END_REF][START_REF] Loshchilov | Self-adaptive surrogate-assisted covariance matrix adaptation evolution strategy[END_REF] and multi-objective optimization [START_REF] Igel | Covariance matrix adaptation for multi-objective optimization[END_REF].

The principle advantage of the CMA-ES, the learning of dependencies between n decision variables, also forms its main practical limitations such as O(n 2 ) memory storage and O(n 2 ) computational complexity per function evaluation [START_REF] Ros | A simple modification in CMA-ES achieving and space complexity[END_REF]. These limitations may preclude the use of the CMA-ES for computationally cheap but large scale optimization problems (e.g., with n > 100) if the internal computational cost of CMA-ES is greater than the cost of one function evaluation. On non-trivial large scale problems with n > 10000 not only the internal computational cost of the CMA-ES becomes significant but it is becoming simply impossible to efficiently store the covariance matrix in memory. One may argue that there are very few known continuous domain real-world problems of that huge dimensionality. This situation probably will not change much before practitioners have a set of tools that are able to efficiently search in such huge search spaces.

Several evolution strategies (ESs) have been proposed to deal with large scale optimization problems: O(n) time and space complexity algorithms such as separable CMA-ES (sep-CMA-ES [START_REF] Ros | A simple modification in CMA-ES achieving and space complexity[END_REF]) and linear time Natural Evolution Strategy (R1-NES [START_REF] Sun | A linear time natural evolution strategy for non-separable functions[END_REF]), L-CMA-ES [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] with O(m 2 n) time and O(mn) space complexity, where only m dominant eigen-pairs of the covariance matrix are computed. The sep-CMA-ES learns only the scaling of variables. The R1-NES learns only the predominant eigen-direction. The L-CMA-ES learns m dominant eigen-pairs, but its O(m 2 n) sampling complexity practically ends up with O(n 2 ) when m = √ n as studied in [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] for non-separable problems where multiple adaptation directions are required.

The problem of growing time and space complexity when optimizing large scale problems is not new. It was addressed in gradient-based optimization community when it became clear that for n > 1000 the storage of the approximate inverse Hessian matrix precludes the use of quasi-Newton methods such as Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [START_REF] Shanno | Conditioning of Quasi-Newton Methods for Function Minimization[END_REF]. As a solution, it was proposed not to store the matrix but to reconstruct it using information from the last m iterations [START_REF]Updating quasi-newton matrices with limited storage[END_REF]. The final algorithm called the limited memory BFGS algorithm (L-BFGS or LM-BFGS) is still considered to be the state-of-the-art of large scale gradientbased optimization [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF]. In this paper, we demonstrate that a very similar idea can be used to reconstruct the covari-ance matrix in the CMA-ES to reduce the time and space complexity to O(mn).

The paper is organized as follows. Section 2 reviews Evolution Strategies (ESs) proposed for large scale optimization. The LM-CMA-ES algorithm is described in section 3. The experimental validation of LM-CMA-ES is reported and discussed in section 4. Section 5 concludes the paper.

EVOLUTION STRATEGIES FOR LARGE SCALE OPTIMIZATION

Historically, first Evolution Strategies [START_REF] Rechenberg | Evolutionsstrategie: optimierung technischer systeme nach prinzipien der biologischen evolution[END_REF] were designed to perform the search without learning dependencies between variables which is a more recent development that gradually led to the CMA-ES algorithm [START_REF] Hansen | Adapting Arbitrary Normal Mutation Distributions in Evolution Strategies: The Covariance Matrix Adaptation[END_REF][START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF]. In this section, we discuss in detail the CMA-ES algorithm and its state-of-the-art derivatives for large scale optimization. For a recent comprehensible overview of Evolution Strategies, the interested reader is referred to [START_REF] Hansen | Evolution Strategies[END_REF].

The CMA-ES

The Covariance Matrix Adaptation Evolution Strategy [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 | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] is probably the most popular and in overall the most efficient Evolution Strategy.

The (µ/µw, λ)-CMA-ES is outlined in Algorithm 1. At iteration t of CMA-ES, a mean m t of the mutation distribution (can be interpreted as an estimation of the optimum) is used to generate its k-th out of λ candidate solution x k ∈ R n (line 5) by adding a random Gaussian mutation defined by a (positive definite) covariance matrix C t ∈ R n×n as

x t k = N m t , σ t 2 C t = m t + σ t N 0, C t , (1)
where σ t is a mutation step-size. These λ solutions then should be evaluated on an objective function f (line 6). The old mean of the mutation distribution is stored in m t and a new mean m t+1 is computed as a weighted sum of the best µ parent individuals selected among λ generated offspring individuals (line 7). The weights w are used to control the impact of selected individuals, weights are usually higher for better ranked individuals (line 1).

The procedure of the adaptation of the step-size σ t in CMA-ES is inherited from the Cumulative Step-Size Adaptation Evolution Strategy (CSA-ES) [START_REF] Hansen | Adapting Arbitrary Normal Mutation Distributions in Evolution Strategies: The Covariance Matrix Adaptation[END_REF] and is controlled by evolution path p t+1 σ . Successful mutation steps m t+1 -m t σ t (line 8) are tracked in the space of sampling, i.e., in the isotropic coordinate system defined by principal components of the covariance matrix C t . To update the evolution path p t+1 σ a decay/relaxation factor cσ is used to decrease the importance of previously performed steps with time. The step-size update rule increases the step-size if the length of the evolution path p t+1 σ is longer than the expected length of the evolution path under random selection E N (0, I ) , and decreases otherwise (line 13). Expectation of N (0, I ) is approximated by

√ n(1 -1 4n + 1 21n 2 )
. A damping parameter dσ controls the change of the step-size.

The covariance matrix update consists of two parts (line 12): rank-one update [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] and rank-µ update [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF]. The rankone update computes evolution path p t+1 c of successful moves of the mean m t+1 -m t σ t of the mutation distribution in the given coordinate system (line 10), in a similar way as the evolution path p t+1 σ of the step-size. To stall the update of p t+1 c when σ increases rapidly, a hσ trigger is used (line 9). 2 )-ln j) for i = 1 . . . µ, µw =

1 µ i=1 w 2 i , cσ = µw +2 n+µw +3 , dσ = 1 + cσ + 2 max(0, µw -1 n+1 -1), cc = 4 n+4 , c1 = 2 min(1,λ/6) (n+1.3) 2 +µw , cµ = 2 (µw -2+1/µw ) (n+2) 2 +µw 2: initialize m t=0 ∈ R n , σ t=0 > 0, p t=0 σ = 0, p t=0 c = 0, C t=0 = I, t ← 0 3: repeat 4:
for k = 1, . . . , λ do 5:

x k = m t + σ t N 0, C t 6:

f k = f (x k ) 7:
m t+1 ← µ i=1 wix i:λ // the symbol i : λ denotes i-th best individual on f 8: The rank-µ update computes a covariance matrix Cµ as a weighted sum of covariances of successful steps of µ best individuals (line 11). The update of C itself is a replace of previously accumulated information by a new one with corresponding weights of importance (line 12): c1 for covariance matrix p t+1 c p t+1 c T of rank-one update and cµ for Cµ of rank-µ update [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] such that c1 + cµ ≤ 1. Recently it was proposed to also take into account unsuccessful mutations in the "active" rank-µ update [START_REF] Hansen | Benchmarking a weighted negative covariance matrix update on the BBOB-2010 noiseless testbed[END_REF][START_REF] Jastrebski | Improving Evolution Strategies through Active Covariance Matrix Adaptation[END_REF].

p t+1 σ ← (1 -cσ)p t σ + cσ(2 -cσ) √ µwC t -1 2 m t+1 -m t σ t 9: hσ = 1 1 p t+1 σ < √ 1-(1-cσ ) 2(t+1) (1.4+ 2 n+1 ) E N(0,I ) 10: p t+1 c ← (1 -cc)p t c + hσ cc(2 -cc) √ µw m t+1 -m t σ t 11: Cµ = µ i=1 wi x i:λ -m t σ t × (x i:λ -m t ) T σ t 12: C t+1 = (1 -c1 -cµ)C t +
In CMA-ES, the factorization of the covariance C into AA T = C is needed to sample the multivariate normal distribution (line 5). The eigendecomposition with O(n 3 ) complexity is used for the factorization. Already in the original CMA-ES it was proposed to perform the eigendecomposition every n/10 generations (not shown in Algorithm 1) to reduce the complexity per function evaluation to O(n 2 )

Large Scale Variants

The original CMA-ES has O(n 2 ) time and space complexity that precludes its applications for large scale optimization with n ≫ 100. To enable the algorithm for large scale optimization, a linear time and space version called sep-CMA-ES was proposed in [START_REF] Ros | A simple modification in CMA-ES achieving and space complexity[END_REF]. The algorithm does not learn dependencies but the scaling of variables by restraining the covariance matrix update to the diagonal elements:

c t+1 jj = (1 -ccov)c t jj + 1 µcov p t+1 c 2 j + cccov 1 - 1 µccov µ i=1 wic t jj z i:λ t+1 2 j , j = 1, . . . , n (2)
where, for j = 1, . . . , n the cjj are the diagonal elements of C t and the z i:λ

t+1 j = x i:λ t+1 j /(σ t (cjj ).
This update reduces the computational complexity to O(n) and allows to exploit problem separability, thus the original property of being rotationally invariant is lost. The algorithm demonstrated good performance on separable problems and even outperformed CMA-ES on non-separable Rosenbrock function for n > 100.

A novel Natural Evolution Strategy (NES) variant, the Rank-One NES (R1-NES), which uses a low rank approximation of the search distribution covariance matrix was proposed recently by [START_REF] Sun | A linear time natural evolution strategy for non-separable functions[END_REF]. The algorithm adapts the search distribution according to the natural gradient with a particular parametrization of the covariance matrix,

C = σ 2 (I + uu T ), (3)
where u and σ are the parameters to be adjusted. The adaptation of the predominant eigen-direction u allows the algorithm to solve highly non-separable problems while maintaining only O(n) time and space complexity.

A version of CMA-ES with a limited memory storage also called limited memory CMA-ES (L-CMA-ES) was proposed by [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF]. The L-CMA-ES uses the m eigen-vectors and eigen-values spanning the m-dimensional dominant subspace of the n-dimensional covariance matrix C . The authors adapted a singular value decomposition updating algorithm developed in [START_REF] Brand | Fast low-rank modifications of the thin singular value decomposition[END_REF] that allowed to avoid the explicit computation and storage of the covariance matrix. For m < n the performance in terms of number of function evaluations gradually decreases while enabling the search in R n for n > 10000. However, the computational complexity of O(m 2 n) practically (for m in order of √ n [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF]) leads to the same limitations as for the original CMA-ES. The (µ/µw, λ)-Cholesky-CMA-ES proposed in [START_REF] Suttorp | Efficient covariance matrix update for variable metric evolution strategies[END_REF] is of special interest in this paper because the LM-CMA-ES is based on this algorithm. The Cholesky-CMA represents a version of CMA-ES with rank-one update where instead of performing the factorization of the covariance matrix C t into A t A t T = C t , the Cholesky factor A t and its inverse A t -1 are iteratively updated. From Theorem 1 [START_REF] Suttorp | Efficient covariance matrix update for variable metric evolution strategies[END_REF] it follows that if C t is updated as

C t+1 = αC t + βv t v t T , (4)
where v ∈ R n is given in the decomposition form v t = A t z t , and α, β ∈ R + , then for z = 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 , (5)
for zt = 0 we have A t+1 = √ αA t . From the Theorem 2 [START_REF] Suttorp | Efficient covariance matrix update for variable metric evolution strategies[END_REF] it follows that if A -1 t is 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 ], (6)
for z t = 0 and by

A -1 t+1 = 1 √ α A -1 t for z t = 0. The (µ/µw, λ)-Cholesky-CMA-ES is outlined in Algorithm 2. As well as in the original CMA-ES, Cholesky-CMA-ES Algorithm 2 The (µ/µw, λ)-Cholesky-CMA-ES 1: given n ∈ N+, λ = 4 + ⌊3ln n⌋, µ = ⌊λ/2⌋, wi = ln(µ+1)-ln(i) µln(µ+1)- µ j=1 ln(j) ; i = 1 . . . µ, µw = 1 µ i=1 w 2 i , cσ = √ µw √ n+ √ µw , dσ = 1 + cσ + 2 max(0, µw -1 n+1 -1), cc = 4 n+4 , c1 = 2 (n+ √ 2) 2 2: initialize m t=0 ∈ R n , σ t=0 > 0, p t=0 σ = 0, p t=0 c = 0, A t=0 = I, A t=0 inv = I, t ← 0 3: repeat 4:
for k = 1, . . . , λ do 5:

z k = N (0, I) 6:

x k = m t + σ t Az k 7:

f k = f (x k ) 8: m t+1 ← µ i=1 wix i:λ 9: zw ← µ i=1 wiz i:λ 10: p t+1 σ ← (1 -cσ)p t σ + cσ(2 -cσ) √ µwzw 11: p t+1 c ← (1 -cc)p t c + cc(2 -cc) √ µwAzw 12: v ← A t inv p c 13: A t+1 = √ 1 -c1A t + √ 1-c 1 v t 2 1 + c 1 1-c 1 v t 2 -1 p c v t T 14: A t+1 inv = 1 √ 1-c 1 A t inv - 1 √ 1-c 1 v t 2 1 - 1 1+ c 1 1-c 1 v t 2 v t [v t T A t inv ],
15:

σ t+1 ← σ t exp( cσ dσ ( p t+1 σ E N(0,I ) -1)) 16:
t ← t + 1 17: until stopping criterion is met proceeds by sampling λ candidate solutions (lines 4 -7) and taking into account the most successful µ out of λ solutions in the evolution paths adaptation (lines 10 and 11). However, the eigen-decomposition procedure is not required anymore because the Cholesky factor and its inverse are updated incrementally (line 13 and 14). This simplifies a lot the implementation of the algorithm and reduces its time complexity to O(n 2 ). A postponed update of the Cholesky factors every O(n) iterations would not reduce the asymptotic complexity further (as it does in the original CMA-ES) because the quadratic complexity will remain due to matrixvector multiplications needed to sample new individuals.

The non-elitist Cholesky-CMA is a good alternative to the original CMA-ES and demonstrates a comparable performance [START_REF] Suttorp | Efficient covariance matrix update for variable metric evolution strategies[END_REF]. While it has the same computational and memory complexity, the lack of rank-µ update may deteriorate its performance on problems where it is essential.

THE LM-CMA-ES

In this section, we first present main components of the computationally cheap limited memory CMA-ES and then introduce the algorithm itself. The components are: a procedure for reconstruction of Cholesky factors of a covariance matrix using stored direction vectors, a procedure to store these vectors and a new procedure for step-size adaptation.

Reconstruction of Cholesky factors

The idea to reconstruct the inverse Hessian matrix in the BFGS method [START_REF]Updating quasi-newton matrices with limited storage[END_REF] enabled its application for large scale gradient-based optimization. While the CMA-ES is a gradientfree algorithm, the two algorithms are indeed similar with a difference that the latter estimates the gradient in a stochas-Algorithm 3 Az(): Cholesky factor -vector update

1: given z ∈ R n , m ∈ Z+, j ∈ Z m + , P ∈ R m×n , V ∈ R m×n , b ∈ R m , a ∈ [0, 1] 2: initialize x ← z 3: for t = 1, . . . , min(m, |j|) do 4:
k ← b j t V (j t ,:) • x 5:

x ← ax + kP (j t ,:) 6: return x Algorithm 4 Ainvz(): inverse Cholesky factor -vector update 1:

given z ∈ R n , m ∈ Z+, j ∈ Z m , V ∈ R m×n , d ∈ R m , c ∈ [0, 1] 2: initialize x ← z 3: for t = 1, . . . , min(m, |j|) do 4: k ← d j t V (j t ,:) • x 5:
x ← cx -kV (j t ,:) 6: return x Algorithm 5 UpdateSet(): direction vectors selection

1: given m ∈ R + , j ∈ Z m + , l ∈ Z m + , t ∈ Z+, Nsteps ∈ Z+ 2: if t < m then 3: j t ← t 4: else 5: imin ← 1 + argmini l j i+1 -l j i , |1 ≤ i ≤ (m -1) 6:
if l j i minl j i min -1 ≥ Nsteps then 7:

imin ← 1 8:

if imin = m then 9: j tmp ← j i min 10:

for i = imin, . . . , m -1 do 11:

j i ← j i+1 12:
j m ← j tmp 13: jcur ← j min(t+1,m) 14: lj cur ← t 15: return: j cur , j, l tic way. This observation inspired us to investigate whether a similar matrix reconstruction procedure can be used in CMA-ES as well to reduce its time and space complexity.

As can be seen, the only use of Cholesky factor A t in Algorithm 2 is for sampling of new solutions after A t z k or for its own update to A t+1 . By setting a = √ 1 -c1 and

b t = √ 1-c 1 v t 2 1 + c 1 1-c 1 v t 2 -1 ,

one can rewrite the line (13) as

A t+1 = aA t + b t p t c v t T , (7)
In the following, we show how the vectors needed to sample new candidate solutions can be obtained without an explicit storage of Cholesky factors. At iteration t = 0, A 0 = I and A 0 z = z in line (6) of Algorithm 2, the new updated Cholesky factor

A 1 = aI + b 0 p 0 c v 0 T . At it- eration t = 1, A 1 z = (aI + b 0 p 0 c v 0 T )z = az + b 0 p 0 c (v 0 T z) and A 2 = a(aI + b 0 p 0 c v 0 T ) + b 1 p 1 c v 1 T
. Thus, a very simple iterative procedure which scales as O(mn) can be used to sample candidate solutions in IR n according to the Cholesky factor A t reconstructed from m pairs of vectors p t c and v t .

The pseudo-code of the procedure to reconstruct x = A t z from m direction vectors1 is given in Algorithm 3. At each iteration of reconstruction of x = A t z (lines 3 -4), x is updated as a sum of a-weighted version of itself and b tweighted evolution path p t c scaled by the dot product of v t and x. As can be seen, the Algorithm uses j(t) indexation instead of t. This is simply a convenient way to have references to matrices P and V which store p t c and v t vectors, respectively. In the next subsection, we will show how to efficiently manipulate these vectors.

A very similar approach can be used to reconstruct x = A t -1 z, for the sake of reproducibility the pseudo-code is given in Algorithm 4 for c = 1/ √ 1 -c1 and

d t = 1 √ 1-c 1 v t 2 × × 1 - 1 1+ c 1 1-c 1 v t 2
. The computational complexity of both procedures scales as O(mn).

Direction Vectors Selection and Storage

It is an open question how to use only m ≪ n direction vectors to obtain a comparable amount of useful information as stored in the covariance matrix of the original CMA-ES. For large n and λ ≪ n, evolution path vectors p t c from the last m iterations are likely to be quite similar and therefore to contain only some local information.

In this paper, we propose a simple approach which forces m selected vectors to be at approximately the same distance from each other in terms of number of iterations, but at most with the distance of Nsteps from each other given that the m-th vector is the one from the last iteration. This selection procedure is outlined in Algorithm 5 which outputs an array of pointers j such that j 1 points out to a row in matrices P and V with the oldest saved vectors p c and v which will be taken into account during the reconstruction procedure. The higher the index i of j i the more recent the corresponding direction vector is. The index jcur points out to the oldest vector which will be replaced by the newest one in the same iteration when the procedure is called. The rule to choose a vector to be replaced is the following: find a pair of consecutively saved vectors with the closest distance (in terms of number of iterations, stored in l) between each other (line 5), if this distance is smaller than Nsteps then the most recent vector will be removed by assigning jcur ← imin, otherwise the oldest vector among m saved vectors should be removed. Thus, the procedure gradually replaces vectors in a way to keep them at approximately the same distance, but at most at distance of Nsteps iterations.

Population Success Rule

An elegant success rule for step-size adaptation called the median success rule was recently proposed in [START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF]. It is applicable to non-elitist multi-recombinant evolution strategies. The median success rule compares the median fitness of the population to a fitness from the previous iteration. The comparison fitness is chosen to achieve a target success rate of 1/2. The empirical validation demonstrated that the median success rule is competitive to CSA [START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF].

In practice, one should count the number Ksucc of individuals in the current population better than some j-th best individual of the previous population, where j depends on n and λ but can be set to be 0.3λ [START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF]. Then, a normalized measurement

z ← 2 λ Ksucc - λ + 1 2 (8)
can be computed such that z ≥ 0 iff the median individual was successful.

The step-size is adapted as

σ ← σ exp s dσ , (9)
where s ← (1cσ)s + cσz and ds = 2(n -1)/n. We suppose that while being quite elegant the median success rule has a potential drawback that we will demonstrate on an example. Let us suppose that fitness values (to be minimized) of the previous population are say f t-1 = [2.1, 3.1, 4.1, 5.1, 6.1, 7.1, 8.1] while the fitness values of the current population are f t = [START_REF] Auger | Benchmarking the local metamodel CMA-ES on the noiseless BBOB'2013 test bed[END_REF][START_REF] Auger | Comparison Tables of All Algorithms on All Noiseless Functions[END_REF][START_REF] Brand | Fast low-rank modifications of the thin singular value decomposition[END_REF][START_REF] Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF][START_REF] García | A study on the use of non-parametric tests for analyzing the evolutionary algorithms' behaviour: a case study on the CEC'2005 Special Session on Real Parameter Optimization[END_REF][START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF][START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF]. According to the median success rule if j is chosen as, e.g., 3, the number of successful individual (with fitness values better than or equal to f t-1 (3) = 4.1) is 4 (as f t (1) = 1, f t (2) = 2, f t (3) = 3 and f t (4) = 4). The computed value of Ksucc is then will be used to adapt the step-size. However, its computation does not take into account the values of f t-1 (i) for 1 ≤ i < j and even if all such f t-1 (i) are better than the best solution f t (1), this information will not be taken into account.

This potential drawback is not the drawback in a sense that the median success rule was designed in this way. However, we suppose that the information omitted in the median success rule can be useful since it can provide a better estimate whether and by how much the new population is more successful than the previous one.

In this paper, we introduce the population success rule (PSR) for step-size adaptation for non-elitist multi-recombinant evolution strategies. To estimate the success of the current population we combine fitness function values from the previous and current population into a mixed set

f mix ← f t-1 ∪ f t (10)
Then, we rank all individual in the mixed set to define two sets rt-1 and rt containing ranks of individuals of the previous and current populations ranked in the mixed set.

We compute a normalized success measurement

zP SR ← λ i=1 rt(i) -rt-1(i) λ 2 -z * , (11)
where z * is a target success ratio. The step-size can be adapted as in [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]. The proposed population success rule takes into account all fitness function values from the previous and current generation. This success rule seems to represent a more general case of the 1/5th-rule which can be obtained when λ = 1.

The Algorithm

In the previous subsection we introduced all necessary components of the (µ/µw, λ)-LM-CMA-ES outlined in Algorithm 6. The algorithm represents a computationally efficient limited memory version of CMA-ES, where the Cholesky factor and its inverse are reconstructed from a set of stored direction vectors (lines 6 and 9). The mutation step-size Algorithm 6 The (µ/µw, λ)-LM-CMA-ES 1:

given n ∈ N+, λ = 4 + ⌊3ln n⌋, µ = ⌊λ/2⌋, wi = ln(µ+1)-ln(i) µln(µ+1)- µ j=1 ln(j) ; i = 1 . . . µ, µw = 1 µ i=1 w 2 i , cσ = 0.3, dσ = 1, m = 4 + ⌊3ln n⌋, Nsteps = m, cc = 1 m , c1 = 1 10ln(n+1) 2: initialize m t=0 ∈ R n , σ t=0 > 0, p t=0 c = 0, s ← 0, t ← 0 3: repeat 4:
for k = 1, . . . , λ do 5:

z k = N (0, I) 6:

x k = m t + σ t Az(z k ) 7:

f t k = f (x k ) 8: m t+1 ← µ i=1 wix i:λ 9: p t+1 c ← (1 -cc)p t c + cc(2 -cc) √ µw(m t+1 -m t )/σ 10: v ← Ainvz(p t+1 c ) 11:
jcur ← U pdateSet() 12:

V (jcur ,:) ← v; P (jcur ,:) ← p t+1 c 13:

bj cur ← √ 1-c 1 v t 2 1 + c 1 1-c 1 v t 2 -1 14: dj cur = 1 √ 1-c 1 v t 2 1 - 1 1+ c 1 1-c 1 v t 2 , 15: r t , r t-1 ← Ranks of f t and f t-1 in f t ∪ f t-1 16: zP SR ← λ i=1 r t (i)-r t-1 (i) λ 2 -z * 17: s ← (1 -cσ)s + cσzP SR 18:
σ t+1 ← σ t exp(s/dσ) 19:

t ← t + 1 20: until stopping criterion is met is adapted using the population success rule (lines 15 -13). The algorithm memory and time complexity scales as O(mn).

SIMULATION RESULTS

In this section, we perform a set of numerical experiments to assess the performance of the proposed LM-CMA-ES on large scale optimization problem with n = 128, 256, . . . , 8096. We investigate the performance on three basic problems: Sphere function f Sphere (x) = n i=1 x 2 i , separable Ellipsoid function f Elli (x) = n i=1 10 6 i-1 n-1 x 2 i and its rotated version f ElliRot (x) = f Elli (Qx), where Q is an orthogonal n × n matrix with each column vector q i being a uniformly distributed unit vector implementing an angle-preserving transformation [START_REF] Ros | A simple modification in CMA-ES achieving and space complexity[END_REF].

Experimental Setting

For the sake of reproducibility, the MATLAB/C++ source code of all tested algorithms is available at https://sites.google.com/site/lmcmaeses/.

In the order to estimate the performance of (µ/µw, λ)-LM-CMA-ES, we compare it with (µ/µw, λ)-Cholesky-CMA-ES and (µ/µw, λ)-Sep-CMA-ES. We use the default parameters for Cholesky-CMA-ES and Sep-CMA-ES as given in [START_REF] Suttorp | Efficient covariance matrix update for variable metric evolution strategies[END_REF] and [START_REF] Ros | A simple modification in CMA-ES achieving and space complexity[END_REF], respectively. The parameters of LM-CMA-ES are given in Algorithm 6. For all problems, the mean m t=0 is initialized in the range [-5, 5] n , the population is sampled with initial step-size σ t=0 = 5 and using the same seed per run. Note that in all cases we use the default population size λ = 4 + ⌊3ln n⌋.

Memory and Computational Complexity of LM-CMA-ES

The LM-CMA-ES has O(mn) memory complexity and more specifically stores Q ∈ R m×n , V ∈ R m×n and λ solution vectors xi. For large n and m = λ def ault = λ = 4 + ⌊3ln n⌋ used in this paper, the algorithm stores approximately 3mn real-valued parameters. If a real-valued parameter requires 8 bytes of memory, then for n = 8192 the LM-CMA-ES will require 5.8 megabytes while the original CMA-ES would start to reach its limit by requiring 1 gigabyte of memory. Using the same amount of memory (more specifically, 1.03 gigabyte), the LM-CMA-ES will able to optimize a 1 million dimensional problem. Indeed, by taking m = 1 even less memory would be needed but the latter possibility makes sense only if the performance stays at a reasonable level.

Figure 1 shows how fast CPU time per evaluation scales for different operations (measured on a 2.0 GHz processor). Scalar-vector multiplication of a vector with n variables scales linearly with ca. 6 • 10 -10 n seconds, evaluation of the separable Ellipsoid is twice more expensive if a temporary data is used. Sampling of n normally distributed variables scales as ca. 100 vectors-scalar multiplications. As can be seen in Figure 1 The computation cost of CMA-ES with full covariance matrix learning limits its applicability for n ≫ 100 and makes it intractable because of memory for n > 10000.

Performance on Sphere and Ellipsoid

The Sphere function is often viewed in Evolutionary Computation to be the first function to look at when benchmarking evolutionary algorithms. Figure 2 demonstrates a comparable performance of LM-CMA-ES with population success rule, sep-CMA-ES with CSA and Cholesky-CMA-ES with CSA. It should be further studied what is the effect of the target population success rate (set to z * = 0.25) whose value was chosen the same for all experiments in order to obtain a reasonable performance on Ellipsoid functions.

Figure 3-Left shows that both LM-CMA-ES and Cholesky-CMA-ES are rotationally invariant and therefore they optimization runs (one per function) are almost coincide (within the algorithm). The sep-CMA-ES is not rotationally invariant and therefore it performs better on the separable Ellipsoid than on its rotated version where the exploitation of the separability is not that useful. Importantly, the LM-CMA-ES often outperforms the Cholesky-CMA-ES in the beginning of optimization, while the adaptation of the full covariance matrix makes Cholesky-CMA-ES faster at later stages. Figure 3-Right shows that the loss of performance of LM-CMA-ES compared to Cholesky-CMA-ES is in order of a factor of 3-4 given that for n = 2048 the LM-CMA-ES uses only m = 26 direction vectors. It is important to keep in mind that for n > 10000 the Cholesky-CMA-ES becomes intractable both due to its memory and computational complexity. Then, the sep-CMA-ES becomes an alternative, however, it does not learn dependencies and might be therefore inefficient (see Figure 3-Left). We discussed several large scale ESs in this paper: L-CMA-ES [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] and R1-NES [START_REF] Sun | A linear time natural evolution strategy for non-separable functions[END_REF]. We compared the LM-CMA-ES indirectly by analyzing the results from [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] and [START_REF] Sun | A linear time natural evolution strategy for non-separable functions[END_REF]. It takes about 6000 seconds for L-CMA-ES to solve 200-dimensional Ellipsoid after about 7e + 6 function evaluations with m = √ n = 14 and 4000 seconds after 4e + 6

evaluations with m = n/2 = 100. The LM-CMA-ES solves the same problem after about 125 seconds and 5.3e + 6 function evaluations with m = 19. The performance is comparable while the LM-CMA-ES is about 32 -48 times faster that is unlikely to be only due to a different processor or implementation used. The L-CMA-ES has O(m 2 n) computational complexity and therefore it is in order of m times computationally slower than LM-CMA-ES. The R1-NES algorithm performs well on non-separable problems but tends to fail on problems where the learning of multiple principal components is essential, e.g., it fails on moderate dimensional rotated Ellipsoid function [START_REF] Sun | A linear time natural evolution strategy for non-separable functions[END_REF]. On Rosenbrock function the LM-CMA-ES is about 5 times faster (not shown) in terms of number of function evaluations for n = 256, 512. The R1-NES also samples from the normal distribution, and therefore the lower bound of its computational complexity is predefined (see Figure 1).

We performed an experiment on 100,000-dimensional separable Ellipsoid problems for 100,000 function evaluations (i.e., n evaluations). The original CMA-ES and Cholesky-CMA-ES cannot be applied due to memory requirements. The applicability of L-CMA-ES is also limited due to its O(m 2 n) computational complexity. The results for sep-CMA-ES specifically designed for large scale optimization and the proposed LM-CMA-ES are shown in Figure 4. While the LM-CMA-ES gradually improves the fitness similarly as in Figure 3-Left, the sep-CMA-ES does not improve it because it diverges from the very first iterations. To investigate whether it is a mistake in our implementation, we launched the same experiment using the sep-CMA-ES author's MATLAB implementation where the divergence was also observed.

It should be noted that the separable Ellipsoid can be easily solved by various Evolutionary Algorithms which implicitly or explicitly exploit its separability, our purpose of its usage is to investigate how the LM-CMA-ES performs on problems with high dependencies between variables. Given that the LM-CMA-ES is rotationally invariant, its performance on both separable and non-separable problems is comparable, but the former is cheaper to compute.

DISCUSSION AND CONCLUSION

This paper presents a new approach to efficiently store and exploit the information about dependencies between decision variables of large scale optimization problems. It allows to reconstruct the Cholesky factor and its inverse using m ≪ n direction vectors that turns out to be sufficient to obtain good performance on large scale problems with highlydepended variables. The implementation of this approach in the LM-CMA-ES algorithm makes it possible to optimize a 1 million dimensional problem while learning dependencies between variables at a cost of about 0.1 second per function evaluation on an ordinary machine. Indeed, one should not plan to easily find a global optimum in such a huge search space, but some local optimization/tuning seems reasonable, e.g., in Machine Learning problems.

The proposed LM-CMA-ES algorithm is based on the population success rule which looks promising and requires further theoretical and empirical investigations. It should be studied as well whether it can be claimed to represent a general case of the 1/5th success rule. More experiments are required to investigate whether and when the lack of rank-µ update is a limitation. All parameters chosen for the algorithm were tuned only moderately and specifically for large n and might require a significant revision to address a wider set of optimization problems commonly used for EAs. However, we suppose that the performance on the Ellipsoid function is already worth a closer scientific investigation. We envision that several directions may further improve the algorithm: i) adaptation of m within a fixed range, the impact of m itself should be studied as well, ii) since the population success rule does not make any assumptions about the sampling distribution, the Gaussian sampling can be removed that would further speed-up the algorithm (e.g., to replace CSA by PSR in CMA-ES).

The speculations about a possibility of having CMA-ES like evolutionary processes going on in nature often end up around a hypothesis that there is no such a thing in natural evolution as a full covariance matrix and its update. One may suppose that only a limited number of direction vectors is stored to adjust the mutation in promising directions.

Algorithm 1

 1 The (µ/µw, λ)-CMA-ES 1: given n ∈ N+, λ = 4 + ⌊3ln n⌋, µ = ⌊λ/2⌋, wi = ln(µ+1 2 )-ln i µ j=1 (ln(µ+ 1

  evaluation in sec.

Figure 1 :

 1 Figure 1: Timing results of LM-CMA-ES on the separable Ellipsoid compared to sep-CMA-ES and Cholesky-CMA-ES. The results were computed using at most 10 5 function evaluations for sep-CMA-ES and LM-CMA-ES and using at most 10 4 for Cholesky-CMA-ES.

  , sampling of z k dominates the computational overhead of sep-CMA already after n = 128. The LM-CMA-ES scales almost linearly for n ≥ 1024 as ca. 1.3 • 10 -7 n or ca. 200 scalar-vector multiplications. Matrix-vector multiplication scale quadratically with n and Cholesky-CMA-ES scales as ca. 1.5-2 matrix-vector multiplications. Practically, the LM-CMA-ES is about 40 times faster (in terms of its internal computation cost per function evaluation) for n = 2048 and about 140 times faster for n = 8192 than Cholesky-CMA-ES. The LM-CMA-ES is only about 2

Figure 2 :

 2 Figure 2: Results of LM-CMA-ES on the Sphere function compared to sep-CMA-ES and Cholesky-CMA-ES. Lines show the median of 11 runs for different problem dimensions to reach the target fitness value of 10 -10 . The dotted line is an extrapolation.

Figure 3 :

 3 Figure 3: Left: Convergence plots of LM-CMA-ES, sep-CMA-ES and Cholesky-CMA-ES on 128-dimensional axis-parallel and rotated Ellipsoid functions. Right: The median of 11 runs on separable Ellipsoid function for different problem dimensions. The dotted lines correspond to extrapolated results by preserving the same scaling as between the last two actual estimations.

Figure 4 :

 4 Figure 4: LM-CMA-ES and sep-CMA-ES on separable 100,000-dimensional Ellipsoid problem. The sep-CMA-ES divergences after the first generation (the best fitness is shown). Note that the LM-CMA-ES is rotationally invariant, therefore a similar performance is expected on 100,000-dimensional rotated Ellipsoid.

more precisely, we mean m evolution paths p c and their inverses v but for brevity we say m direction vectors