Ilya Loshchilov
email: ilya.loshchilov@gmail.com

LM-CMA: an Alternative to L-BFGS for Large Scale Black-box Optimization

Keywords: LM-CMA, L-BFGS, CMA-ES, large scale optimization, black-box optimization. Evolutionary Computation

 is often considered to be the method of choice for continuous optimization when first-and/or second-order information is available. However, the use of L-BFGS can be complicated in a black-box scenario where gradient information is not available and therefore should be numerically estimated. The accuracy of this estimation, obtained by finite difference methods, is often problem-dependent that may lead to premature convergence of the algorithm.

). The LM-CMA is a stochastic derivative-free algorithm for numerical optimization of non-linear, non-convex optimization problems. Inspired by the L-BFGS, the LM-CMA 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 memory complexity to O(mn), where n is the number of decision variables. The time complexity of sampling one candidate solution is also O(mn), but scales as only about 25 scalar-vector multiplications in practice. The algorithm has an important property of invariance w.r.t. strictly increasing transformations of the objective function, such transformations do not compromise its ability to approach the optimum. The LM-CMA outperforms the original CMA-ES and its large scale versions on non-separable ill-conditioned problems with a factor increasing with problem dimension. Invariance properties of the algorithm do not prevent it from demonstrating a comparable performance to L-BFGS on non-trivial large scale smooth and nonsmooth optimization problems.

Introduction

In a black-box scenario, knowledge about an objective function f : X → R, to be optimized on some space X, is restricted to the handling of a device that delivers the value of f (x) for any input x ∈ X. The goal of black-box optimization is to find solutions with small (in the case of minimization) value f (x), using the least number of calls to the function f [START_REF] Ollivier | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]. In continuous domain, f is defined as a mapping R n → R, where n is the number of variables. The increasing typical number of variables involved in everyday optimization problems makes it harder to supply the search with useful problem-specific knowledge, e.g., gradient information, valid assumptions about problem properties. The use of large scale black-box optimization approaches would seem attractive providing that a comparable performance can be achieved.

to appear in Evolutionary Computation journal, MIT Press

The use of well recognized gradient-based approaches such as the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [START_REF] Shanno | Conditioning of Quasi-Newton Methods for Function Minimization[END_REF] is complicated in the blackbox scenario since gradient information is not available and therefore should be estimated by costly finite difference methods (e.g., n + 1 function evaluations per gradient estimation for forward difference and 2n + 1 for central difference). The latter procedures are problem-sensitive and may require a priori knowledge about the problem at hand, e.g., scaling of f , decision variables and expected condition number [START_REF] Li | Effective condition number for finite difference method[END_REF].

By the 1980s, another difficulty has become evident: the use of quasi-Newton methods such as BFGS is limited to small and medium scale optimization problems for which the approximate inverse Hessian matrix can be stored in memory. As a solution, it was proposed not to store the matrix but to reconstruct it using information from the last m iterations [START_REF] Nocedal | Updating quasi-newton matrices with limited storage[END_REF]. The final algorithm called the limited memory BFGS algorithm (L-BFGS or LM-BFGS) proposed by [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] is still considered to be the state-of-the-art of large scale gradient-based optimization [START_REF] Becker | A quasi-newton proximal splitting method[END_REF]. However, when a large scale black-box function is considered, the L-BFGS is forced to deal both with a scarce information coming from only m recent gradients and potentially numerically imprecise estimations of these gradients which scale up the run-time in the number of function evaluations by a factor of n. It is reasonable to wonder whether the L-BFGS and other derivative-based algorithms are still competitive in these settings or better performance and robustness can be achieved by derivative-free algorithms.

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) seems to be a reasonable alternative, it is a derivative-free algorithm 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 the Black-Box Optimization Benchmarking (BBOB) workshop [START_REF] Finck | Real-parameter black-box optimization benchmarking 2010: Experimental setup[END_REF][START_REF] Auger | BBOB 2010: 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]Loshchilov, 2013a). 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] Kern | Local Meta-Models for Optimization Using Evolution Strategies[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 CMA-ES, the learning of dependencies between n decision variables, also forms its main practical limitations such as O(n 2 ) memory storage required to run the algorithm and O(n 2 ) computational time complexity per function evaluation [START_REF] Hansen | A simple modification in CMA-ES achieving and space complexity[END_REF]. These limitations may preclude the use of CMA-ES for computationally cheap but large scale optimization problems if the internal computational cost of CMA-ES is greater than the cost of one function evaluation. On nontrivial large scale problems with n > 10, 000 not only the internal computational cost of CMA-ES becomes enormous but it is becoming simply impossible to efficiently store the covariance matrix in memory. An open problem is how to extend efficient black-box approaches such as CMA-ES to n ≫ 1000 while keeping a reasonable trade-off between the performance in terms of the number of function evaluations and the internal time and space complexity. The low complexity methods such as separable CMA-ES (sep-CMA-ES by [START_REF] Hansen | A simple modification in CMA-ES achieving and space complexity[END_REF]), linear time Natural Evolution Strategy (R1-NES by [START_REF] Sun | A linear time natural evolution strategy for non-separable functions[END_REF]) and VD-CMA by [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF] are useful when the large scale optimization problem at hand is separable or decision variables are weakly correlated, otherwise the performance of these algorithms w.r.t. the original CMA-ES may deteriorate significantly.

In this paper, we present a greatly improved version of the recently proposed extension of CMA-ES to large scale optimization called the limited memory CMA-ES (LM-CMA) by [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF]. Instead of storing the covariance matrix, the LM-CMA stores m direction vectors in memory and uses them to generate solutions. The algorithm has O(mn) space complexity, where m can be a function of n. The time complexity is linear in practice with the smallest constant factor among the presented evolutionary algorithms.

The paper is organized as follows. First, we briefly describe L-BFGS in Section 2 and CMA-ES with its large scale alternatives in Section 3. Then, we present the improved version of LM-CMA in Section 4, investigate its performance w.r.t. large scale alternatives in Section 5 and conclude the paper in Section 6.

The L-BFGS

An early version of the L-BFGS method, at that time called the SQN method, was proposed by [START_REF] Nocedal | Updating quasi-newton matrices with limited storage[END_REF]. During the first m iterations the L-BFGS is identical to the BFGS method, but stores BFGS corrections separately until the maximum number of them m is used up. Then, the oldest corrections are replaced by the newest ones. The approximate of the inverse Hessian of f at iteration k, H k is obtained by applying m BFGS updates to a sparse symmetric and positive definite matrix H 0 provided by the user [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF].

Let us denote iterates by x k , s k = x k+1x k and y k = g k+1g k , where g denotes gradient. The method uses the inverse BFGS formula in the form

H k+1 = V T k H k V k + ρ k s k s T k , (1)
where ρ k = 1/y T k s k , and [START_REF] Dennis | Numerical methods for unconstrained optimization and nonlinear equations[END_REF][START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF].

V k = I -ρ k y k s T k
The L-BFGS method works as follows [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF]:

Step 1. Choose x 0 , m, 0 < β ′ < 1/2, β ′ < β < 1, and a symmetric and positive definite starting matrix H 0 . Set k = 0.

Step 2. Compute

d k = -H k g k , (2)
x k+1 = x k + α k d k , (3)
where α k satisfies the Wolfe conditions [START_REF] Wolfe | Convergence conditions for ascent methods[END_REF]:

f (x k + α k d k ) ≤ f (x k ) + β ′ α k g T k d k , (4)
g(x k + α k d k ) T d k ≥ βg T k d k . (5)
The novelty introduced by [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] w.r.t. the version given in [START_REF] Nocedal | Updating quasi-newton matrices with limited storage[END_REF] is that the line search is not forced to perform at least one cubic interpolation, but the unit steplength α k = 1 is always tried first, and if it satisfies the Wolfe conditions, it is accepted.

Evolutionary Computation

Step 3. Let m = min(k, m-1). Update H 0 m+1 times using the pairs y j , s j k j=k-m as follows:

H k+1 = (V T k • • • V T k-m)H 0 (V k-m • • • V k ) +ρ k-m(V T k • • • V T k-m+1 )s k-ms T k-m(V k-m+1 • • • V k ) +ρ k-m+1 (V T k • • • V T k-m+2 )s k-m+1 s T k-m+1 (V k-m+2 • • • V k ) . . . +ρ k s k s T k
Step 4. Set k = k + 1 and go to Step 2. The algorithm space and time complexity scales as O(mn) per iteration (not per function evaluation), where m in order of 5-40 suggested in the original paper is still the most common setting. An extension to bound constrained optimization called L-BFGS-B has the efficiency of the original algorithm, however at the cost of a significantly more complex implementation [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. Extensions to optimization with arbitrary constraints are currently not available. Satisfactory and computationally tractable handling of noise is at least problematic, often impossible.

Nevertheless, as already mentioned above, when gradient information is available, L-BFGS is competitive to other techniques [START_REF] Becker | A quasi-newton proximal splitting method[END_REF][START_REF] Ngiam | On optimization methods for deep learning[END_REF] and often can be viewed as a method of choice [START_REF] Andrew | Scalable training of l 1-regularized log-linear models[END_REF] for large scale continuous optimization. However, in the black-box scenario when gradientinformation is not available (direct search settings), the advantages of L-BFGS are becoming less obvious and derivative-free algorithms can potentially perform comparable. In this paper, we investigate this scenario in detail.

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 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]. More specifically, the analysis of theoretical foundations of Evolution Strategies is provided by [START_REF] Wierstra | Natural evolution strategies[END_REF]; [START_REF] Ollivier | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]; [START_REF] Akimoto | Objective improvement in information-geometric optimization[END_REF]; [START_REF] Glasmachers | Convergence of the IGO-Flow of isotropic gaussian distributions on convex quadratic problems[END_REF]; [START_REF] Auger | Linear Convergence of Comparison-based Step-size Adaptive Randomized Search via Stability of Markov Chains[END_REF]; Hansen and Auger (2014); [START_REF] Arnold | On the behaviour of the (1,λ)-ES for conically constrained linear problems[END_REF]; [START_REF] Beyer | Convergence analysis of evolutionary algorithms that are based on the paradigm of information geometry[END_REF].

