arXiv:1511.00221v1 [cs.NE] 1 Nov 2015

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

ilya.loshchilov@gmail.com
Laboratory of Intelligent Systems (LIS)
École Polytechnique Fédéral de Lausanne (EPFL), Lausanne, Switzerland

Abstract
The limited memory BFGS method (L-BFGS) of Liu and Nocedal (1989) 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.
In this paper, we demonstrate an alternative to L-BFGS, the limited memory Covariance Matrix Adaptation Evolution Strategy (LM-CMA) proposed by Loshchilov (2014).
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.
Keywords
LM-CMA, L-BFGS, CMA-ES, large scale optimization, black-box optimization.

1 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 (Ollivier et al., 2011). In continuous domain, f is defined as a mapping
Rn → 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

I. Loshchilov

The use of well recognized gradient-based approaches such as the BroydenFletcher-Goldfarb-Shanno (BFGS) algorithm (Shanno, 1970) 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 (Li et al.,
2007).
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 (Nocedal, 1980). The final algorithm called the limited
memory BFGS algorithm (L-BFGS or LM-BFGS) proposed by Liu and Nocedal (1989)
is still considered to be the state-of-the-art of large scale gradient-based optimization
(Becker and Fadili, 2012). 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 (Hansen et al., 2003). This algorithm constantly demonstrates good performance at various platforms for comparing continuous optimizers such as the Black-Box Optimization Benchmarking (BBOB) workshop
(Finck et al., 2010; Auger et al., 2010; Loshchilov et al., 2013) and the Special Session at
Congress on Evolutionary Computation (Garcı́a et al., 2009; Loshchilov, 2013a). The
CMA-ES was also extended to noisy (Hansen et al., 2009), expensive (Kern et al., 2006;
Loshchilov et al., 2012) and multi-objective optimization (Igel et al., 2007).
The principle advantage of CMA-ES, the learning of dependencies between n decision variables, also forms its main practical limitations such as O(n2 ) memory storage
required to run the algorithm and O(n2 ) computational time complexity per function
evaluation (Ros and Hansen, 2008). 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 (sepCMA-ES by Ros and Hansen (2008)), linear time Natural Evolution Strategy (R1-NES
by Sun et al. (2011)) and VD-CMA by Akimoto et al. (2014) 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.
2

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

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 Loshchilov (2014). Instead of storing the covariance matrix, the LMCMA 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.

2 The L-BFGS
An early version of the L-BFGS method, at that time called the SQN method, was proposed by Nocedal (1980). 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, Hk is obtained by applying m
BFGS updates to a sparse symmetric and positive definite matrix H0 provided by the
user (Liu and Nocedal, 1989).
Let us denote iterates by xk , sk = xk+1 − xk and yk = gk+1 − gk , where g denotes
gradient. The method uses the inverse BFGS formula in the form
Hk+1 = VTk Hk Vk + ρk sk sTk ,
1/yTk sk ,

(1)

ρk yk sTk

where ρk =
and Vk = I −
(Dennis Jr and Schnabel, 1996;
Liu and Nocedal, 1989).
The L-BFGS method works as follows (Liu and Nocedal, 1989):
Step 1. Choose x0 , m, 0 < β ′ < 1/2, β ′ < β < 1, and a symmetric and positive
definite starting matrix H0 . Set k = 0.
Step 2. Compute
dk = −Hk gk ,

(2)

xk+1 = xk + αk dk ,

(3)

where αk satisfies the Wolfe conditions (Wolfe, 1969):
f (xk + αk dk ) ≤ f (xk ) + β ′ αk gTk dk ,

(4)

g(xk + αk dk )T dk ≥ βgTk dk .

(5)

The novelty introduced by Liu and Nocedal (1989) w.r.t. the version given in
Nocedal (1980) 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

3

I. Loshchilov

ok
n
Step 3. Let m̂ = min(k, m−1). Update H0 m̂+1 times using the pairs yj , sj

j=k−m̂

as follows:

Hk+1 = (VTk · · · VTk−m̂ )H0 (Vk−m̂ · · · Vk )
+ρk−m̂ (VTk · · · VTk−m̂+1 )sk−m̂ sTk−m̂ (Vk−m̂+1 · · · Vk )
+ρk−m̂+1 (VTk · · · VTk−m̂+2 )sk−m̂+1 sTk−m̂+1 (Vk−m̂+2 · · · Vk )
..
.
+ρk sk sTk
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 LBFGS-B has the efficiency of the original algorithm, however at the cost of a significantly more complex implementation (Byrd et al., 1995). 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 (Becker and Fadili, 2012; Ngiam et al.,
2011) and often can be viewed as a method of choice (Andrew and Gao, 2007) 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.

3 Evolution Strategies for Large Scale Optimization
Historically, first Evolution Strategies (Rechenberg, 1973) 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 (Hansen and Ostermeier,
1996; Hansen et al., 2003). 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 Hansen et al.
(2015). More specifically, the analysis of theoretical foundations of Evolution Strategies
is provided by Wierstra et al. (2014); Ollivier et al. (2011); Akimoto and Ollivier (2013);
Glasmachers (2012); Auger and Hansen (2013); Hansen and Auger (2014); Arnold
(2014); Beyer (2014).
3.1 The CMA-ES
The Covariance Matrix Adaptation Evolution Strategy (Hansen and Ostermeier,
1996, 2001; Hansen et al., 2003) 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 mt of the mutation distribution (can be interpreted as an estimation of the optimum) is used to generate its k-th out of λ candidate solution xk ∈ Rn (line 5) by adding
a random Gaussian mutation defined by a (positive definite) covariance matrix

4

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

Algorithm 1 The (µ/µw , λ)-CMA-ES
1:

given n ∈ N+ , λ = 4 + ⌊3 ln n⌋, µ = ⌊λ/2⌋, wi =

ln(µ+ 12 )−ln i
1
j=1 (ln(µ+ 2 )−ln

Pµ

q

j)

µw +2
µw −1
n+µw +3 , dσ = 1 + cσ + 2 max(0,
n+1 − 1),
2 (µw −2+1/µw )
2 min(1,λ/6)
(n+1.3)2 +µw , cµ =
(n+2)2 +µw
t=0
n
2: initialize m
∈ R , σ t=0 > 0, pt=0
= 0, pt=0
= 0, Ct=0 = I, t ← 0
σ
c

µw =

3:

4:
5:
6:
7:
8:
9:
10:
11:
12:

Pµ 1

i=1

wi2

, cσ =

for i = 1 . . . µ,

cc =

4
n+4 ,

c1 =

repeat
for k ← 1, . . . , λ do
xk ← mt + σ t N (0, Ct )
fk ← fP
(xk )
µ
mt+1 ← i=1 wi xi:λ // the symbol i : λ denotes i-th best individual on f
p
t
√
− 1 t+1
pt+1
← (1 − cσ )ptσ + cσ (2 − cσ ) µw Ct 2 m σt−m
σ
hσ ← 11 pt+1 <√1−(1−c )2(t+1) (1.4+ 2 ) EkN(0,I )k
k σ k
σ
n+1
p
√ mt+1 −mt
t
c
µw σ t
pt+1
←
(1
−
c
)p
+
h
(2
−
c
)
c
σ
c
c
c
c
Pµ
(xi:λ −mt )T
xi:λ −mt
Cµ ← i=1 wi σt ×
σt
T

+cµ Cµ
Ct+1 ← (1 − c1 − cµ )Ct + c1 pt+1
pt+1
|{z}
| c {zc }

rank−one update rank−µ update

13:
14:
15:

t+1

t

k
kpt+1
σ
exp( dcσσ ( EkN(0,I
)k

σ
←σ
− 1))
t←t+1
until stopping criterion is met

Ct ∈ Rn×n and a mutation step-size σ t as follows:



2
xtk ← N mt , σ t Ct ← mt + σ t N 0, Ct

(6)

These λ solutions then should be evaluated with an objective function f (line 6).
The old mean of the mutation distribution is stored in mt and a new mean mt+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)
(Hansen and Ostermeier, 1996) and is controlled by evolution path pt+1
σ . Successful
mt+1 −mt
(line 8) are tracked in the space of sampling, i.e., in the isotropic
mutation steps
σt
coordinate system defined by principal components of the covariance matrix Ct . To
update the evolution path pt+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 pt+1
is longer than the expected
σ
length of the evolution path under random selection E kN (0, I )k,√and decreases other1
1
wise (line 13). The expectation of kN (0, I )k is approximated by n(1 − 4n
+ 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
(Hansen and Ostermeier, 2001) and rank-µ update (Hansen et al., 2003). The rank-one
t+1
t
update computes evolution path pt+1
of successful moves of the mean m σt−m of the
c
Evolutionary Computation

5

I. Loshchilov

mutation distribution in the given coordinate system (line 10), in a similar way as for
the evolution path pt+1
of the step-size. To stall the update of pt+1
when σ increases
σ
c
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
t+1 T
of rank-one update
weights of importance (line 12): c1 for covariance matrix pt+1
c pc
and cµ for Cµ of rank-µ update (Hansen et al., 2003) such that c1 + cµ ≤ 1. It was
also proposed to take into account unsuccessful mutations in the ”active” rank-µ update
(Hansen and Ros, 2010; Jastrebski and Arnold, 2006).
In CMA-ES, the factorization of the covariance C into AAT = C is needed to
sample the multivariate normal distribution (line 5). The eigendecomposition with
O(n3 ) 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(n2 ).
3.2 Large Scale Variants
The original O(n2 ) 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 Hansen (2008). The algorithm does not learn dependencies but the scaling of
variables by restraining the covariance matrix update to the diagonal elements:

t
ct+1
jj = (1 − ccov )cjj +

1
µcov

pt+1
c

2
j


+ cccov 1 −

1
µccov

X
µ
i=1

wi ctjj zi:λ t+1

2
j

, j = 1, . . . , n (7)

where, for j = 1, . . . , n the cjj are the diagonal elements of Ct and the zi:λ t+1
q

xi:λ t+1 j /(σ t ctjj ).



j

=

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
Sun et al. (2011), adapts the search distribution according to the natural gradient with
a particular low rank parametrization of the covariance matrix,
C = σ 2 (I + uuT ),

(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
Ollivier et al. (2011).
The IGO framework was used to derive a similar to R1-NES algorithm called VDCMA (Akimoto et al., 2014) with the sampling distribution parametrized by a

6

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

Gaussian model with the covariance matrix restricted as follows:
C = D(I + uuT )D,

(9)

where D is a diagonal matrix of dimension n and u is a vector in Rn . 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 (Akimoto et al., 2014). It has O(n) time and
O(µn) space complexity but i) typically demonstrates a better performance than sepCMA-ES and R1-NES and ii) can solve a larger class of functions (Akimoto et al., 2014).
A version of CMA-ES with a limited memory storage also called limited memory
CMA-ES (L-CMA-ES) was proposed by Knight and Lunacek (2007). 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 Brand (2006) 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 Rn for n > 10, √
000. However, the computational complexity of
O(m2 n) practically (for m in order of n as suggested by Knight and Lunacek (2007))
leads to the same limitations of O(n2 ) time complexity as in the original CMA-ES.
The (µ/µw , λ)-Cholesky-CMA-ES proposed by Suttorp et al. (2009) is of special interest in this paper because the LM-CMA is based on this algorithm. The CholeskyCMA represents a version of CMA-ES with rank-one update where instead of performT
ing the factorization of the covariance matrix Ct into At At = Ct , the Cholesky factor
−1
At and its inverse At are iteratively updated. From Theorem 1 (Suttorp et al., 2009)
it follows that if Ct is updated as
T

Ct+1 = αCt + βvt vt ,

(10)

n

t

t t

+

where v ∈ R is given in the decomposition form v = A z , and α, β ∈ R , then
for z 6= 0 a Cholesky factor of the matrix Ct+1 can be computed by
!
r
√
√
α
β t 2
T
t
t+1
(11)
A
= αA +
1 + kz k − 1 [At zt ]zt ,
α
kzt k2
√
for zt = 0 we have At+1 = αAt . From the Theorem 2 (Suttorp et al., 2009) it
t
follows that if A−1 is the inverse of At , then the inverse of At+1 can be computed by


1
1
1
t
t+1
1 − q
 zt [zt T A−1 t ],
(12)
A−1
= √ A−1 − √
2
t
α
αkz k
1 + β kzt k2
α

−1 t+1

−1 t

for zt = 0.
for zt 6= 0 and by A
= √1α A
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(n2 ). 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

7

I. Loshchilov

Algorithm 2 The (µ/µw , λ)-Cholesky-CMA-ES
1:

given n ∈ N+ , λ = 4 + ⌊3 ln n⌋, µ = ⌊λ/2⌋, wi =
µw =
c1 =

2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

Pµ 1

i=1

wi2

, cσ =

√

µw
√
√
,
n+ µw

ln(µ+1)−ln(i)
P
;i
µ ln(µ+1)− µ
j=1 ln(j)

dσ = 1 + cσ + 2 max(0,

2
√ 2
(n+ 2)
t=0

q

µw −1
n+1

= 1 . . . µ,

− 1), cc =

4
n+4 ,

initialize m
∈ Rn , σ t=0 > 0, pσt=−1 = 0, pct=−1 = 0, At=0 = I, At=0
inv = I, t ← 0
repeat
for k ← 1, . . . , λ do
zk ← N (0, I)
xk ← mt + σ t Azk
fk ← fP
(xk )
µ
mt+1 ←
Pµ i=1 wi xi:λ
zw ← i=1 wi zi:λ p
√
ptσ ← (1 − cσ )pt−1
+ cσ (2 − cσ ) µw zw
σ
p
√
ptc ← (1 − cc )pt−1
+ cc (2 − cc ) µw At zw
c
t
t
t
v ← Ainv pc
q

√
√
c1
1
t k2 − 1 pt vt T
kv
1
+
At+1 ← 1 − c1 At + kv1−c
2
tk
1−c1
!c
At+1
inv ←
t+1

√ 1
At
1−c1 inv
t

−

1
√
1−c1 kvt k2

1−

q

1
c1
1+ 1−c
kvt k2

T

vt [vt Atinv ],

1

kptσ k
exp( dcσσ ( EkN(0,I
)k

σ
←σ
− 1))
t←t+1
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 (Suttorp et al., 2009). While it has the same
computational and memory complexity, the lack of rank-µ update may deteriorate its
performance on problems where it is essential.

4 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 Loshchilov (2014) 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.
8

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

Algorithm 3 Az(): Cholesky factor - vector update
1:
2:
3:
4:
5:
6:

|i|

m×n
given z ∈ Rn , m ∈ Z+ , j ∈ Zm
, V ∈ Rm×n , b ∈ Rm , a ∈ [0, 1]
+ , i ∈ Z+ , P ∈ R
initialize x ← z
for t ← 1, . . . , |i| do
k ← bjit V(jit ,:) · z
x ← ax + kP(jit ,:)
return x

Algorithm 4 Ainvz(): inverse Cholesky factor - vector update
1:
2:
3:
4:
5:
6:

|i|

m
given z ∈ Rn , m ∈ Z+ , j ∈ Zm
+ , i ∈ Z+ , d ∈ R , c ∈ [0, 1]
initialize x ← z
for t ← 1, . . . , |i| do
k ← djit V(jit ,:) · x
x ← cx − kV(jit ,:)
return x

4.1 Reconstruction of Cholesky Factors
q
√
√
1
By setting a = 1 − c1 , bt = kv1−c
1+
t k2

c1
t 2
1−c1 kv k


− 1 and considering the evolu-

tion path ptc (a change of optimum estimate m smoothed over iterations, see line 12 of
t
Algorithm 7) together with vt = A−1 ptc , one can rewrite Equation (11) as
T

At+1 = aAt + bt ptc vt ,

(13)

The product of a random vector z and the Cholesky factor At thus can be directly
computed. At iteration t = 0, A0 = I and A0 z = z, the new updated Cholesky factor
T
T
T
A1 = aI + b0 p0c v0 . At iteration t = 1, A1 z = (aI + b0 p0c v0 )z = az + b0 p0c (v0 z) and
T
T
A2 = a(aI + b0 p0c v0 ) + b1 p1c v1 . Thus, a very simple iterative procedure which scales
as O(mn) can be used to sample candidate solutions in IRn according to the Cholesky
factor At reconstructed from m pairs of vectors ptc and vt .
The pseudo-code of the procedure to reconstruct x = At z from m direction vectors
is given in Algorithm 3. At each iteration of reconstruction of x = At z (lines 3 - 4), x
is updated as a sum of a-weighted version of itself and bt -weighted evolution path ptc
(accessed from a matrix P ∈ Rm×n as P(it ,:) ) scaled by the dot product of vt and x. As
can be seen, Algorithms 3 and 4 use jit indexation instead of t. This is a convenient
way to have references to matrices P and V which store pc and v vectors, respectively.
In the next subsections, we will show how to efficiently manipulate these vectors.
−1
A very similar approach can be used to reconstruct x = At √
z, for the sake of
reproducibility the pseudo-code!is given in Algorithm 4 for c = 1/ 1 − c1 and dt =
√

1
1−c1 kvt k2

×

1−

1
q
c1
1+ 1−c
kvt k2

. The computational complexity of both procedures

1

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 pt+1 (see line 15 in Algorithm 7), all inverse vectors
vt for t = ℓ, . . . , m should be iteratively recomputed (Krause, 2014). This procedure
corresponds to line 16 in Algorithm 7.
Evolutionary Computation

9

I. Loshchilov

Algorithm 5 UpdateSet(): direction vectors update
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:

m
m
m×n
given m ∈ R+ , j ∈ Zm
, pc ∈ Rn , T ∈ Z+
+ , l ∈ Z+ , t ∈ Z+ , N ∈ Z+ , P ∈ R
t ← ⌊t/T ⌋
if t ≤ m then
jt ← t
else


imin ← 1 + argmini lji+1 − lji − Ni , |1 ≤ i ≤ (m − 1)
if lji
− lji −1 − Ni ≥ 0 then
min
min
imin ← 1
jtmp ← jimin
for i ← imin , . . . , m − 1 do
ji ← ji+1
jm ← jtmp
jcur ← jmin(t+1,m)
ljcur ← tT
P(jcur ,:) ← pc
return: jcur , j, l