The CMA-ES

The Covariance Matrix Adaptation Evolution Strategy (Hansen andOstermeier, 1996, 2001;[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 Algorithm 1 The (µ/µ w , λ)-CMA-ES

1: given n ∈ N + , λ = 4 + ⌊3 ln n⌋, µ = ⌊λ/2⌋, w i = ln(µ+ 1 2 )-ln i µ j=1 (ln(µ+ 1 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), c c = 4 n+4 , c 1 = 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 w i x i:λ // the symbol i : λ denotes i-th best individual on f 8: p t+1 σ ← (1 -c σ )p t σ + c σ (2 -c σ ) √ µ w C 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 -c c )p t c + h σ c c (2 -c c ) √ µ w m t+1 -m t σ t 11: C µ ← µ i=1 w i x i:λ -m t σ t × (x i:λ -m t ) T σ t 12: C t+1 ← (1 -c 1 -c µ )C t + c 1 p t+1 c p t+1 c T rank-one update +c µ C µ rank-µ update 13: σ t+1 ← σ t exp( cσ dσ ( p t+1 σ E N(0,I ) -1))
14:

t ← t + 1 15: until stopping criterion is met C t ∈ R n×n and a mutation step-size σ t as follows:

x t k ← N m t , σ t 2 C t ← m t + σ t N 0, C t (6)
These λ solutions then should be evaluated with 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 the selected individuals, the 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 the 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). The 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 rank-one update computes evolution path p t+1 c of successful moves of the mean m t+1 -m t σ t of the Evolutionary Computation mutation distribution in the given coordinate system (line 10), in a similar way as for 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).

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 the previously accumulated information by a new one with corresponding weights of importance (line 12): c 1 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 c 1 + c µ ≤ 1. It was also proposed to 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].

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 O(n 2 ) time and space complexity of CMA-ES precludes its applications to large scale optimization with n ≫ 1000. To enable the algorithm for large scale optimization, a linear time and space version called sep-CMA-ES was proposed by Ros and [START_REF] Hansen | Adaptive encoding: How to render search coordinate system invariant[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 -c cov )c t jj + 1 µ cov p t+1 c 2 j + c ccov 1 - 1 µ ccov µ i=1 w i c t jj z i:λ t+1 2 j , j = 1, . . . , n (7)
where, for j = 1, . . . , n the c jj are the diagonal elements of C t and the z i:λ t+1 j = x i:λ t+1 j /(σ t c t jj ). This update reduces the computational complexity to O(n) and allows to exploit problem separability. The algorithm demonstrated good performance on separable problems and even outperformed CMA-ES on non-separable Rosenbrock function for n > 100.

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

C = σ 2 (I + uu T ), (8)
where u and σ are the parameters to be adjusted. The adaptation of the predominant eigendirection u allows the algorithm to solve highly non-separable problems while maintaining only O(n) time and O(µn) space complexity. The use of the natural gradient in the derivation of the NES algorithm motivated a research which led to the formulation of the Information Geometric Optimization (IGO) framework by [START_REF] Ollivier | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF].

The IGO framework was used to derive a similar to R1-NES algorithm called VD-CMA [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF] with the sampling distribution parametrized by a Gaussian model with the covariance matrix restricted as follows:

C = D(I + uu T )D, ( 9
)
where D is a diagonal matrix of dimension n and u is a vector in R n . This model is able to represent a scaling for each variable by D and a principal component, which is generally not parallel to an axis, by Dv [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF]. It has O(n) time and O(µn) space complexity but i) typically demonstrates a better performance than sep-CMA-ES and R1-NES and ii) can solve a larger class of functions [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF].

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 eigenvectors and eigenvalues spanning the m-dimensional dominant subspace of the n × n-dimensional covariance matrix C . The authors adapted a singular value decomposition updating algorithm developed by [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 the number of function evaluations gradually decreases while enabling the search in R n for n > 10, 000. However, the computational complexity of O(m 2 n) practically (for m in order of √ n as suggested by [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF])

leads to the same limitations of O(n 2 ) time complexity as in the original CMA-ES. The (µ/µ w , λ)-Cholesky-CMA-ES proposed by [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 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 , (10)
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 , (11)
for z t = 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 ], (12)
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 proceeds by sampling λ candidate solutions (lines 4 -7) and taking into account the most successful µ out of λ solutions in the evolution path adaptation (lines 10 and 11). However, the eigendecomposition 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 keeps its time complexity as 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 Evolutionary Computation Algorithm 2 The (µ/µ w , λ)-Cholesky-CMA-ES

1: given n ∈ N + , λ = 4 + ⌊3 ln n⌋, µ = ⌊λ/2⌋, w i = 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), c c = 4 n+4 , c 1 = 2 (n+ √ 2) 2 2: initialize m t=0 ∈ R n , σ t=0 > 0, p t=-1 σ = 0, p t=-1 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 w i x i:λ 9: z w ← µ i=1 w i z i:λ 10: p t σ ← (1 -c σ )p t-1 σ + c σ (2 -c σ ) √ µ w z w 11: p t c ← (1 -c c )p t-1 c + c c (2 -c c ) √ µ w A t z w 12: v t ← A t inv p t c 13: A t+1 ← √ 1 -c 1 A t + √ 1-c1 v t 2 1 + c1 1-c1 v t 2 -1 p t c v t T 14: A t+1 inv ← 1 √ 1-c1 A t inv - 1 √ 1-c1 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 σ E N(0,I ) -1))
16:

t ← t + 1 17: until stopping criterion is met original CMA-ES) because the quadratic complexity will remain due to matrix-vector 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

The LM-CMA is inspired by the L-BFGS algorithm but instead of storing m gradients for performing inverse Hessian requiring operations it stores m direction vectors to reproduce the Cholesky factor A and generate candidate solutions with a limited time and space cost O(mn) (see Section 4.1). These m vectors are estimates of descent directions provided by evolution path vectors and should be stored with a particular temporal distance to enrich A (see Section 4.2). An important novelty introduced w.r.t. the original LM-CMA proposed by [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF] is a procedure for sampling from a family of search representations defined by Cholesky factors reconstructed from m * ≤ m vectors (see Section 4.3) and according to the Rademacher distribution (see Section 4.4). These novelties allow to simultaneously reduce the internal time complexity of the algorithm and improve its performance in terms of the number of function evaluations. Before describing the algorithm itself, we gradually introduce all the necessary components.

Algorithm 3 Az(): Cholesky factor -vector update

1: given z ∈ R n , m ∈ Z + , j ∈ Z m + , i ∈ Z |i| + , P ∈ R m×n , V ∈ R m×n , b ∈ R m , a ∈ [0, 1] 2: initialize x ← z 3: for t ← 1, . . . , |i| do 4: k ← b j i t V (j i t ,:) • z 5:
x ← ax + kP (j i t ,:) 6: return x Algorithm 4 Ainvz(): inverse Cholesky factor -vector update

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

Reconstruction of Cholesky Factors

By setting a = √ 1 -c 1 , b t = √ 1-c1 v t 2 1 + c1 1-c1 v t 2
-1 and considering the evolution path p t c (a change of optimum estimate m smoothed over iterations, see line 12 of Algorithm 7) together with v t = A -1 t p t c , one can rewrite Equation ( 11) as

A t+1 = aA t + b t p t c v t T , (13)
The product of a random vector z and the Cholesky factor A t thus can be directly computed. At iteration t = 0, A 0 = I and A 0 z = z, the new updated Cholesky factor

A 1 = aI + b 0 p 0 c v 0 T . At iteration 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 vectors 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 t -weighted evolution path p t c (accessed from a matrix P ∈ R m×n as P (it,:) ) scaled by the dot product of v t and x. As can be seen, Algorithms 3 and 4 use j i t indexation instead of t. This is a convenient way to have references to matrices P and V which store p c and v vectors, respectively. In the next subsections, 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 -c 1 and d t = 1 √ 1-c1 v t 2 × 1 - 1 1+ c 1 1-c 1 v t 2 .

The computational complexity of both procedures scales as O(mn).

It is important to note that when a vector p ℓ from a set of m vectors stored in P is replaced by a new vector p t+1 (see line 15 in Algorithm 7), all inverse vectors v t for t = ℓ, . . . , m should be iteratively recomputed [START_REF] Krause | [END_REF]. This procedure corresponds to line 16 in Algorithm 7.

Evolutionary Computation

Algorithm 5 UpdateSet(): direction vectors update

1: given m ∈ R + , j ∈ Z m + , l ∈ Z m + , t ∈ Z + , N ∈ Z m + , P ∈ R m×n , p c ∈ R n , T ∈ Z + 2: t ← ⌊t/T ⌋ 3: if t ≤ m then 4: j t ← t 5: else 6: i min ← 1 + argmin i l j i+1 -l j i -N i , |1 ≤ i ≤ (m -1) 7: if l j i min -l j i min -1 -N i ≥ 0 then 8: i min ← 1 9: j tmp ← j imin 10:
for i ← i min , . . . , m -1 do 11:

j i ← j i+1 12:
j m ← j tmp 13: j cur ← j min(t+1,m) 14: l jcur ← tT 15: P (j cur ,:) ← p c 16: return: j cur , j, l

Direction Vectors Selection and Storage

The choice to store only m ≪ n direction vectors p c to obtain a comparable amount of useful information as stored in the covariance matrix of the original CMA-ES requires a careful procedure of selection and storage. A simple yet powerful procedure proposed by [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF] is to preserve a certain temporal distance in terms of number of iterations between the stored direction vectors. The procedure tends to store a more unique information in contrast to the case if the latest m evolution path vectors would be stored. The latter case is different from the storage of m gradients as in L-BFGS since the evolution path is gradually updated at each iteration with a relatively small learning rate c c and from µ ≪ n selected vectors.

The 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 j cur points out to the 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 (P (j i min -1 ,:) , P

(j i min ,:
) ) with the distance between them (in terms of indexes of iterations, stored in l) closest to a target distance N i (line 6). If this distance is smaller than N i then the index j imin will be swapped with last index of j (lines 9-12) and the corresponding vector P (j i min

,:) will be replaced by the new vector p c (line 15), otherwise the oldest vector among m saved vectors will be removed (as a result of line 8). Thus, the procedure gradually replaces vectors in a way to keep them at the distance N i and with the overall time horizon for all vectors of at most m-1 i N i iterations. The procedure can be called periodically every T ∈ Z + iterations of the algorithm. The values of N i are to be defined, e.g., as a function of problem dimension n and direction vector index i. Here, however, we set N i to n for all i, i.e., the target distance equals to the problem dimension.

Algorithm 6 SelectSubset(): direction vectors selection

1: given m ∈ Z + , m σ = 4, k ∈ Z + 2: if k = 1 then 3: m σ ← 10m σ 4: m * ← min(⌊m σ |N (0, 1)|⌋ , m) 5: i ← (m + 1 -m * ), . . . , m 6: return i

Sampling from a Family of Search Representations

At iteration t, a new k-th solution can be generated as

x k ← m t + σ t Az(z k , i), (14)
where z k ∈ R n is a vector drawn from some distribution and transformed by a Cholesky factor by calling Az(z k , i). The Az() procedure (see Algorithm 3) has an input i which defines indexes of direction vectors used to reconstruct the Cholesky factor. It is important to note that P (1,:) refers to the first vector physically stored in matrix P, P (j 1 ,:) refers to the oldest vector, P (j i t ,:) refers to the i t -th oldest vector according to an array i with indexes of vectors of interest. Thus, by setting i = 1, . . . , m all m vectors will be used in the reconstruction. Importantly, omission of some vector in i can be viewed as setting of the corresponding learning rate in Equation ( 13) to 0.