4.2 Direction Vectors Selection and Storage
The choice to store only m ≪ n direction vectors pc 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 Loshchilov (2014) 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 cc and from µ ≪ n selected vectors.
The selection procedure is outlined in Algorithm 5 which outputs an array of
pointers j such that j1 points out to a row in matrices P and V with the oldest saved
vectors pc and v which will be taken into account during the reconstruction procedure.
The higher the index i of ji the more recent the corresponding direction vector is. The
index jcur 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
(j
,:)
(j
,:)
is the following. Find a pair of consecutively saved vectors (P imin −1 , P imin ) with
the distance between them (in terms of indexes of iterations, stored in l) closest to a
target distance Ni (line 6). If this distance is smaller than Ni then the index jimin will
(j
,:)
be swapped with last index of j (lines 9-12) and the corresponding vector P imin will
be replaced by the new vector pc (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 Ni and with the overall time horizon
Pm−1
for all vectors of at most i
Ni iterations. The procedure can be called periodically
every T ∈ Z+ iterations of the algorithm. The values of Ni are to be defined, e.g., as a
function of problem dimension n and direction vector index i. Here, however, we set
Ni to n for all i, i.e., the target distance equals to the problem dimension.
10

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

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

4.3 Sampling from a Family of Search Representations
At iteration t, a new k-th solution can be generated as
xk ← mt + σ t Az(zk , i),

(14)

n

where zk ∈ R is a vector drawn from some distribution and transformed by a
Cholesky factor by calling Az(zk , 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(j1 ,:) refers to the oldest vector, P(jit ,:) refers to the it -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 (Loshchilov, 2014) 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(n3/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

2
half-normal distribution N 0, mσ (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 xk appears to be sampled from an isotropic normal
distribution, the computation of v inverses is performed using all m vectors.
4.4 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 (Yao and Liu, 1997; Hansen, 2008). Moreover, the
Adaptive Encoding procedure proposed by Hansen (2008) can be coupled with any
sampling distribution as in Loshchilov et al. (2011), 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

11

I. Loshchilov

In this paper, inspired by Loshchilov et al. (2011), 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 Gaus√
−cβ 2
sian, for any positive
of the mass lies
√ real number β ≤√ n, all but at most 3 exp
within the annulus n − 1 − β ≤ r ≤ n − 1 + β, where, c is a fixed positive constant
(Hopcroft and Kannan, 2015). 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.
4.5 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 Loshchilov (2014). 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 (Ait Elhara et al., 2013). To estimate
the success of the current population we combine fitness function values from the previous and current population into a mixed set
fmix ← f t−1 ∪ f t

(15)

Then, all individuals in the mixed set are ranked to define two sets rt−1 and rt (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
zP SR ←

Pλ

i=1

rt−1 (i) − rt (i)
− z ∗,
λ2

(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 ← (1 − cσ )s + cσ zP 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.
4.6 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 zk sampled from
the Rademacher distribution and a Cholesky factor At reconstructed from m∗ out of
12

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

Algorithm 7 The (µ/µw , λ)-LM-CMA
1:

given n ∈ N+ , λ = 4 + ⌊3 ln n⌋, µ = ⌊λ/2⌋, wi =
µw =

Pµ 1

i=1

1
10 ln(n+1) ,

2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

wi2

20:
21:
22:

, cσ = 0.3, z = 0.25, m = 4 + ⌊3 ln n⌋, Nsteps = n, cc =

= 1 . . . µ,
0.5
√
,
n

c1 =

dσ = 1, T = ⌊log(n)⌋
initialize mt=0 ∈ Rn , σ t=0 > 0, pt=0
= 0, s ← 0, t ← 0
c
repeat
for k ← 1, . . . , λ do
if k (mod 2) = 1 then
zk ← Rademacher()
i ← SelectSubset(k)
xk ← mt + σ t Az(zk , i)
else
xk ← mt − (xk−1 − mt )
t
fk ← fP
(xk )
mt+1 ← µi=1 wi xi:λ p
√
pt+1
← (1 − cc )ptc + cc (2 − cc ) µw (mt+1 − mt )/σ t
c
if t (mod T ) = 0 then
U pdateSet(pt+1
c )
U pdateInverses()
rt , rt−1 ←PRanks of f t and f t−1 in f t ∪ f t−1
λ

18:
19:

∗

ln(µ+1)−ln(i)
P
;i
µ ln(µ+1)− µ
j=1 ln(j)

rt−1 (i)−rt (i)

− z∗
zP SR ← i=1 λ2
s ← (1 − cσ )s + cσ zP SR
σ t+1 ← σ t exp(s/dσ )
t←t+1
until stopping criterion is met

m vectors (line 4-11) as described in Sections 4.1-4.4. We introduce the mirrored sampling (Brockhoff et al., 2010) to generate the actual xk only every second time and thus
decrease the computation cost per function evaluation by a factor of two by evaluating
mt + σ t Az(zk ) and then its mirrored version mt − (xk−1 − mt ). The latter approach
sometimes also improves the convergence rate.
The best µ out of λ solutions are selected to compute the new mean mt+1 of the
mutation distribution in line 12. The new evolution path pt+1
is updated (line 13) from
c
√
the change of the mean vector µw (mt+1 −mt )/σ t and represents an estimate of descent
direction. One vector among m vectors is selected and replaced by the new pt+1
in
c
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.

5 Experimental Validation
The performance of the LM-CMA is investigated comparatively to the L-BFGS
(Liu and Nocedal, 1989), the active CMA-ES by Hansen and Ros (2010) and the VDCMA by Akimoto et al. (2014). 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 Akimoto et al. (2014).
Evolutionary Computation

13

I. Loshchilov

Table 1: Test functions, initialization intervals and initial standard deviation (when
applied). R is an orthogonal n × n matrix with each column vector qi being a
uniformly distributed unit vector implementing an angle-preserving transformation
(Ros and Hansen, 2008)
.
0
Name
Sphere

Ellipsoid
Rosenbrock
Discus
Cigar
Different Powers

Function
Pn
2
fSphere (x)=
i=1 xi
Pn
6 i−1 2
n−1 x
fElli (x)=
i
i=1 10

Pn−1
100.(x2i − xi+1 )2 + (xi − 1)2
fRosen (x)=
i=1
P
2
fDiscus (x)= 106 x21 + n
i=2 xi
P
2
fCigar (x)= x21 + 106 n
x
i=2 i
Pn
2+4(i−1)/(n−1)
fDif f P ow (x)=
i=1 |xi |

Target f (x)
10−10

Init
[−5, 5]n

σ
3

10−10

[−5, 5]n

3

−10

10

[−5, 5]n

3

10−10

[−5, 5]n

3

10−10

[−5, 5]n

3

10−10

[−5, 5]n

3

Rotated Ellipsoid

fRotElli (x)=fElli (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

fRotDif f P ow (x)=fDif f P ow (Rx)

10−10

[−5, 5]n

3

We use the L-BFGS implemented in MinFunc library by Schmidt (2005) in its default parameter settings 1 , the active CMA-ES (aCMA) without restarts in its default
parametrization of CMA-ES MATLAB code version 3.61 2 . 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 online 3 . 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 (Finck et al., 2010). 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.
1 http://www.cs.ubc.ca/

˜schmidtm/Software/minFunc.html
˜hansen/cmaes.m
3 http://sites.google.com/site/ecjlmcma/

2 http://www.lri.fr/

14

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

10

10

1 solution
VDCMA
L−BFGS
LM−CMA
n solutions
aCMA

9

Number of memory slots

10

8

10

7

10

6

10

5

10

4

10

3

10

2

10

3

10

4

10

5

10

6

10

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

5.1 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 VDCMA 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 × 108 slots (respectively, 8 × 108 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,
Evolutionary Computation

15

I. Loshchilov

−1

10

Scalar−Vector
Rademacher Vector
Ellipsoid Evaluation
LM−CMA
Gaussian Vector
VDCMA
Matrix−Vector

−2

CPU cost in sec.

10

−3

10

−4

10

−5

10

−6

10

−7

10

−8

10

2

10

3

10

4

10

Dimension
Figure 2: 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.

the algorithm automatically reduces m to a feasible mf .
5.2 Time Complexity
The average amount of CPU time internally required by an algorithm per evaluation
of some objective function f ∈ Rn (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 seconds, 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 orig16

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

inal version by Loshchilov (2014) 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)).
5.3 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 (Ollivier et al., 2011). Moreover, this invariance property provides robustness to noise as far as this noise does not impact a comparison of solutions
of interest (Auger and Hansen, 2013).
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.
5.4 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
(Hansen et al., 2011). 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 particular 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 (Loshchilov et al., 2013). 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 Ct=0 = R−1 (R−1 )T , and the initial
search point(s) are transformed accordingly Hansen (2006). 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,
Evolutionary Computation

17

I. Loshchilov

Separable Ellipsoid

10

10
LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

8

10

6

10

LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

8

10

6

10

4

Objective Function

Objective Function

Rotated Ellipsoid

10

10

10

2

10

0

10

−2

10

−4

10

−6

4

10

2

10

0

10

−2

10

−4

10

−6

10

10

−8

−8

10

10

−10

−10

10

2

10

3

4

10

5

10

6

10

7

10

10

8

10

10

2

10

3

4

10

Rosenbrock

7

10

8

10

10

Rotated Rosenbrock
LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

8

6

10

4

LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

8

10

6

10

4

10

Objective Function

10

Objective Function

6

10

Function Evaluations

10

2

10

0

10

−2

10

−4

10

−6

2

10

0

10

−2

10

−4

10

−6

10

10

−8

−8

10

10

−10

10

5

10

Function Evaluations

−10

2

10

3

10

4

10

5

10

6

10

7

10

Function Evaluations

8

10

9

10

10

2

10

3

10

4

10

5

10

6

10

7

10

8

10

9

10

Function Evaluations

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.

differences in performance on axis-aligned and rotated ill-conditioned functions were
marginal.
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
3-4 show 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 fRotElli , fRotDiscus and fRotDif f P ow functions due to the restricted form
of the adapted covariance matrix of Equation (9). Both fRosen , fCigar and their rotated
18

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

Separable Discus

8

10
LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

6

10

4

4

10

2

10

0

10

−2

10

−4

10

−6

2

10

0

10

−2

10

−4

10

−6

10

10

−8

−8

10

10

−10

10

LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

6

10

Objective Function

10

Objective Function

Rotated Discus

8

10

−10

0

10

1

10

2

10

3

4

10

10

5

6

10

10

7

10

8

10

0

10

10

1

10

2

10

Function Evaluations

Separable Cigar
10

6

10

7

8

10

10

LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

10

Objective Function

Objective Function

5

10

10

5

10

0

10

−5

5

10

0

10

−5

10

10

−10

−10

10

2

3

10

4

10

5

10

10

6

10

10

2

3

10

4

10

5

10

Function Evaluations

6

10

10

Function Evaluations

Separable Sum of Different Powers

Rotated Sum of Different Powers
LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

6

10

4

10

LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

6

10

4

10

2

10

Objective Function

Objective Function

4

10

Rotated Cigar
LM−CMA
L−BFGS
aCMA
VDCMA
CL−BFGS

10

0

10

−2

10

−4

10

−6

2

10

0

10

−2

10

−4

10

−6

10

10

−8

−8

10

10

−10

10

3

10

Function Evaluations

−10

0

10

1

10

2

10

3

10

4

10

5

10

Function Evaluations

6

10

7

10

8

10

10

0

10

1

10

2

10

3

10

4

10

5

10

6

10

7

10

8

10

Function Evaluations

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 512Separable (Left Column) and Rotated (Right Column) functions.
versions can be solved efficiently since they have a Hessian matrix whose inverse can
be well approximated by Equation (9) (Akimoto et al., 2014).
An important observation from Figures 3-4 is that even the exact gradient inforEvolutionary Computation

19

I. Loshchilov

Sphere

Rosenbrock

6

Function Evaluations / Dimension

Function Evaluations / Dimension

10
2

10

1

10

0

10

LM−CMA
VDCMA
aCMA
L−BFGS
CL−BFGS
2

10

5

10

LM−CMA
VDCMA
aCMA
L−BFGS
CL−BFGS

4

10

3

10

2

3

10

10

4

10

2

3

10

10

LM−CMA
aCMA
L−BFGS
CL−BFGS

Function Evaluations / Dimension

Function Evaluations / Dimension

Cigar

3

10

4

10

3

10

LM−CMA
VDCMA
aCMA
L−BFGS
CL−BFGS

2

10

1

2

10

3

10

10

4

10

2

3

10

10

Dimension
Rotated Discus

4

Rotated Sum of Different Powers

5

10

LM−CMA
aCMA
L−BFGS
CL−BFGS

Function Evaluations / Dimension

Function Evaluations / Dimension

4

10

Dimension

10

3

10

2

10

1

10

10

Dimension

Rotated Ellipsoid

5

4

10

Dimension

4

10

LM−CMA
aCMA
L−BFGS
CL−BFGS

3

10

2

10

1

2

10

3

10

Dimension

4

10

10

2

10

3

10

4

10

Dimension

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.

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 fElli
and fRotElli . The L-BFGS with the central difference method (CL-BFGS) experiences
an early triggering of stopping criteria on fRosen and fDif 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.
20

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

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 fElli .
5.5 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 3-4, 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 3-4), 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 3-4). 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 fRotDif f P ow , where the loss is increasing with problem
dimension.
Quite surprisingly, the LM-CMA outperforms VD-CMA and aCMA on fSphere .
This performance is close to the one obtained for (1+1) Evolution Strategy with optimal step-size. Bad performance on fSphere 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 ∗ 6= 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 zP 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 (Schumer and Steiglitz, 1968;
Rechenberg, 1973). We consider the PSR to be its population-based version.
The performance of LM-CMA on fElli 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 fElli since the number of parameters of the full covariance matrix
to learn is (n2 + n)/2 (Hansen and Ostermeier, 2001). 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
Evolutionary Computation

21

I. Loshchilov

Rosenbrock

10

Separable Ellipsoid

11

10

10
LM−CMA
L−BFGS

LM−CMA
L−BFGS

9

10

10

Objective Function

Objective Function

10

8

10

7

10

9

10

8

10
6

10

5

10

7

3

10

4

10

5

10

Function Evaluations

6

10

7

10

10

3

10

4

10

5

10

6

10

7

10

Function Evaluations

Figure 6: 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.
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.
The LM-CMA is better or is as good as VD-CMA on fRosen and fCigar 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 fRosen suggests that the
problem is more difficult (e.g., more difficult than fElli ) 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.
5.6 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 Overton (2015) demonstrated the difficulties encountered by BFGS on some nonsmooth functions. We selected one of the test functions
from Overton (2015) called the second nonsmooth variant of Nesterov-ChebyshevRosenbrock function defined as follows:
22

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

128−D

4

10

LM−CMA default
CL−BFGS default
LM−CMA tuned

2

LM−CMA default
CL−BFGS default
LM−CMA tuned

2

10

Objective Function

Objective Function

10

0

10

0

10

−2

10

−2

10

−4

10

−4

10

−6

10

−6

10

−8

−8

10

10

−10

10

2048−D

4

10

−10

0

2000

4000

6000

8000

10000

10

Function Evaluations / Dimension

0

2000

4000

6000

8000

10000

Function Evaluations / Dimension

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.

n−1

N̂ (x) =

X
1
|xi+1 − 2 |xi | + 1|
|x1 − 1| +
4
i=1

(18)

This function is nonsmooth (though locally Lipschitz) as well as nonconvex, it has
2n−1 Clarke stationary points (Overton, 2015). Overton (2015) 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 Abramson and Audet (2006))
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 (Overton,
2015). 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 c1 = 15/(10 ln(n + 1)), iii) an extremely
small learning rate for step-size adaptation cσ = 0.3/n2 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.
Evolutionary Computation

23

I. Loshchilov

Rotated Ellipsoid

5

10

Function Evaluations / Dimension

Function Evaluations / Dimension

m=2
m=5
m=10
m=3+4log(n)
m=2sqrt(n)

4

10

Rotated Discus

5

10

m=2
m=5
m=10
m=3+4log(n)
m=2sqrt(n)

4

10

3

2

3

10

10

10

4

10

2

3

10

10

Dimension
Cigar

4

4

10

Dimension
Rotated Sum of Different Powers

m=2
m=5
m=10
m=3+4log(n)
m=2sqrt(n)

Function Evaluations / Dimension

Function Evaluations / Dimension

10

3

10

2

10

4

10

m=2
m=5
m=10
m=3+4log(n)
m=2sqrt(n)

3

2

10

3

10

Dimension

4

10

10

2

10

3

10

4

10

Dimension

Figure 8: Sensitivity of LM-CMA to different settings of m.
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.
5.7 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 Nsteps , the learning rate cc for the
evolution path and learning rate c1 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 Nsteps = n on
problems where a global coordinate system is expected to be constant, e.g., on a class
24

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

3

3

Rosenbrock, m=5, 4866x10 evals.

4

4

10

3

1

2

10

1

0

−1

10
−3
10

0

10

10

Ellipsoid, m=5, 14311x10 evals.

10
−3
10

0

10

6

5

Ellipsoid, m=24, 12943x10 evals.

6

5

2

10

4

10

Eigenvalue

Eigenvalue

10

3

10

2

10

1

10
−3
10

−2

10

−1

10
−3
10

0

10

10

0

−2

10

−1

10
−3
10

0

10

10

Eigenvalue Index / Dimension

3

7

10

6

7

6

6

5

3

10

2

10

Eigenvalue

Eigenvalue

10

4

10

3

10

2

10

1

1

0

1

−2

10

−1

0

10

10

Eigenvalue Index / Dimension

0

−2

10

−1

10
−3
10

0

10

10

Eigenvalue Index / Dimension

3

6

10

Discus, m=24, 4974x10 evals.

6

10

4

2

10

4

10

Eigenvalue

Eigenvalue

3

10

3

10

2

10

1

−2

10

−1

0

10

10

Eigenvalue Index / Dimension

0

−2

10

−1

10
−3
10

0

10

10

Eigenvalue Index / Dimension

3

Sum of Diff. Pow., m=5, 2728x10 evals.

Sum of Diff. Pow., m=24, 2124x10 evals.

5

6

10

5

5

10

10

4

Eigenvalue

3

10

2

10

4

10

Eigenvalue

4

3

10

2

10

1

−2

10

−1

10

0

10

10
−3
10

2

10

0

Eigenvalue Index / Dimension

3

10

1

10

0

10

10

1

10

0

10

Sum of Diff. Pow., m=64, 1972x10 evals.

6

10

−1

10

3

10

10

−2

10

Eigenvalue Index / Dimension

3

10

10
−3
10

2

10

0

10
−3
10

3

10

1

10

0

10
−3
10

10

10

1

10

Discus, m=64, 4186x10 evals.

5

10

4

10

0

10

10

5

10

−1

10

3

10

5

−2

10

Eigenvalue Index / Dimension

3

Discus, m=5, 6712x10 evals.

6

10

0

10
−3
10

3

10
10

10

10
−3
10

4

10

2

10

10

Cigar, m=64, 180x10 evals.

10

5

10

4

0

10

10

10

5

10

−1

10

3

Cigar, m=24, 196x10 evals.

10

10

−2

10

Eigenvalue Index / Dimension

3

Cigar, m=5, 473x10 evals.

7

2

10

0

Eigenvalue Index / Dimension

3

10

1

10

0

10

10

1

10

Ellipsoid, m=64, 12357x10 evals.

10

4

3

0

10

10

5

10

−1

10

3

10

4

−2

10

Eigenvalue Index / Dimension

10

10

Eigenvalue

−1

10

3

10

Eigenvalue

0

−2

10

Eigenvalue Index / Dimension

3

6

10

0

−2

10

Eigenvalue Index / Dimension

2

10

1

10

10
−3
10

Eigenvalue

10

Eigenvalue

Eigenvalue

Eigenvalue

2

10

Rosenbrock, m=64, 6180x10 evals.

3

10

10

Eigenvalue

4

10

3

10

6

3

Rosenbrock, m=24, 5832x10 evals.

10

0

−2

10

−1

10

Eigenvalue Index / Dimension

0

10

10
−3
10

−2

10

−1

10

0

10

Eigenvalue Index / Dimension

T

Figure 9: Eigenspectrums of Ct = At At 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

25

I. Loshchilov

of problems described by the general ellipsoid model (Beyer, 2014). The learning rate
0.5
which is
for the evolution path is related to both m and n, here, we set it to cc = √
√ n
roughly inversely proportional to the (if affordable) suggested m = ⌊2 n⌋. We found
that the chosen cc is still valid for the default m = 4 + ⌊3 ln n⌋. We do not have a
1
good interpretation for the learning rate c1 = 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 Loshchilov et al. (2014).
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 > 106 . 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 fDiscus 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 fElli , fDiscus and fDif f P ow . The latter
is not the case for fCigar , 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 Ct constructed as
t tT
A A from the Cholesky factor At . 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 fCigar , 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 fElli , fDiscus and fDif f P ow , that in turn leads to a better performance.
However, this is not always the case as can be seen for fRosen that again demonstrates
that optimal hyper-parameter settings are problem-dependent.

6 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 (Knight and Lunacek, 2007) but with a modified sampling procedure to obtain linear time complexity as proposed for the Adaptive Coordinate Descent by Loshchilov
(2013b).
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 sam26

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

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 LMCMA 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 (Loshchilov, 2014). This would ensure an additional
level of invariance and robustness on large scale black-box optimization problems.

Acknowledgments
I am grateful to Michè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.

References
Abramson, M. A. and Audet, C. (2006). Convergence of mesh adaptive direct search to secondorder stationary points. SIAM Journal on Optimization, 17(2):606–619.
Ait Elhara, O., Auger, A., and Hansen, N. (2013). A median success rule for non-elitist evolution
strategies: Study of feasibility. In Genetic and Evolutionary Computation Conference, pages 415–
422. ACM.
Akimoto, Y., Auger, A., and Hansen, N. (2014). Comparison-based natural gradient optimization
in high dimension. In Genetic and Evolutionary Computation Conference, pages 373–380. ACM.
Akimoto, Y. and Ollivier, Y. (2013). Objective improvement in information-geometric optimization. In Twelfth workshop on Foundations of Genetic Algorithms XII, pages 1–10. ACM.
Andrew, G. and Gao, J. (2007). Scalable training of l 1-regularized log-linear models. In 24th
International Conference on Machine Learning, pages 33–40. ACM.
Arnold, D. V. (2014). On the behaviour of the (1,λ)-ES for conically constrained linear problems.
Evolutionary Computation, 22(3):503–523.
Auger, A., Finck, S., Hansen, N., and Ros, R. (2010). BBOB 2010: Comparison Tables of All
Algorithms on All Noiseless Functions. Technical Report RR-7215, INRIA.
Auger, A. and Hansen, N. (2013). Linear Convergence of Comparison-based Step-size Adaptive
Randomized Search via Stability of Markov Chains. arXiv preprint arXiv:1310.7697.
Becker, S. and Fadili, J. (2012). A quasi-newton proximal splitting method. In Advances in Neural
Information Processing Systems, pages 2618–2626.
Evolutionary Computation

27

I. Loshchilov

Beyer, H.-G. (2014). Convergence analysis of evolutionary algorithms that are based on the
paradigm of information geometry. Evolutionary Computation, 22(4):679–709.
Brand, M. (2006). Fast low-rank modifications of the thin singular value decomposition. Linear
algebra and its applications, 415(1):20–30.
Brockhoff, D., Auger, A., Hansen, N., Arnold, D. V., and Hohm, T. (2010). Mirrored sampling
and sequential selection for evolution strategies. In Parallel Problem Solving from Nature–PPSN,
pages 11–21. Springer.
Byrd, R. H., Lu, P., Nocedal, J., and Zhu, C. (1995). A limited memory algorithm for bound
constrained optimization. SIAM Journal on Scientific Computing, 16(5):1190–1208.
Dennis Jr, J. E. and Schnabel, R. B. (1996). Numerical methods for unconstrained optimization and
nonlinear equations, volume 16. SIAM.
Finck, S., Hansen, N., Ros, R., and Auger, A. (2010). Real-parameter black-box optimization
benchmarking 2010: Experimental setup. Technical Report 2009/21, Research Center PPE.
Garcı́a, S., Molina, D., Lozano, M., and Herrera, F. (2009). 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. Journal of Heuristics, 15:617–644.
Glasmachers, T. (2012). Convergence of the IGO-Flow of isotropic gaussian distributions on
convex quadratic problems. In Parallel Problem Solving from Nature–PPSN, pages 1–10. Springer
Berlin Heidelberg.
Hansen, N. (2006). The cma evolution strategy: a comparing review. In Towards a new evolutionary
computation, pages 75–102. Springer.
Hansen, N. (2008). Adaptive encoding: How to render search coordinate system invariant. In
Parallel Problem Solving from Nature–PPSN, pages 205–214. Springer.
Hansen, N., Arnold, D. V., and Auger, A. (2015). Evolution strategies. In Springer Handbook of
Computational Intelligence, pages 871–898. Springer.
Hansen, N., Atamna, A., and Auger, A. (2014). How to Assess Step-Size Adaptation Mechanisms
in Randomised Search. In Parallel Problem Solving from Nature–PPSN, pages 60–69. Springer.
Hansen, N. and Auger, A. (2014). Principled design of continuous stochastic search: From theory
to practice. In Theory and Principled Methods for the Design of Metaheuristics, pages 145–180.
Springer Berlin Heidelberg.
Hansen, N., Müller, S., and Koumoutsakos, P. (2003). Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES). Evolutionary
Computation, 11(1):1–18.
Hansen, N., Niederberger, A. S., Guzzella, L., and Koumoutsakos, P. (2009). A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion. Evolutionary Computation, IEEE Transactions on, 13(1):180–197.
Hansen, N. and Ostermeier, A. (1996). Adapting Arbitrary Normal Mutation Distributions in
Evolution Strategies: The Covariance Matrix Adaptation. In International Conference on Evolutionary Computation, pages 312–317.
Hansen, N. and Ostermeier, A. (2001). Completely Derandomized Self-Adaptation in Evolution
Strategies. Evolutionary Computation, 9(2):159–195.
Hansen, N. and Ros, R. (2010). Benchmarking a weighted negative covariance matrix update
on the BBOB-2010 noiseless testbed. In Genetic and Evolutionary Computation Conference, pages
1673–1680. ACM.

28

Evolutionary Computation

LM-CMA for Large Scale Black-box Optimization

Hansen, N., Ros, R., Mauny, N., Schoenauer, M., and Auger, A. (2011). Impacts of invariance in
search: When CMA-ES and PSO face ill-conditioned and non-separable problems. Applied Soft
Computing, 11(8):5755–5769.
Hopcroft, J. and Kannan, R. (2015). Foundations of Data Science. to appear.
Igel, C., Hansen, N., and Roth, S. (2007). Covariance matrix adaptation for multi-objective optimization. Evolutionary computation, 15(1):1–28.
Jastrebski, G. A. and Arnold, D. V. (2006). Improving Evolution Strategies through Active Covariance Matrix Adaptation. In IEEE Congress on Evolutionary Computation, pages 2814–2821.
Kern, S., Hansen, N., and Koumoutsakos, P. (2006). Local Meta-Models for Optimization Using
Evolution Strategies. In PPSN, pages 939–948. LNCS 4193, Springer Verlag.
Knight, J. N. and Lunacek, M. (2007). Reducing the space-time complexity of the CMA-ES. In
Genetic and Evolutionary Computation Conference, pages 658–665. ACM.
Krause, O. (2014). Personal communication.
Li, Z.-C., Chien, C.-S., and Huang, H.-T. (2007). Effective condition number for finite difference
method. Journal of computational and applied mathematics, 198(1):208–235.
Liu, D. C. and Nocedal, J. (1989). On the limited memory BFGS method for large scale optimization. Mathematical programming, 45(1-3):503–528.
Loshchilov, I. (2013a). CMA-ES with restarts for solving CEC 2013 benchmark problems. In
Evolutionary Computation (CEC), 2013 IEEE Congress on, pages 369–376. IEEE.
Loshchilov, I. (2013b). Surrogate-Assisted Evolutionary Algorithms. PhD thesis, Université Paris
Sud-Paris XI.
Loshchilov, I. (2014). A Computationally Efficient Limited Memory CMA-ES for Large Scale
Optimization. In Genetic and Evolutionary Computation Conference, pages 397–404. ACM.
Loshchilov, I., Schoenauer, M., and Sebag, M. (2011). Adaptive Coordinate Descent. In Genetic
and Evolutionary Computation Conference, pages 885–892. ACM.
Loshchilov, I., Schoenauer, M., and Sebag, M. (2012). Self-adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy. In Genetic and Evolutionary Computation Conference,
pages 321–328. ACM.
Loshchilov, I., Schoenauer, M., and Sebag, M. (2013). Bi-population CMA-ES agorithms with
surrogate models and line searches. In Genetic and Evolutionary Computation Conference, pages
1177–1184. ACM.
Loshchilov, I., Schoenauer, M., Sebag, M., and Hansen, N. (2014). Maximum Likelihood-based
Online Adaptation of Hyper-parameters in CMA-ES. In Parallel Problem Solving from Nature–
PPSN, pages 70–79. Springer International Publishing.
Ngiam, J., Coates, A., Lahiri, A., Prochnow, B., Le, Q. V., and Ng, A. Y. (2011). On optimization
methods for deep learning. In International Conference on Machine Learning, pages 265–272.
Nocedal, J. (1980). Updating quasi-newton matrices with limited storage. Math. of computation,
35(151):773–782.
Ollivier, Y., Arnold, L., Auger, A., and Hansen, N. (2011). Information-geometric optimization
algorithms: A unifying picture via invariance principles. arXiv preprint arXiv:1106.3708.
Overton, M. L. (2015). Nonsmooth, nonconvex optimization algorithms and examples. Slides of
a talk given at Challenges in Optimization for Data Science, July 1–2, 2015, Universite Pierre
et Marie Curie – Paris 6, Paris, France.
Evolutionary Computation

29

I. Loshchilov

Rechenberg, I. (1973). Evolutionsstrategie: optimierung technischer systeme nach prinzipien der biologischen evolution. Frommann-Holzboog.
Ros, R. and Hansen, N. (2008). A simple modification in CMA-ES achieving and space complexity. In Parallel Problem Solving from Nature–PPSN, pages 296–305. Springer.
Schmidt, M. (2005).
Supérieure, Paris.

Minfunc.

Technical report, Laboratoire d’Informatique de l’Ecole Normale

Schumer, M. and Steiglitz, K. (1968). Adaptive step size random search. Automatic Control, IEEE
Transactions on, 13(3):270–276.
Shanno, D. F. (1970). Conditioning of Quasi-Newton Methods for Function Minimization. Math.
of Computation, 24(111):647–656.
Sun, Y., Gomez, F., Schaul, T., and Schmidhuber, J. (2011). A linear time natural evolution strategy
for non-separable functions. arXiv preprint arXiv:1106.1998.
Suttorp, T., Hansen, N., and Igel, C. (2009). Efficient covariance matrix update for variable metric
evolution strategies. Machine Learning, 75(2):167–197.
Wierstra, D., Schaul, T., Glasmachers, T., Sun, Y., Peters, J., and Schmidhuber, J. (2014). Natural
evolution strategies. The Journal of Machine Learning Research, 15(1):949–980.
Wolfe, P. (1969). Convergence conditions for ascent methods. SIAM review, 11(2):226–235.
Yao, X. and Liu, Y. (1997). Fast evolution strategies. In Evolutionary Programming VI, pages 149–
161. Springer.

30

Evolutionary Computation