By varying i, one can control the reconstruction of the Cholesky factor used for sampling and in this way explore a family of possible transformations of the coordinate system. The maximum number of vectors defined by m can be associated with the number of degrees of freedom of this exploration.

While in the original LM-CMA [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF] the value of m is set to 4 + ⌊3 ln n⌋ to allow the algorithm scale up to millions of variables, we found that greater values of m, e.g., √ n often lead to better performance (see Section 5.7 for a detailed analysis). Thus, when memory allows, a gain in performance can be achieved. However, due to an internal cost O(mn) of Az(), the time cost then would scale as O(n 3/2 ) which is undesirable for n ≫ 1000. This is where the use of m * out of m vectors can drastically reduce the time complexity. We propose to sample m * from a truncated half-normal distribution N 0, m 2 σ (see line 4 of Algorithm 6) and set i to the latest m * vectors (line 5). For a constant m σ = 4, the time complexity of Az() scales as O(n). New value of m * is generated for each new individual. Additionally, to exploit the oldest information, we force m * to be generated with 10m σ for one out of λ individuals. While for m * = 0 the new solution x k appears to be sampled from an isotropic normal distribution, the computation of v inverses is performed using all m vectors.

Sampling from the Rademacher Distribution

Evolution Strategies are mainly associated with the multivariate normal distribution used to sample candidate solutions. However, alternative distributions such as the Cauchy distribution can be used [START_REF] Yao | Fast evolution strategies[END_REF][START_REF] Hansen | Adaptive encoding: How to render search coordinate system invariant[END_REF]. Moreover, the Adaptive Encoding procedure proposed by [START_REF] Hansen | Adaptive encoding: How to render search coordinate system invariant[END_REF] can be coupled with any sampling distribution as in [START_REF] Loshchilov | Adaptive Coordinate Descent[END_REF], where it was shown that completely deterministic adaptive coordinate descent on principal components obtained with the Adaptive Encoding procedure allows to obtain the performance comparable to the one of CMA-ES.

Evolutionary Computation

In this paper, inspired by [START_REF] Loshchilov | Adaptive Coordinate Descent[END_REF], we replace the original multivariate normal distribution used in LM-CMA by the Rademacher distribution, where a random variable has 50% chance to be either -1 or +1 (also can be viewed as a Bernoulli distribution). Thus, a pre-image vector of candidate solution z contains n values which are either -1 or +1. Our intention to use this distribution is three-fold: i) to reduce the computation complexity by a constant but rather significant factor, ii) to demonstrate that the Rademacher distribution can potentially be an alternative to the Gaussian distribution at least in large scale settings, iii) to demonstrate that our new step-size adaptation rule (see next section), which does not make assumptions about the sampling distribution, can work well when used with non-Gaussian distributions. As a support for this substitution, we recall that for a n-dimensional unit-variance spherical Gaussian, for any positive real number β ≤ √ n, all but at most 3 exp -cβ 2 of the mass lies within the annulus

√ n -1 -β ≤ r ≤ √ n -1 + β,
where, c is a fixed positive constant [START_REF] Hopcroft | Foundations of Data Science[END_REF]. Thus, when n is large, the mass is concentrated in a thin annulus of width O(1) at radius √ n. Interestingly, the sampling from the Rademacher distribution reproduces this effect of large-dimensional Gaussian sampling since the distance from the center of a n-dimensional hypercube to its corners is √ n.

Population Success Rule

The step-size used to define the scale of deviation of a sampled candidate solution from the mean of the mutation distribution can be adapted by the Population Success Rule (PSR) proposed for LM-CMA by [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF]. This procedure does not assume that candidate solutions should come from the multivariate normal distribution as it is often assumed in Evolution Strategies including CMA-ES. Therefore, PSR procedure is well suited for the Rademacher distribution.

The PSR is inspired by the median success rule [START_REF] Ait Elhara | A median success rule for non-elitist evolution strategies: Study of feasibility[END_REF]. 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 (15)
Then, all individuals in the mixed set are ranked to define two sets r t-1 and r t (the lower the rank the better the individual) containing ranks of individuals of the previous and current populations ranked in the mixed set. A normalized success measurement is computed as

z P SR ← λ i=1 r t-1 (i) -r t (i) λ 2 -z * , (16)
where z * is a target success ratio and λ 2 accounts for the normalization of the sum term and for different possible population size λ. Then, for s ← (1c σ )s + c σ z P SR , the step-size is adapted as

σ t+1 ← σ t exp (s/ d σ ), (17)
where d σ is a damping factor which we set here to 1.

The Algorithm

The improved LM-CMA is given in Algorithm 7. At each iteration t, λ candidate solutions are generated by mutation defined as a product of a vector z k sampled from the Rademacher distribution and a Cholesky factor A t reconstructed from m * out of Algorithm 7 The (µ/µ w , λ)-LM-CMA

1: given n ∈ N + , λ = 4 + ⌊3 ln n⌋, µ = ⌊λ/2⌋, w i = ln(µ+1)-ln(i) µ ln(µ+1)-µ j=1 ln(j) ; i = 1 . . . µ, µ w = 1 µ i=1 w 2 i , c σ = 0.3, z * = 0.25, m = 4 + ⌊3 ln n⌋, N steps = n, c c = 0.5 √ n , c 1 = 1 10 ln(n+1) , d σ = 1, T = ⌊log(n)⌋ 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: if k (mod 2) = 1 then 6: z k ← Rademacher() 7: i ← SelectSubset(k) 8: x k ← m t + σ t Az(z k , i) 9: else 10: x k ← m t -(x k-1 -m t ) 11: f t k ← f (x k ) 12: m t+1 ← µ i=1 w i x i:λ 13: p t+1 c ← (1 -c c )p t c + c c (2 -c c ) √ µ w (m t+1 -m t )/σ t 14: if t (mod T ) = 0 then 15:
U pdateSet(p t+1 c )

16:

U pdateInverses()

17: r t , r t-1 ← Ranks of f t and f t-1 in f t ∪ f t-1 18: z P SR ← λ i=1 r t-1 (i)-r t (i) λ 2 -z * 19: s ← (1 -c σ )s + c σ z P SR 20: σ t+1 ← σ t exp(s/d σ )
21:

t ← t + 1 22: until stopping criterion is met m vectors (line 4-11) as described in Sections 4.1-4.4. We introduce the mirrored sampling [START_REF] Brockhoff | Mirrored sampling and sequential selection for evolution strategies[END_REF] to generate the actual x k only every second time and thus decrease the computation cost per function evaluation by a factor of two by evaluating m t + σ t Az(z k ) and then its mirrored version m t -(x k-1m t ). The latter approach sometimes also improves the convergence rate.

The best µ out of λ solutions are selected to compute the new mean m t+1 of the mutation distribution in line 12. The new evolution path p t+1 c is updated (line 13) from the change of the mean vector √ µ w (m t+1 -m t )/σ t and represents an estimate of descent direction. One vector among m vectors is selected and replaced by the new p t+1 c in UpdateSet() procedure described in Section 4.2. All inverses v of evolution path vectors which are at least as recent as the direction vector to be replaced should be recomputed in the UpdateInverses() procedure as described in Section 4.1. The step-size is updated according to the PSR rule described in Section 4.5.

Experimental Validation

The performance of the LM-CMA is investigated comparatively to the L-BFGS [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF], the active CMA-ES by [START_REF] Hansen | Benchmarking a weighted negative covariance matrix update on the BBOB-2010 noiseless testbed[END_REF] and the VD-CMA by [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF]. The sep-CMA-ES is removed from the comparison due to its similar but inferior performance w.r.t. the VD-CMA observed both in our study and by [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF].

Evolutionary Computation

Table 1: Test functions, initialization intervals and initial standard deviation (when applied). R is an orthogonal n × n matrix with each column vector q i being a uniformly distributed unit vector implementing an angle-preserving transformation (Ros and [START_REF] Hansen | Adaptive encoding: How to render search coordinate system invariant[END_REF] .

Name Function Target f (x) Init σ 0 Sphere f Sphere (x)= n i=1 x 2 i 10 -10 [-5, 5] n 3 Ellipsoid f Elli (x)= n i=1 10 6 i-1 n-1 x 2 i 10 -10 [-5, 5] n 3 Rosenbrock fRosen(x)= n-1 i=1 100.(x 2 i -xi+1) 2 + (xi -1) 2 10 -10 [-5, 5] n 3 Discus fDiscus(x)= 10 6 x 2 1 + n i=2 x 2 i 10 -10 [-5, 5] n 3 Cigar fCigar (x)= x 2 1 + 10 6 n i=2 x 2 i 10 -10 [-5, 5] n 3 Different Powers f Dif f P ow (x)= n i=1 |xi| 2+4(i-1)/(n-1) 10 -10 [-5, 5] n 3 Rotated Ellipsoid f RotElli (x)=f Elli (Rx) 10 -10 [-5, 5] n 3 Rotated Rosenbrock fRotRosen(x)=fRosen(Rx) 10 -10 [-5, 5] n 3 Rotated Discus fRotDiscus(x)=fDiscus (Rx) 10 -10 [-5, 5] n 3 Rotated Cigar fRotCigar (x)=fCigar (Rx) 10 -10 [-5, 5] n 3 Rotated Different Powers f RotDif f P ow (x)=f Dif f P ow (Rx) 10 -10 [-5, 5] n 3
We use the L-BFGS implemented in MinFunc library by [START_REF] Schmidt | Minfunc. Technical report[END_REF] in its default parameter settings1 , the active CMA-ES (aCMA) without restarts in its default parametrization of CMA-ES MATLAB code version 3.612 . For faster performance in terms of CPU time, the VD-CMA was (exactly) reimplemented in C language from the MATLAB code provided by the authors. For the sake of reproducibility, the source code of all algorithms is available online3 . The default parameters of LM-CMA are given in Algorithm 7.

We use a set of benchmark problems (see Table 1) commonly used in Evolutionary Computation, more specifically in the BBOB framework [START_REF] Finck | Real-parameter black-box optimization benchmarking 2010: Experimental setup[END_REF]. Indeed, many problems are missing including the ones where tested methods and LM-CMA fail to timely demonstrate reasonable performance in large scale settings. We focus on algorithm performance w.r.t. both the number of function evaluations used to reach target values of f , CPU time spent per function evaluation and the number of memory slots required to run algorithms. Any subset of these metrics can dominate search cost in large scale settings, while in low scale settings memory is typically of a lesser importance.

In this section, we first investigate the scalability of the proposed algorithm w.r.t. the existing alternatives. While both the computational time and space complexities scale moderately with problem dimension, the algorithm is capable to preserve certain invariance properties of the original CMA-ES. Moreover, we obtain unexpectedly good results on some well-known benchmark problems, e.g., linear scaling of the budget of function evaluations to solve Separable and Rotated Ellipsoid problems. We demonstrate that the performance of LM-CMA is comparable to the one of L-BFGS with exact estimation of gradient information. Importantly, we show that LM-CMA is competitive to L-BFGS in very large scale (for derivative-free optimization) settings with 100,000 variables. Finally, we investigate the sensitivity of LM-CMA to some key internal parameters such as the number of stored direction vectors m.

Space Complexity

The number of memory slots required to run optimization algorithms versus problem dimension n referred to as space complexity can limit applicability of certain algorithms to large scale optimization problems. Here, we list the number of slots up to constant and asymptotically constant terms.

The presented algorithms store λ generated solutions (LM-CMA, VD-CMA and aCMA with the default λ = 4 + ⌊3 ln n⌋ and µ = λ/2, L-BFGS with λ = 1) and some intermediate information (LM-CMA and L-BFGS with m pairs of vectors, aCMA with at least two matrices of size n × n) to perform the search. Our implementation of VD-CMA requires (max(µ + 15, λ) + 7)n slots compared to (2.5λ + 21)n of the original MATLAB code. The LM-CMA requires (2m + λ + 6)n + 5m slots, the L-BFGS requires (2m + 3)n slots and aCMA requires (2n + λ + 3)n slots.

Figure 1 shows that due to its quadratic space complexity aCMA requires about 2 × 10 8 slots (respectively, 8 × 10 8 slots) for 10,000-dimensional (respectively, 20,000dimensional) problems which with 8 bytes per double-precision floating point number would correspond to about 1.6 GB (respectively, 6.4 GB) of computer memory. This simply precludes the use of CMA-ES and its variants with explicit storage of the full covariance matrix or Cholesky factors to large scale optimization problems with n > 10, 000. LM-CMA stores m pairs of vectors as well as the L-BFGS. For m = 4 + ⌊3 ln n⌋ (as the default population size in CMA-ES), L-BFGS is 2 times and LM-CMA is 3 times more expensive in memory than VD-CMA, but they all basically can be run for millions of variables.

In this paper, we argue that additional memory can be used while it is allowed and is at no cost. Thus, while the default m = 4 + ⌊3 ln n⌋, we suggest to use m = ⌊2 √ n⌋ if memory allows (see Section 5.7). In general, the user can provide a threshold on memory and if, e.g., by using m = ⌊2 √ n⌋ this memory threshold would be violated, the algorithm automatically reduces m to a feasible m f .

Time Complexity

The average amount of CPU time internally required by an algorithm per evaluation of some objective function f ∈ R n (not per algorithm iteration) referred to as time complexity also can limit applicability of certain algorithms to large scale optimization problems. They simply can be too expensive to run, e.g., much more expensive than to perform function evaluations.

Figure 2 shows how fast CPU time per evaluation scales for different operations measured on one 1.8 GHz processor of an Intel Core i7-4500U. Scalar-vector multiplication of a vector with n variables scales linearly with ca. 4 • 10 -10 n evaluation of the separable Ellipsoid function is 2-3 times more expensive if a temporary data is used. Sampling of n normally distributed variables scales as ca. 60 vectors-scalar multiplications that defines the cost of sampling of unique candidate solutions of many Evolution Strategies such as separable CMA-ES and VD-CMA. However, the sampling of variables according to the Rademacher distribution is about 10 times cheaper. The use of mirrored sampling also decreases the computational burden without worsening the convergence. Finally, the internal computation cost of LM-CMA scales linearly as about 25 scalar-vector multiplications per function evaluation. It is much faster than the lower bound for the original CMA-ES defined by one matrix-vector multiplication required to sample one individual. We observe that the cost of one matrix-vector multiplications costs about 2n scalar-vector multiplications, the overhead is probably due to access to matrix members.

The proposed version of LM-CMA is about 10 times faster internally than the orig-inal version by [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF] due to the use of mirrored sampling, the Rademacher sampling distribution and sampling with m * instead of m direction vectors both for m = 4 + ⌊3 ln n⌋ and m = ⌊2 √ n⌋. For 8192-dimensional problems it is about 1000 times faster internally than CMA-ES algorithms with the full covariance matrix update (the cost of Cholesky-CMA-ES is given in Loshchilov ( 2014)).

Invariance under Rank-preserving Transformations of the Objective Function

The LM-CMA belongs to a family of so-called comparison-based algorithms. The performance of these algorithms is unaffected by rank-preserving (strictly monotonically increasing) transformations of the objective function, e.g., whether the function f , f 3 or f × |f | -2/3 is minimized [START_REF] Ollivier | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]. Moreover, this invariance property provides robustness to noise as far as this noise does not impact a comparison of solutions of interest [START_REF] Auger | Linear Convergence of Comparison-based Step-size Adaptive Randomized Search via Stability of Markov Chains[END_REF].

In contrast, gradient-based algorithms are sensitive to rank-preserving transformations of f . While the availability of gradient information may mitigate the problem that objective functions with the same contours can be solved with a different number of functions evaluations, the lack of gradient information forces the user to estimate it with approaches whose numerical stability is subject to scaling of f . Here, we simulate L-BFGS in an idealistic black-box scenario when gradient information is estimated perfectly (we provide exact gradients) but at the cost of n + 1 function evaluations per gradient that corresponds to the cost of the forward difference method. Additionally, we investigate the performance of L-BFGS with the central difference method (2n + 1 evaluations per gradient) which is twice more expensive but numerically more stable. We denote this method as CL-BFGS.

Invariance under Search Space Transformations

Invariance properties under different search space transformations include translation invariance, scale invariance, rotational invariance and general linear invariance under any full rank matrix R when the algorithm performance on f (x) and f (Rx) is the same given that the initial conditions of the algorithm are chosen appropriately [START_REF] Hansen | Impacts of invariance in search: When CMA-ES and PSO face ill-conditioned and non-separable problems[END_REF]. Thus, the lack of the latter invariance is associated with a better algorithm performance for some R and worse for the others. In practice, it often appears to be relatively simple to design an algorithm specifically for a set of problems with a R, e.g., identity matrix, and then demonstrate its good performance. If this set contains separable problems, the problems where the optimum can be found with a coordinate-wise search, then even on highly multi-modal functions great results can be easily achieved [START_REF] Loshchilov | Bi-population CMA-ES agorithms with surrogate models and line searches[END_REF]. Many derivative-free search algorithms in one or another way exploit problem separability and fail to demonstrate a comparable performance on, e.g., rotated versions of the same problems. This would not be an issue if most of real-world problems are separable, this is, however, unlikely to be the case and some partial-separability or full non-separability are more likely to be present.

The original CMA-ES is invariant w.r.t. any invertible linear transformation of the search space, R, if the initial covariance matrix C t=0 = R -1 (R -1 ) T , and the initial search point(s) are transformed accordingly [START_REF] Hansen | The cma evolution strategy: a comparing review[END_REF]. However, R matrix is often unknown (otherwise, one could directly transform the objective function) and cannot be stored in memory in large scale settings with n ≫ 10, 000. Thus, the covariance matrix adapted by LM-CMA has at most rank m and so the intrinsic coordinate system cannot capture some full rank matrix R entirely. Therefore, the performance of the algorithm on f (Rx) compared to f (x) depends on R. However, in our experiments, differences in performance on axis-aligned and rotated ill-conditioned functions were marginal.

Evolutionary Computation

Here, we test LM-CMA, aCMA, L-BFGS, CL-BFGS both on separable problems and their rotated versions (see Table 1). It is simply intractable to run algorithms on large scale rotated problems with n > 1000 due to the quadratic cost of involved matrixvectors multiplications (see Figure 2). Fortunately, there is no need to do it for algorithms that are invariant to rotations of the search space since their performance is the same as on the separable problems whose evaluation is cheap (linear in time). Figures 34show that the performance of aCMA on 512-dimensional (the dimensionality still feasible to perform full runs of aCMA) separable (left column) and rotated (right column) problems is very similar and the difference (if any) is likely due to a non-invariant initialization. The invariance in performance is not guaranteed but rather observed for LM-CMA, L-BFGS and CL-BFGS. However, the performance of the VD-CMA degrades significantly on f RotElli , f RotDiscus and f RotDif f P ow functions due to the restricted form of the adapted covariance matrix of Equation ( 9). Both f Rosen , f Cigar and their rotated versions can be solved efficiently since they have a Hessian matrix whose inverse can be well approximated by Equation ( 9) [START_REF] Akimoto | Comparison-based natural gradient optimization in high dimension[END_REF].

An important observation from Figures 34is that even the exact gradient infor- mation is not sufficient for L-BFGS to avoid numerical problems which lead to an imprecise estimation of the inverse Hessian matrix and premature convergence on f Elli and f RotElli . The L-BFGS with the central difference method (CL-BFGS) experiences an early triggering of stopping criteria on f Rosen and f Dif f P ow . While numerical problems due to imprecise derivative estimations are quite natural for L-BFGS especially on ill-conditioned problems, we assume that with a better implementation of the algorithm (e.g., with high-precision arithmetic) one could obtain a more stable convergence. Therefore, we extrapolate the convergence curves of L-BFGS and CL-BFGS towards the target f = 10 -10 after removing the part of the curve which clearly belongs to the stagnation, e.g., f < 10 -7 on f Elli .

Evolutionary Computation

Scaling with Problem Dimension

The performance versus the increasing number of problem variables is given in Figure 5. We exclude the results of VD-CMA on some problems because, as can be seen from Figures 34, the algorithm does not find the optimum with a reasonable number of function evaluations or/and convergences prematurely. For algorithms demonstrating the same performance on separable and non-separable problems (see Figures 34), we plot some results obtained on separable problems as obtained on rotated problems in Figure 5 to avoid possible misunderstanding from designers of separability-oriented algorithms.

The results suggest that L-BFGS shows the best performance, this is not surprising given the form of the selected objective functions (see Table 1). We should also keep in mind that the exact gradients were provided and this still led to premature convergence on some functions (see Figures 34). In the black-box scenario, one would probably use L-BFGS with the forward or central (CL-BFGS) difference methods. The latter is often found to lead to a loss by a factor of 2 (as expected due to 2n + 1 versus n + 1 cost per gradient), except for the f RotDif f P ow , where the loss is increasing with problem dimension.

Quite surprisingly, the LM-CMA outperforms VD-CMA and aCMA on f Sphere . This performance is close to the one obtained for (1+1) Evolution Strategy with optimal step-size. Bad performance on f Sphere anyway would not directly mean that an algorithm is useless, but could illustrate its performance in vicinity of local optima when variable-metric algorithms (e.g., CMA-like algorithms) may perform an isotropic search w.r.t. an adapted internal coordinate system. The obtained results are mainly due to the Population Success Rule which deserves an independent study similar to the one by Hansen et al. (2014). Nevertheless, we would like to mention a few key points of the PSR. By design, depending on the target success ratio z * , one can get either biased (for z * = 0) or unbiased (for z * = 0) random walk on random functions. It would be a bias to say that either biased or unbiased change of σ "is better" on random functions, since the latter depends on the context. Due to the fact that the (weighted) mean of each new population is computed from the best µ out of λ individuals, the λ individuals of the new generation are typically as good as the (weighted) best µ individuals of the previous one, and, thus, if z * = 0 one may expect z P SR > 0 from Equation ( 16). Typically, it is reasonable to choose z * ∈ (0, 0.5) lower-bounded by 0 due to random functions and upper-bounded by 0.5 due to linear functions. In this study, we choose 0.3 which lies roughly in the middle of the interval. It is important to mention a striking similarity with the 1/5th success rule [START_REF] Schumer | Adaptive step size random search. Automatic Control[END_REF][START_REF] Rechenberg | Evolutionsstrategie: optimierung technischer systeme nach prinzipien der biologischen evolution[END_REF]. We consider the PSR to be its population-based version.

The performance of LM-CMA on f Elli is probably the most surprising result of this work. In general, the scaling of CMA-ES is expected to be from super-linear to quadratic with n on f Elli since the number of parameters of the full covariance matrix to learn is (n 2 + n)/2 [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]. While aCMA demonstrates this scaling, LM-CMA scales linearly albeit with a significant constant factor. The performance of both algorithms coincides at n ≈ 1000, then, LM-CMA outperforms aCMA (given that our extrapolation is reasonable) with a factor increasing with n. It should be noted that aCMA is slower in terms of CPU time per function evaluation by a factor of n/10 (see Figure 2). Another interesting observation is that the L-BFGS is only slightly faster than LM-CMA, while CL-BFGS is actually outperformed by the latter. An insight to these observations can be found in Figure 3 where both LM-CMA and L-BFGS outperform aCMA by a factor of 10 in the initial part of the search, while aCMA compensates this loss by having the covariance matrix well adapted that allows to accelerate convergence close to the optimum. This might be explained as follows: a smaller number of internal parameters defining the intrinsic coordinate system can be learned faster and with greater learning rates, this allows a faster convergence but may slow down the search in vicinity of the optimum if the condition number cannot be captured by the reduced intrinsic coordinate system.

Evolutionary Computation

The LM-CMA is better or is as good as VD-CMA on f Rosen and f Cigar where it is expected to be outperformed by the latter due a presumably few principal components needed to be learnt to solve these problems. The scaling on f Rosen suggests that the problem is more difficult (e.g., more difficult than f Elli ) than one could expect, mainly due to an adaptation of the intrinsic coordinate system required while following the banana shape valley of this function.

The results on 100,000-dimensional problems (see Figure 6) show that LM-CMA outperforms L-BFGS on the first 10n -20n function evaluations which corresponds to the first 10-20 iterations of L-BFGS. This observation suggests that LM-CMA can be viewed as an alternative to L-BFGS when n is large and the available number of function evaluations is limited. While it can provide a competitive performance in the beginning, it is also able to learn dependencies between variables to approach the optimum.

Performance on a nonsmooth variant of Nesterov's function

While designed for smooth optimization, BFGS is known to work well for nonsmooth optimization too. A recent study by [START_REF] Overton | Nonsmooth, nonconvex optimization algorithms and examples[END_REF] demonstrated the difficulties encountered by BFGS on some nonsmooth functions. We selected one of the test functions from [START_REF] Overton | Nonsmooth, nonconvex optimization algorithms and examples[END_REF] called the second nonsmooth variant of Nesterov-Chebyshev-Rosenbrock function defined as follows:

N (x) = 1 4 |x 1 -1| + n-1 i=1 |x i+1 -2 |x i | + 1| (18)
This function is nonsmooth (though locally Lipschitz) as well as nonconvex, it has 2 n-1 Clarke stationary points [START_REF] Overton | Nonsmooth, nonconvex optimization algorithms and examples[END_REF]. [START_REF] Overton | Nonsmooth, nonconvex optimization algorithms and examples[END_REF] showed that for n = 5 BFGS starting from 1000 randomly generated points finds all 16 Clarke stationary points (for the definition of Clarke stationary points see [START_REF] Abramson | Convergence of mesh adaptive direct search to secondorder stationary points[END_REF]) and the probability to find the global minimizer is only by about a factor of 2 greater than to find any of the Clarke stationary points. This probability dropped by a factor of 2 for n = 6 while and since the number of Clarke stationary points doubled [START_REF] Overton | Nonsmooth, nonconvex optimization algorithms and examples[END_REF]. Clearly, the problem becomes extremely difficult for BFGS when n is large.

We launched LM-CMA and CL-BFGS (L-BFGS performed worse) on N (x) for n = 128 and n = 2048. Figure 7 shows that CL-BFGS performs better than LM-CMA, however, both algorithms in default settings and with restarts do not perform well. We tuned both LM-CMA and CL-BFGS but report the results only for LM-CMA since we failed to improve the performance of CL-BFGS by more than one order of magnitude of the objective function value. The tuned parameters for LM-CMA are: i) doubled population size λ, ii) increased learning rate by 15 to c 1 = 15/(10 ln(n + 1)), iii) an extremely small learning rate for step-size adaptation c σ = 0.3/n 2 instead of c σ = 0.3. The last modification is probably the most important, practically, it defines the schedule how step-size decreases. A similar effect can be achieved by reducing z * or increasing d σ . Faster learning of dependencies between variables and slower step-size decrease drastically improve the convergence and the problem can be solved both for n = 128 and n = 2048 (Figure 7). Interestingly, the number of function evaluations scales almost linearly with problem dimension.

We expected that tuning of CL-BFGS will lead to similar improvements. Surprisingly, our attempts to modify its parameters, often in order to slow down the convergence (e.g., type and number of line-search steps, Wolfe conditions parameters) failed. We still expect that certain modifications should improve CL-BFGS and thus we leave this question open. The settings tuned for N (x) function differ significantly from the default ones. It is of great interest to find an online procedure to adapt them. The next section is aimed at gaining some intuition on parameters importance in LM-CMA.

Evolutionary Computation

Rotated Ellipsoid

Sensitivity to Parameters

The black-box scenario implies that the optimization problem at hand is not known, it is therefore hard if even possible to suggest a "right" parametrization of our algorithm that works best on all problems. Offline tuning in large scale optimization is also computationally expensive. It is rather optimistic to believe that one always can afford enough computational resources to run algorithms till the optimum on very large realworld optimization problems. Nevertheless, we tend to focus on this scenario in order to gain an understanding about scalability on benchmark problems.

Our experience with parameter selection by exclusion of non-viable settings suggests that there exists a dependency between the population size λ, number of stored vectors m, the target temporal distance between them N steps , the learning rate c c for the evolution path and learning rate c 1 for the Cholesky factor update. The main reason for this is that all of them impact how well the intrinsic coordinate system defined by the Cholesky factor reflects the current optimization landscape. A posteriori, if m ≪ n, it seems reasonable to store vectors with a temporal distance in order of N steps = n on problems where a global coordinate system is expected to be constant, e.g., on a class Rosenbrock, m=5, 4866x10 3 evals. Evolutionary Computation of problems described by the general ellipsoid model [START_REF] Beyer | Convergence analysis of evolutionary algorithms that are based on the paradigm of information geometry[END_REF]. The learning rate for the evolution path is related to both m and n, here, we set it to c c = 0.5 √ n which is roughly inversely proportional to the (if affordable) suggested m = ⌊2 √ n⌋. We found that the chosen c c is still valid for the default m = 4 + ⌊3 ln n⌋. We do not have a good interpretation for the learning rate c 1 = 1 10 ln(n+1) . In general, we are not in favor of strongly arguing for some parameters settings against the others since as already mentioned above they are problem-dependent. A more appropriate approach would be to perform online adaptation of hyper-parameters as implemented for the original CMA-ES by [START_REF] Loshchilov | Maximum Likelihood-based Online Adaptation of Hyper-parameters in CMA-ES[END_REF].

We present an analysis for m which directly affects the amount of memory required to run the algorithm, and, thus, is of special interest since the user might be restricted in memory on very large scale optimization problems with n > 10 6 . Figure 8 shows that the greater the m the better the performance. The results obtained for the default m = 4+⌊3 ln n⌋, i.e., the results demonstrated in the previous sections can be improved with m = ⌊2 √ n⌋. The improvements are especially pronounced on f Discus functions, where the factor is increasing with n and the overall cost to solve the function reaches the one extrapolated for aCMA at n = 8192 (see Figure 5). It is surprising to observe that m = 5 and even m = 2 are sufficient to solve f Elli , f Discus and f Dif f P ow . The latter is not the case for f Cigar , where small values of m lead to an almost quadratic growth of run-time. The overall conclusion would be that on certain problems the choice of m is not critical, while greater values of m are preferable in general.

We investigated the eigenspectrum of the covariance matrix C t constructed as A t A t T from the Cholesky factor A t . The results for single runs on different 1024dimensional functions and for different m are shown in Figure 9. The evolution of the eigenspectrum during the run is shown by gradually darkening (blue) lines with increasing t. Clearly, the number of eigenvalues is determined by m. The profiles, e.g., the one of f Cigar , also reflect the structure of the problems (see Table 1). The greater the m, the greater condition number can be captured by the intrinsic coordinate system as can be see for f Elli , f Discus and f Dif f P ow , that in turn leads to a better performance. However, this is not always the case as can be seen for f Rosen that again demonstrates that optimal hyper-parameter settings are problem-dependent.

Conclusions

We adapt an idea from derivative-based optimization to extend best performing evolutionary algorithms such CMA-ES to large scale optimization. This allows to reduce the cost of optimization in terms of time by a factor of n/10 and memory by a factor between √ n and n. Importantly, it also often reduces the number of function evaluations required to find the optimum. The idea to store a limited number of vectors and use them to adapt an intrinsic coordinate system is not the only but one of probably very few ways to efficiently search in large scale continuous domains. We propose two quite similar alternatives: i) the storage of points and a later estimation of descent directions from differences of these points, and ii) the use of a reduced matrix m × n as in [START_REF] Knight | Reducing the space-time complexity of the CMA-ES[END_REF] but with a modified sampling procedure to obtain linear time complexity as proposed for the Adaptive Coordinate Descent by [START_REF] Loshchilov | Surrogate-Assisted Evolutionary Algorithms[END_REF].

The use of the Population Success Rule is rather optional and alternative step-size adaptation procedures can be applied. However, we find its similarity with the 1/5-th rule quite interesting. The procedure does not make any assumption about the sam-pling distribution, this allowed to use the Rademacher distribution. When n is large, the sampling from a n-dimensional Rademacher distribution resembles the sampling from a n-dimensional Gaussian distribution since the probability mass of the latter is concentrated in a thin annulus of width O(1) at radius √ n. The presented comparison shows that LM-CMA outperforms other evolutionary algorithms and is comparable to L-BFGS on non-trivial large scale optimization problems when the black-box (derivative-free) scenario is considered. Clearly, the black-box scenario is a pessimistic scenario but a substantial part of works that use finite difference methods for optimization deal with this scenario, and, thus, can consider LM-CMA as an alternative. Importantly, LM-CMA is invariant to rank-preserving transformations of the objective function and therefore is potentially more robust than L-BFGS. The results shown in Figure 7 suggest that the use of a smaller number of direction vectors m can be still efficient, i.e., more efficient algorithms, e.g., with adaptive m (or an adaptive m × n transformation matrix) can be designed. It seems both promising and feasible to extend the algorithm to constrained, noisy and/or multi-objective optimization, the domains, which are both hardly accessible for L-BFGS and keenly demanded by practitioners. As an important contribution to the success in this direction, it would be helpful to implement online adaptation of internal hyper-parameters as already implemented in the original CMA-ES [START_REF] Loshchilov | A Computationally Efficient Limited Memory CMA-ES for Large Scale Optimization[END_REF]. This would ensure an additional level of invariance and robustness on large scale black-box optimization problems.

Figure 1 :

 1 Figure 1: Number of memory slots of floating point variables required to run different optimization algorithms versus problem dimension n.

Figure 2 :

 2 Figure2: Timing results of LM-CMA and VD-CMA averaged over the whole run on the separable Ellipsoid compared to timing results of simple operations averaged over 100 seconds of experiments. The results for L-BFGS are not shown but for an optimized implementation would be comparable to one scalar-vector multiplication.

Figure 3 :

 3 Figure 3: The trajectories show the median of 11 runs of LM-CMA, L-BFGS with exact gradients provided at the cost of n + 1 evaluations per gradient, CL-BFGS with central differencing, active CMA-ES and VD-CMA on 512-Separable/Original (Left Column) and Rotated (Right Column) functions.

Figure 4 :

 4 Figure 4: The trajectories show the median of 11 runs of LM-CMA, L-BFGS with exact gradients, CL-BFGS with central differencing, active CMA-ES and VD-CMA on 512-Separable (Left Column) and Rotated (Right Column) functions.

Figure 5 :

 5 Figure 5: Median (out of 11 runs) number of function evaluations required to find f (x) = 10 -10 for LM-CMA, L-BFGS with exact gradients, CL-BFGS with central differencing, active CMA-ES and VD-CMA. Dotted lines depict extrapolated results.

Figure 6 :

 6 Figure6: The trajectories show the median of 11 runs of LM-CMA, L-BFGS (with exact gradients provided at the cost of n+1 evaluations per gradient) on 100,000-dimensional Rosenbrock and Ellipsoid functions.

Figure 7 :

 7 Figure 7: The trajectories show the median of 11 runs of LM-CMA in default settings, CL-BFGS in default settings and tuned LM-CMA (all three algorithms are with restarts) on the second nonsmooth variant of Nesterov-Chebyshev-Rosenbrock function in dimensions 128 and 2048.

Figure 8 :

 8 Figure 8: Sensitivity of LM-CMA to different settings of m.

Figure 9 :

 9 Figure 9: Eigenspectrums of C t = A t A t T for t denoting iteration of LM-CMA with m direction vectors (m = 5, m = 4 + ⌊3 ln 1024⌋ = 24, m = ⌊2 √ n⌋ = 64) on 1024dimensional problems. Darker (blue) lines correspond to later iterations. The number of function evaluations to reach f (x) = 10 -10 is given in the title of each sub-figure.

Evolutionary Computation

Evolutionary Computation

http://www.cs.ubc.ca/ ˜schmidtm/Software/minFunc.html

http://www.lri.fr/ ˜hansen/cmaes.m

http://sites.google.com/site/ecjlmcma/

Acknowledgments

I am grateful to Michèle Sebag and Marc Schoenauer for many valuable discussions and insights. I also would like to thank Oswin Krause, Youhei Akimoto and the anonymous reviewers whose interest and valuable comments helped to improve this work.