The CMA Evolution Strategy: A Tutorial

Nikolaus Hansen
June 28, 2011

Contents
Nomenclature 3

0 Preliminaries 4
0.1 Eigendecomposition of a Positive Definite Matrix . . . . . . . . . . . . . . . 5
0.2 The Multivariate Normal Distribution . . . . . . . . . . . . . . . . . . . . . 6
0.3 Randomized Black Box Optimization . . . . . . . . . . . . . . . . . . . . . 7
0.4 Hessian and Covariance Matrices . . . . . . . . . . . . . . . . . . . . . . . . 8

1 Basic Equation: Sampling 8

2 Selection and Recombination: Moving the Mean 9

3 Adapting the Covariance Matrix 10

3.1 Estimating the Covariance Matrix From Scratch . . . . . . . . . . . . . . . . 10
3.2 Rank-µ-Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3.3 Rank-One-Update . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3.1 A Different Viewpoint . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.3.2 Cumulation: Utilizing the Evolution Path . . . . . . . . . . . . . . . 16
3.4 Combining Rank-µ-Update and Cumulation . . . . . . . . . . . . . . . . . . 17

4 Step-Size Control 17

5 Discussion 21

A Algorithm Summary: The CMA-ES 25

B Implementational Concerns 27
B.1 Multivariate normal distribution . . . . . . . . . . . . . . . . . . . . . . . . 28
B.2 Strategy internal numerical effort . . . . . . . . . . . . . . . . . . . . . . . . 28
B.3 Termination criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
B.4 Flat fitness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
B.5 Boundaries and Constraints. . . . . . . . . . . . . . . . . . . . . . . . . . . 29

C MATLAB Source Code 31

1
D Reformulation of Learning Parameter ccov 33

2
Nomenclature
We adopt the usual vector notation, where bold letters, v, are column vectors, capital bold
letters, A, are matrices, and a transpose is denoted by v T . A list of used abbreviations and
symbols is given in alphabetical order.

Abbreviations
CMA Covariance Matrix Adaptation
EMNA Estimation of Multivariate Normal Algorithm
ES Evolution Strategy
(µ/µ{I,W} , λ)-ES, Evolution Strategy with µ parents, with recombination of all µ parents,
either Intermediate or Weighted, and λ offspring.
RHS Right Hand Side.

Greek symbols
λ ≥ 2, population size, sample size, number of offspring, see (5).
µ ≤ λ parent number, number of selected search points in the population, see (6).
2 −1
Pµ

µeff = i=1 wi = 1/kwk2 , the variance effective selection mass, see (8).
σ (g) ∈ R+ , step-size.

Latin symbols
B ∈ Rn , an orthogonal matrix. Columns of B are eigenvectors of C with unit length and
correspond to the diagonal elements of D.
C (g) ∈ Rn×n , covariance matrix at generation g.
cii , diagonal elements of C.
cc ≤ 1, learning rate for cumulation for the rank-one update of the covariance matrix, see
(22) and (42), and Table 1.
c1 ≤ 1 − cµ , learning rate for the rank-one update of the covariance matrix update, see (26),
(27), and (43), and Table 1.
cµ ≤ 1 − c1 , learning rate for the rank-µ update of the covariance matrix update, see (14),
(27), and (43), and Table 1.
cσ < 1, learning rate for the cumulation for the step-size control, see (28) and (40), and
Table 1.
D ∈ Rn , a diagonal matrix. The diagonal elements of D are square roots of eigenvalues of
C and correspond to the respective columns of B.
di > 0, diagonal elements of diagonal matrix D, d2i are eigenvalues of C.
dσ ≈ 1, damping parameter for step-size update, see (29), (34), and (41).
E Expectation value

3
 f : Rn → R, x 7→ f (x), objective function (fitness function) to be minimized.
Pn
fsphere : Rn → R, x 7→ kxk2 = xT x = i=1 x2i .
g ∈ N0 , generation counter, iteration number.
I ∈ Rn×n , Identity matrix, unity matrix.
m(g) ∈ Rn , mean value of the search distribution at generation g.
n ∈ N, search space dimension, see f .
N (0, I), multivariate normal distribution with zero mean and unity covariance matrix. A
vector distributed according to N (0, I) has independent, (0, 1)-normally distributed
components.
N (m, C) ∼ m + N (0, C), multivariate normal distribution with mean m ∈ Rn and
covariance matrix C ∈ Rn×n . The matrix C is symmetric and positive definite.
p ∈ Rn , evolution path, a sequence of successive (normalized) steps, the strategy takes over
a number of generations.
wi , where i = 1, . . . , µ, recombination weights, see (6).
(g+1)
xk ∈ Rn , k-th offspring/individual from generation g + 1. We also refer to x(g+1) , as
search point, or object parameters/variables, commonly used synonyms are candidate
solution, or design variables.
(g+1) (g+1) (g+1)
xi:λ , i-th best individual out of x1 , . . . , xλ , see (5). The index i : λ denotes the
(g+1) (g+1) (g+1)
index of the i-th ranked individual and f (x1:λ ) ≤ f (x2:λ ) ≤ · · · ≤ f (xλ:λ ),
where f is the objective function to be minimized.
(g+1) (g+1)
yk = (xk − m(g) )/σ (g) corresponding to xk = m + σy k .

0 Preliminaries
This tutorial introduces the CMA Evolution Strategy (ES), where CMA stands for Covariance
Matrix Adaptation. The CMA-ES is a stochastic method for real-parameter (continuous do-
main) optimization of non-linear, non-convex functions (see also Section 0.3 below).1 We try
to motivate and derive the algorithm from intuitive concepts and from requirements of non-
linear, non-convex search in continuous domain. In order to refer to the described algorithm,
also cite [9]. For finding a concise algorithm description go directly to Appendix A. The
respective Matlab source code is given in Appendix C.
Before we start to introduce the algorithm in Sect. 1, a few required fundamentals are
summed up.
1 While CMA variants for multi-objective optimization and elitistic variants have been proposed, this tutorial is

solely dedicated to single objective optimization and to non-elitistic truncation selection, also referred to as comma-
selection.

4
0.1 Eigendecomposition of a Positive Definite Matrix
A symmetric, positive definite matrix, C ∈ Rn×n , is characterized in that for all x ∈ Rn \{0}
holds xT Cx > 0. The matrix C has an orthonormal basis of eigenvectors, B = [b1 , . . . , bn ],
with corresponding eigenvalues, d21 , . . . , d2n > 0.
That means for each bi holds
Cbi = d2i bi . (1)
The important message from (1) is that eigenvectors are not rotated by C. This feature
uniquely distinguishes eigenvectors.
 Because we assume the orthogonal eigenvectors to be
T 1 if i = j
of unit length, bi bj = δij = , and B T B = I (obviously this means
0 otherwise
B −1 = B T , and it follows BB T = I). An basis of eigenvectors is practical, because
v ∈ Rn we can find coefficients αi , such that v = i αi bi , and then we have
P
for anyP
Cv = i d2i αi bi .
The eigendecomposition of C obeys
C = BD 2 B T , (2)
where
B is an orthogonal matrix, B T B = BB T = I. Columns of B form an orthonormal basis
of eigenvectors.
D 2 = DD = diag(d1 , . . . , dn )2 = diag(d21 , . . . , d2n ) is a diagonal matrix with eigenvalues
of C as diagonal elements.
D = diag(d1 , . . . , dn ) is a diagonal matrix with square roots of eigenvalues of C as diagonal
elements.

The matrix decomposition (2) is unique, apart from signs of columns of B and permutations
of columns in B and D 2 respectively, given all eigenvalues are different.2
Given the eigendecomposition (2), the inverse C −1 can be computed via

−1
C −1 = BD 2 B T
−1
= BT D −2 B −1
= B D −2BT 
1 1
= B diag , . . . , BT .
d21 d2n
From (2) we naturally define the square root of C as
1
C 2 = BDB T (3)
and therefore
1
C− 2 = BD −1 B
T

1 1
= B diag ,..., BT
d1 dn
2 Given m eigenvalues are equal, any orthonormal basis of their m-dimensional subspace can be used as column

vectors. For m > 1 there are infinitely many such bases.

5



N 0, σ 2 I N 0, D 2 N (0, C)

Figure 1: Ellipsoids depicting one-σ lines of equal density of six different normal distribu-
tions, where σ ∈ R+ , D is a diagonal matrix, and C is a positive definite full covariance
matrix. Thin lines depict possible objective function contour lines

0.2 The Multivariate Normal Distribution

A multivariate normal distribution, N (m, C), has a unimodal, “bell-shaped” density, where
the top of the bell (the modal value) corresponds to the distribution mean, m. The distribution
N (m, C) is uniquely determined by its mean m ∈ Rn and its symmetric and positive definite
covariance matrix C ∈ Rn×n . Covariance (positive definite) matrices have an appealing
geometrical interpretation: they can be uniquely identified with the (hyper-)ellipsoid {x ∈
Rn | xT C −1 x = 1}, as shown in Fig. 1. The ellipsoid is a surface of equal density of the
distribution. The principal axes of the ellipsoid correspond to the eigenvectors of C, the
squared axes lengths correspond to the eigenvalues. The eigendecomposition is denoted by
2
C = B (D) B T (see Sect. 0.1). If D = σI, where σ ∈ R+ and I denotes the identity
matrix, C = σ 2 I and the ellipsoid is isotropic (Fig. 1, left). If B = I, then C = D 2 is a
diagonal matrix and the ellipsoid is axis parallel oriented (middle). In the coordinate system
given by the columns of B, the distribution N (0, C) is always uncorrelated.
The normal distribution N (m, C) can be written in different ways.

N (m, C) ∼ m + N (0, C)
1
∼ m + C 2 N (0, I)
∼ m + BD B T N (0, I)
| {z }
∼ N(0,I)

∼ m + B DN (0, I) , (4)
| {z }
∼ N(0,D 2 )

1
where “∼” denotes equality in distribution, and C 2 = BDB T . The last row can be well
interpreted, from right to left
N (0, I) produces an spherical (isotropic) distribution as in Fig. 1, left.

6
Initialize distribution parameters θ (0)
For generation g = 0, 1, 2, . . .

Sample λ independent points from distribution P x|θ (g) → x1 , . . . , xλ
Evaluate the sample x1 , . . . , xλ on f
Update parameters θ (g+1) = Fθ (θ (g) , (x1 , f (x1 )), . . . , (xλ , f (xλ )))
break, if termination criterion met

Figure 2: Randomized black box search. f : Rn → R is the objective function

D scales the spherical

distribution within the coordinate axes as in Fig. 1, middle. DN (0, I) ∼
N 0, D 2 has n independent components. The matrix D can be interpreted as (indi-
vidual) step-size matrix and its diagonal entries are the standard deviations of the com-
ponents.
B defines a new orientation for the ellipsoid, where the new principal axes of the ellipsoid
2
correspond to the columns of B. Note that B has n 2−n degrees of freedom.

Equation (4) is useful to compute N (m, C) distributed vectors, because N (0, I) is a vector
of independent (0, 1)-normally distributed numbers that can easily be realized on a computer.

0.3 Randomized Black Box Optimization

We consider the black box search scenario, where we want to minimize an objective function
(or cost function or fitness function)

f : Rn → R
x 7→ f (x) .

The objective is to find one or more search points (candidate solutions), x ∈ Rn , with a func-
tion value, f (x), as small as possible. We do not state the objective of searching for a global
optimum, as this is often neither feasible nor relevant in practice. Black box optimization
refers to the situation, where function values of evaluated search points are the only accessible
information on f .3 The search points to be evaluated can be freely chosen. We define the
search costs as the number of executed function evaluations, in other words the amount of
information we needed to aquire from f 4 . Any performance measure must consider the search
costs together with the achieved objective function value.5
A randomized black box search algorithm is outlined in Fig. 2. In the CMA Evolution
Strategy the search distribution, P , is a multivariate normal distribution. Given all variances
and covariances, the normal distribution has the largest entropy of all distributions in Rn .
3 Knowledge about the underlying optimization problem might well enter the composition of f and the chosen

problem encoding.
4 Also f is sometimes denoted as cost function, but it should not to be confused with the search costs.
5 A performance measure can be obtained from a number of trials as, for example, the mean number of function

evaluations to reach a given function value, or the median best function value obtained after a given number of
function evaluations.

7
Furthermore, coordinate directions are not distinguished in any way. Both makes the normal
distribution a particularly attractive candidate for randomized search.
Randomized search algorithms are regarded to be robust in a rugged search landscape,
which can comprise discontinuities, (sharp) ridges, or local optima. The covariance matrix
adaptation (CMA) in particular is designed to tackle, additionally, ill-conditioned and non-
separable6 problems.

0.4 Hessian and Covariance Matrices

We consider the convex-quadratic objective function fH : x 7→ 12 xT Hx, where the Hessian
matrix H ∈ Rn×n is a positive definite matrix. Given a search distribution N (m, C), there is
−1
a close relation between H and C: Setting C P= H on fH is equivalent to optimizing the
1 T 1 2
isotropic function fsphere (x) = 2 x x = 2 i xi (where H = I) with C = I.7 That is, on
convex-quadratic objective functions, setting the covariance matrix of the search distribution
to the inverse Hessian matrix is equivalent to rescaling the ellipsoid function into a spherical
one. Consequently, we assume that the optimal covariance matrix equals to the inverse Hessian
matrix, up to a constant factor.8 Furthermore, choosing a covariance matrix or choosing a
respective affine linear transformation of the search space (i.e. of x) is equivalent [8], because
for any full rank n × n-matrix A we find a positive definite Hessian such that 12 (Ax)T Ax =
1 T T 1 T
2 x A Ax = 2 x Hx.
The final objective of covariance matrix adaptation is to closely approximate the contour
lines of the objective function f . On convex-quadratic functions this amounts to approximating
the inverse Hessian matrix, similar to a quasi-Newton method.
In Fig. 1 the solid-line distribution in the right figure follows the objective function con-
tours most suitably, and it is easy to foresee that it will aid to approach the optimum the most.
The condition number of a positive definite matrix A is defined via the Euclidean norm:
def
cond(A) = kAk × kA−1 k, where kAk = supkxk=1 kAxk. For a positive definite (Hessian
or covariance) matrix A holds kAk = λmax and cond(A) = λλmax min
≥ 1, where λmax and
λmin are the largest and smallest eigenvalue of A.

1 Basic Equation: Sampling

In the CMA Evolution Strategy, a population of new search points (individuals, offspring) is
generated by sampling a multivariate normal distribution.9 The basic equation for sampling
the search points, for generation number g = 0, 1, 2, . . . , reads10
 
(g+1)
xk ∼ m(g) + σ (g) N 0, C (g) for k = 1, . . . , λ (5)

6 An n-dimensional separable problem can be solved by solving n 1-dimensional problems separately, which is a

far easier task.

7 Also the initial mean value m has to be transformed accordingly.
8 Even though there is good intuition and strong empirical evidence for this statement, a rigorous proof is missing.
9 Recall that, given all (co-)variances, the normal distribution has the largest entropy of all distributions in Rn .
10 Framed equations belong to the final algorithm of a CMA Evolution Strategy.

8
where
∼ denotes the same distribution on the left and right side.
N(0, C (g) ) is a multivariate normal distribution with zero mean and covariance (g)

matrix C ,
(g) (g) (g) (g) (g) 2 (g)
see Sect. 0.2. It holds m + σ N(0, C ) ∼ N m , (σ ) C .
(g+1)
xk ∈ Rn , k-th offspring (individual, search point) from generation g + 1.
m(g) ∈ Rn , mean value of the search distribution at generation g.
σ (g) ∈ R+ , “overall” standard deviation, step-size, at generation g.
2
C (g) ∈ Rn×n , covariance matrix at generation g. Up to the scalar factor σ (g) , C (g) is the
covariance matrix of the search distribution.
λ ≥ 2, population size, sample size, number of offspring.

To define the complete iteration step, the remaining question is, how to calculate m(g+1) ,
C (g+1) , and σ (g+1) for the next generation g + 1. The next three sections will answer
these questions, respectively. An algorithm summary with all parameter settings and MAT-
LAB source code is given in Appendix A and C, respectively.

2 Selection and Recombination: Moving the Mean

The new mean m(g+1) of the search distribution is a weighted average of µ selected points
(g+1) (g+1)
from the sample x1 , . . . , xλ :

µ
(g+1)
X
(g+1)
m = wi xi:λ (6)
i=1

µ
X
wi = 1, w1 ≥ w2 ≥ · · · ≥ wµ > 0 (7)
i=1

where
µ ≤ λ is the parent population size, i.e. the number of selected points.
wi=1...µ ∈ R+ , positive weight coefficients for recombination. For wi=1...µ = 1/µ, Equation
(6) calculates the mean value of µ selected points.
(g+1) (g+1) (g+1)
xi:λ , i-th best individual out of x1 , . . . , xλ from (5). The index i : λ denotes the
(g+1) (g+1) (g+1)
index of the i-th ranked individual and f (x1:λ ) ≤ f (x2:λ ) ≤ · · · ≤ f (xλ:λ ),
where f is the objective function to be minimized.

9
Equation (6) implements truncation selection by choosing µ < λ out of λ offspring points.
Assigning different weights wi must also be interpreted as a selection mechanism. Equation
(6) implements weighted intermediate recombination by taking µ > 1 individuals into account
for a weighted average.
The measure
2 µ
!−1
kwk21

kwk1 1 X
2
µeff = = = = wi (8)
kwk2 kwk22 kwk22 i=1

will be repeatedly used in the following and can be paraphrased as variance effective selection
mass. From the definition of wi in (7) we derive 1 ≤ µeff ≤ µ, and µeff = µ for equal
recombination weights, i.e. wi = 1/µ for all i = 1 . . . µ. Usually, µeff ≈ λ/4 indicates a
reasonable setting of wi . A typical setting could be wi ∝ µ − i + 1, and µ ≈ λ/2.

3 Adapting the Covariance Matrix

In this section, the update of the covariance matrix, C, is derived. We will start out estimating
the covariance matrix from a single population of one generation (Sect. 3.1). For small pop-
ulations this estimation is unreliable and an adaptation procedure has to be invented (rank-µ-
update, Sect. 3.2). In the limit case only a single point can be used to update (adapt) the covari-
ance matrix at each generation (rank-one-update, Sect. 3.3). The adaptation can be enhanced
by exploiting dependencies between successive steps applying cumulation (Sect. 3.3.2). Fi-
nally we combine the rank-µ and rank-one updating methods (Sect. 3.4).

3.1 Estimating the Covariance Matrix From Scratch

For the moment we assume that the population contains enough information to reliably es-
timate a covariance matrix from the population.11 For the sake of convenience we assume
σ (g) = 1 (see (5)) in this section. For σ (g) 6= 1 the formulae hold except for a constant factor.
We can (re-)estimate the original covariance matrix C (g) using the sampled population
(g+1) (g+1)
from (5), x1 . . . xλ , via the empirical covariance matrix
  T
λ λ λ
1 X
x(g+1) 1 X (g+1)   (g+1) 1 X (g+1) 
C (g+1)
emp = i − x xi − x . (9)
λ − 1 i=1 λ j=1 j λ j=1 j

(g+1)
The empirical covariance matrix C emp is an unbiased estimator of C (g) : assuming the
(g+1)
xi , i = 1 . . . λ, to be random variables (rather than a realized sample), we have that
 (g+1)  (g) 
E C emp  C = C (g) . Consider now a slightly different approach to get an estimator for
C (g) .
λ T
(g+1) 1 X  (g+1) 
(g+1)
Cλ = xi − m(g) xi − m(g) (10)
λ i=1
11 To re-estimate the covariance matrix, C, from a N (0, I) distributed sample such that cond(C) < 10 a sample

size λ ≥ 4n is needed, as can be observed in numerical experiments.

10
 (g+1)
Also the matrix C λ is an unbiased estimator of C (g) . The remarkable difference between
(g+1)
(9) and (10) is the reference mean value. For C emp it is the mean of the actually realized
(g+1)
sample. For C λ it is the true mean value, m(g) , of the sampled distribution (see (5)).
(g+1) (g+1) (g+1)
Therefore, the estimators C emp and C λ can be interpreted differently: while C emp
(g+1)
estimates the distribution variance within the sampled points, C λ estimates variances of
(g+1)
sampled steps, xi − m(g) .
1
A minor difference between (9) and (10) is the different normalizations λ−1 versus λ1 ,
necessary to get an unbiased estimator in both cases. In (9) one degree of freedom is
already taken by the inner summand. In order to get a maximum likelihood estimator in
both cases λ1 must be used.
Equation (10) re-estimates the original covariance matrix. To “estimate” a “better” co-
variance matrix (10) is modified and the same, weighted selection mechanism as in (6) is used
[11].
µ   T
(g+1) (g+1)
X
C (g+1)
µ = wi xi:λ − m (g)
x i:λ − m (g)
(11)
i=1
(g+1) (g+1)
The matrix Cµ an estimator for the distribution of selected steps, just as C λ
is is an
(g+1)
estimator of the original distribution of steps before selection. Sampling from C µ tends to
reproduce selected, i.e. successful steps, giving a justification for what a “better” covariance
matrix means.
Following [9], we compare (11) with the Estimation of Multivariate Normal Algorithm
EMNAglobal [18, 19]. The covariance matrix in EMNAglobal reads, similar to (9),
µ
(g+1) 1 X  (g+1) 
(g+1)
T
C EMNAglobal = xi:λ − m(g+1) xi:λ − m(g+1) , (12)
µ i=1
Pµ (g+1)
where m(g+1) = 1
µ i=1 xi:λ . Similarly, applying the so-called Cross-Entropy method
µ (g+1)
to continuous domain optimization [21] yields the covariance matrix µ−1 C EMNAglobal ,
i.e. the unbiased empirical covariance matrix of the µ best points. In both cases the subtle,
but most important difference to (11) is, again, the choice of the reference mean value.12
Equation (12) estimates the variance within the selected population while (11) estimates
selected steps. Equation (12) reveals always smaller variances than (11), because its ref-
erence mean value is the minimizer for the variances. Moreover, in most conceivable
selection situations (12) decreases the variances compared to C (g) .
Figure 3 demonstrates the estimation results on a linear objective function for λ = 150,
µ = 50, and wi = 1/µ. Equation (11) geometrically increases the expected variance
in direction of the gradient (where the selection takes place, here the diagonal), given
ordinary settings for parent number µ and recombination weights w1 , . . . , wµ . Equation
(12) always decreases the variance in gradient direction geometrically fast! Therefore, (12)
is highly susceptible to premature convergence, in particular with small parent populations,
where the population cannot be expected to bracket the optimum at any time. However,
for large values of µ in large populations with large initial variances, the impact of the
different reference mean value can become marginal.
12 Taking
Pµ 1 Pµ
a weighted sum, i=1 wi . . . , instead of the mean, µ i=1 . . . , is an appealing, but less important,
difference.

11
 (g+1)
Cµ

(g+1)
C EMNAglobal

sampling estimation new distribution

P2
Figure 3: Estimation of the covariance matrix on flinear (x) = − i=1 xi to be minimized.
Contour lines (dotted) indicate that the strategy should move toward the upper right corner.
(g+1)
Above: estimation of C µ according to (11), where wi = 1/µ. Below: estimation of
(g+1)
C EMNAglobal according to (12). Left: sample of λ = 150 N (0, I) distributed points. Middle:
the µ = 50 selected points (dots) determining the entries for the estimation equation (solid
straight lines). Right: search distribution of the next generation (solid ellipsoids). Given wi =
(g+1)
1/µ, estimation via C µ increases the expected variance in gradient direction for all µ <
(g+1)
λ/2, while estimation via C EMNAglobal decreases this variance for any µ < λ geometrically
fast

(g+1)
In order to ensure with (5), (6), and (11), that C µ is a reliable estimator, the variance
effective selection mass µeff (cf. (8)) must be large enough: getting condition numbers (cf.
(g) Pn
Sect. 0.4) smaller than ten for C µ on fsphere (x) = i=1 x2i , requires µeff ≈ 10n. The next
step is to circumvent this restriction on µeff .

3.2 Rank-µ-Update
To achieve fast search (opposite to more robust or more global search), e.g. competitive per-
formance on fsphere , the population size λ must be small. Because µeff ≈ λ/4 also µeff must
be small and we may assume, e.g., µeff ≤ 1 + ln n. Then, it is not possible to get a reliable
estimator for a good covariance matrix from (11). As a remedy, information from previous
generations is used additionally. For example, after a sufficient number of generations, the

12
mean of the estimated covariance matrices from all generations,
g
(g+1) 1 X 1
C = C (i+1) (13)
g + 1 i=0 σ (i) 2 µ

(g)
becomes a reliable estimator for the selected steps. To make C µ from different generations
comparable, the different σ (i) are incorporated. (Assuming σ (i) = 1, (13) resembles the
covariance matrix from the Estimation of Multivariate Normal Algorithm EMNAi [19].)
In (13), all generation steps have the same weight. To assign recent generations a higher
weight, exponential smoothing is introduced. Choosing C (0) = I to be the unity matrix and a
learning rate 0 < cµ ≤ 1, then C (g+1) reads
1
C (g+1) = (1 − cµ )C (g) + cµ 2 C (g+1)
µ
σ (g)
µ
(g+1) (g+1) T
X
= (1 − cµ )C (g) + cµ wi y i:λ y i:λ (14)
i=1
µ  !
(g) 1/2 (g+1) (g+1) T 1/2
X
= C I + cµ wi z i:λ z i:λ −I C (g)
i=1

where
cµ ≤ 1 learning rate for updating the covariance matrix. For cµ = 1, no prior information is
1 (g+1)
retained and C (g+1) = σ(g) 2 Cµ . For cµ = 0, no learning takes place and C (g+1) =
C (0) . Here, cµ ≈ min(1, µeff /n2 ) is a reasonably choice.
(g+1) (g+1)
y i:λ = (xi:λ − m(g) )/σ (g) .
(g+1) −1/2 (g+1)
z i:λ = C (g) y i:λ is the mutation vector expressed in the unique coordinate system
where the sampling is isotropic and the respective coordinate system transformation
does not rotate the original principal axes of the distribution.

This covariance matrix update is called rank-µ-update [13], because the sum of outer products
in (14) is of rank min(µ, n) (with probability one). Note that this sum can even consist of a
single term, if µ = 1.
The number 1/cµ is the backward time horizon that contains roughly 63% of the overall
weight.
Because (14) expands to the weighted sum
g
X 1
C (g+1) = (1 − cµ )g+1 C (0) + cµ (1 − cµ )g−i 2 C (i+1)
µ , (15)
i=0 σ (i)

the backward time horizon, ∆g, where about 63% of the overall weight is summed up, is
defined by
g
X 1
cµ (1 − cµ )g−i ≈ 0.63 ≈ 1 − . (16)
i=g+1−∆g
e

13
 Resolving the sum yields
1
(1 − cµ )∆g ≈
, (17)
e
and resolving for ∆g, using the Taylor approximation for ln, yields
1
∆g ≈ . (18)
cµ

That is, approximately 37% of the information in C (g+1) is older than 1/cµ generations,
and, according to (17), the original weight is reduced by a factor of 0.37 after approxi-
mately 1/cµ generations.13

The choice of cµ is crucial. Small values lead to slow learning, too large values lead to a
failure, because the covariance matrix degenerates. Fortunately, a good setting seems to be
largely independent of the function to be optimized.14 A first order approximation for a good
choice is cµ ≈ µeff /n2 . Therefore, the characteristic time horizon for (14) is roughly n2 /µeff .
Even for the learning rate cµ = 1, adapting the covariance matrix cannot be accomplished
within one generation. The effect of the original sample distribution does not vanish until a
sufficient number of generations. Assuming fixed search costs (number of function evalua-
tions), a small population size λ allows a larger number of generations and therefore usually
leads to a faster adaptation of the covariance matrix.

3.3 Rank-One-Update
In Section 3.1 we estimated the complete covariance matrix from scratch, using all selected
steps from a single generation. We now take precisely the opposite viewpoint. We will re-
peatedly update the covariance matrix in the generation sequence using a single selected step
only. First, this perspective will give a justification of the adaptation rule (14). Second, we
will introduce the so-called evolution path that is finally used for a rank-one update of the
covariance matrix.

3.3.1 A Different Viewpoint

We consider a specific method to produce n-dimensional normal distributions with zero mean.
Let the vectors y 1 , . . . , y g0 ∈ Rn , g0 ≥ n, span Rn and let N (0, 1) denote independent (0, 1)-
normally distributed random numbers, then
g0
!
X
T
N (0, 1) y 1 + · · · + N (0, 1) y g0 ∼ N 0, yi yi (19)
i=1
Pg0
is a normally distributed random vector with zero mean and covariance matrix i=1 yi yT
i .
The random vector (19) is generated by adding “line-distributions” N (0, 1) y i . The singu-
lar distribution N (0, 1) y i ∼ N(0, y i y T
i ) generates the vector y i with maximum likelihood
considering all normal distributions with zero mean.
13 This can be shown more easily, because (1 − c )g = exp ln(1 − c )g = exp(g ln(1 − c )) ≈ exp(−gc )
µ µ µ µ
for small cµ , and for g ≈ 1/cµ we get immediately (1 − cµ )g ≈ exp(−1).
14 We use the sphere model f 2
P
sphere (x) = i xi to empirically find a good setting for the parameter cµ , dependent
on n and µeff . The found setting was applicable to any non-noisy objective function we tried so far.

14




N 0, C (0) N 0, C (1) N 0, C (2)

Figure 4: Change of the distribution according to the covariance matrix update (20). Left:
vectors e1 and e2 , and C (0) = I = e1 eT T
1 + e2 e2 . Middle: vectors 0.91 e1 , 0.91 e2 , and
0.41 y 1 (the
coefficients deduce from c1 = 0.17), and C (1) = (1 − c1 ) I + c1 y 1 y T 1 , where
y 1 = −0.59
−2.2 . The distribution ellipsoid is elongated into the direction of y 1 , and therefore
0.97


increases the likelihood of y 1 . Right: C (2) = (1 − c1 ) C (1) + c1 y 2 y T
2 , where y 2 = 1.5 .

The line distribution that generates a vector y with the maximum likelihood must “live” on
a line that includes y, and therefore the distribution must obey N(0, 1)σy ∼ N(0, σ 2 yy T ).
Any other line distribution with zero mean cannot generate y at all. Choosing σ reduces to
choosing the maximum likelihood of kyk for the one-dimensional gaussian N(0, σ 2 kyk2 ),
which is σ = 1.
The covariance matrix yy T has rank one, its only eigenvectors are R\0 × y with
eigenvalue kyk2 . Using equation (19), any normal distribution can be realized if y i are
chosen appropriately. For example, (19) resembles (4) with m = 0, using the orthogonal
eigenvectors y i = dii bi , for i = 1, . . . , n, where bi are the columns of B. In general, the
vectors y i need not to be eigenvectors of the covariance matrix (and usually are not).
Considering (19) and a slight simplification of (14), we try to gain insight into the adapta-
tion rule for the covariance matrix. Let the sum in (14) consist of a single summand only (e.g.
(g+1)
x1:λ −m(g)
µ = 1), and let y g+1 = σ (g)
. Then, the rank-one update for the covariance matrix
reads
C (g+1) = (1 − c1 )C (g) + c1 y g+1 y g+1 T (20)
The right summand is of rank one and adds the maximum likelihood term for y g+1 into the
covariance matrix C (g) . Therefore the probability to generate y g+1 in the next generation
increases.
An example of the first two iteration steps of (20) is shown in Figure 4. The distribution
N(0, C (1) ) tends to reproduce y 1 with a larger probability than the initial distribution N(0, I);
the distribution N(0, C (2) ) tends to reproduce y 2 with a larger probability than N(0, C (1) ),
and so forth. When y 1 , . . . , y g denote the formerly selected, favorable steps, N(0, C (g) )
tends to reproduce these steps. The process leads to an alignment of the search distribution
N(0, C (g) ) to the distribution of the selected steps. If both distributions become alike, as
under random selection, in expectation no further change of the covariance matrix takes place
[7].

15
3.3.2 Cumulation: Utilizing the Evolution Path
(g+1) (g+1)
We have used the selected steps, y i:λ = (xi:λ − m(g) )/σ (g) , to update the covariance
matrix in (14) and (20). Because yy T = −y(−y)T , the sign of the steps is irrelevant for
the update of the covariance matrix — that is, the sign information is not used for calculating
C (g+1) . To exploit the sign information, the so-called evolution path is introduced [14, 16].
We call a sequence of successive steps, the strategy takes over a number of generations,
an evolution path. An evolution path can be expressed by a sum of consecutive steps. This
summation is referred to as cumulation. To construct an evolution path, the step-size σ is
disregarded. For example, an evolution path of three steps of the distribution mean m can be
constructed by the sum

m(g+1) − m(g) m(g) − m(g−1) m(g−1) − m(g−2)

(g)
+ (g−1)
+ . (21)
σ σ σ (g−2)
In practice, to construct the evolution path, pc ∈ Rn , we use exponential smoothing as in (14),
(0)
and start with pc = 0.15
p m(g+1) − m(g)
p(g+1)
c = (1 − cc )p(g)
c + cc (2 − cc )µeff (22)
σ (g)
where
(g)
pc ∈ Rn , evolution path at generation g.
cc ≤ 1. Again, 1/cc is the backward time horizon of the evolution path pc that contains
roughly
√ 63% of the overall weight (compare derivation of (18)). A time horizon between
n and n is reasonable.

p
The factor cc (2 − cc )µeff is a normalization constant for pc . For cc = 1 and µeff = 1, the
(g+1) (g+1)
factor reduces to one, and pc = (x1:λ − m(g) )/σ (g) .
p
The factor cc (2 − cc )µeff is chosen, such that

pc(g+1) ∼ N (0, C) (23)

if
(g+1)
xi:λ − m(g)
p(g)
c ∼ ∼ N (0, C) for all i = 1, . . . , µ . (24)
σ (g)
To derive (23) from (24) and (22) remark that
µ
p 2 X 1
(1 − cc )2 + cc (2 − cc ) = 1 and wi Ni (0, C) ∼ √ N(0, C) . (25)
i=1
µ eff

(g+1)
The (rank-one) update of the covariance matrix C (g) via the evolution path pc reads [14]
T
C (g+1) = (1 − c1 )C (g) + c1 p(g+1)
c p(g+1)
c . (26)
15 In the final algorithm (22) is still slightly modified, compare (42).

16
An empirically validated choice for the learning rate in (26) is c1 ≈ 2/n2 . For cc = 1 and
µ = 1, Equations (26), (20), and (14) are identical.
Using the evolution path for the update of C is a significant improvement of (14) for
small µeff , because correlations between consecutive steps are exploited. The leading signs of
steps, and the dependencies between consecutive steps play a significant role for the resulting
(g+1)
evolution path pc . For cc ≈ 3/n the number of function evaluations needed to adapt a
nearly optimal covariance matrix on cigar-like objective functions becomes O(n).
As a last step, we combine (14) and (26).

3.4 Combining Rank-µ-Update and Cumulation

The final CMA update of the covariance matrix combines (14) and (26).

C (g+1) = (1 − c1 − cµ ) C (g)
µ T
(g+1) T

(g+1) (g+1)
X
+ c1 p(g+1) p + cµ w i y y (27)
| c {z c } i:λ i:λ
i=1
rank-one update | {z }
rank-µ update

where
c1 ≈ 2/n2 .
cµ ≈ min(µeff /n2 , 1 − c1 ).
(g+1) (g+1)
y i:λ = (xi:λ − m(g) )/σ (g) .

Equation (27) reduces to (14) for c1 = 0 and to (26) for cµ = 0. The equation combines
the advantages of (14) and (26). On the one hand, the information within the population of
one generation is used efficiently by the rank-µ update. On the other hand, information of
correlations between generations is exploited by using the evolution path for the rank-one
update. The former is important in large populations, the latter is in particular important in
small populations.

4 Step-Size Control
The covariance matrix adaptation, introduced in the last section, does not explicitly control the
“overall scale” of the distribution, the step-size. The covariance matrix adaptation increases
the scale only in one direction for each selected step, and it decreases the scale only implicitly
by fading out old information via the factor 1 − c1 − cµ . Less informally, we find two specific
reasons to introduce a step-size control in addition to the adaptation rule (27) for C (g) .
1. The optimal overall step length cannot be well approximated by (27), in particular if
µeff is chosen larger than one.

17
Figure 5: Three evolution paths of respectively six steps from different selection situations
(idealized). The lengths of the single steps are all comparable. The length of the evolution
paths (sum of steps) is remarkably different and is exploited for step-size control

Pn 2
For example, on fsphere (x) = i=1 xi , the optimal step-size σ equals approxi-
p (g)
mately µ fsphere (x)/n, given C ≈ I and µeff = µ  n [3, 20]. This depen-
dency on µ cannot be realized by (14), and is also not well approximated by (27).

2. The largest reliable learning rate for the covariance matrix update in (27) is too slow to
achieve competitive change rates for the overall step length.
To achieve optimal performance on fsphere with an Evolution Strategy, the overall
step length must decrease by a factor of approximately exp(0.202) ≈ 1.22 within n
function evaluations, as can be derived from progress formulas [3, p. 229]. That is,
the time horizon for the step length change must be proportional to n or shorter. From
the learning rates c1 and cµ in (27) follows that the adaptation is too slow to perform
competitive on fsphere whenever µeff  n. This can be validated by simulations
even for moderate dimensions, say, n ≥ 10 and small µeff , say, µeff ≤ 1 + ln n.

To control the step-size σ (g) we utilize an evolution path, i.e. a sum of successive steps (see
page 16). The method can be applied independently of the covariance matrix update and is
denoted as cumulative path length control, cumulative step-size control, or cumulative step
length adaptation (CSA). The length of an evolution path is exploited, based on the following
reasoning (compare also Fig. 5).

• Whenever the evolution path is short, single steps cancel each other out (Fig. 5, left).
Loosely speaking, they are anti-correlated. If steps annihilate each other, the step-size
should be decreased.
• Whenever the evolution path is long, the single steps are pointing to similar directions
(Fig. 5, right). Loosely speaking, they are correlated. Because the steps are similar, the
same distance can be covered by fewer but longer steps into the same directions. In the
limit case, where consecutive steps have identical direction, they can be replaced by an
enlarged single step. Consequently, the step-size should be increased.

18
 • Subsuming, in the desired situation the steps are (approximately) perpendicular in ex-
pectation and therefore uncorrelated (Fig. 5, middle).
To decide whether the evolution path is “long” or “short”, we compare the length of the path
with its expected length under random selection.16 Under random selection consecutive steps
are independent and therefore uncorrelated (we just realized that “uncorrelated” steps are the
desired situation). If selection biases the evolution path to be longer then expected, σ is in-
creased, and, vice versa, if selection biases the evolution path to be shorter than expected, σ is
decreased. In the ideal situation, selection does not bias the length of the evolution path and
the length equals its expected length under random selection.
In practice, to construct the evolution path, pσ , the same techniques as in (22) are applied.
In contrast to (22), a conjugate evolution path is constructed, because the expected length
of the evolution path pc from (22) depends on its direction (compare (23)). Initialized with
(0)
pσ = 0, the conjugate evolution path reads

p − 12 m(g+1) − m(g)
p(g+1)
σ = (1 − cσ )p(g)
σ + cσ (2 − cσ )µeff C (g) (28)
σ (g)

where
(g)
pσ ∈ Rn is the conjugate evolution path at generation g.
cσ < 1. Again, 1/cσ is the backward time √ horizon of the evolution path (compare (18)). For
small µeff , a time horizon between n and n is reasonable.
p
cσ (2 − cσ )µeff is a normalization constant, see (22).
− 1 def −1 T
2 T
C (g) 2 = B (g) D (g) B (g) , where C (g) = B (g) D (g) B (g) is an eigendecompo-
sition of C (g) , where B (g) is an orthonormal basis of eigenvectors, and the diagonal
elements of the diagonal matrix D (g) are square roots of the corresponding positive
eigenvalues (cf. Sect. 0.1).

−1 − 21
For C (g) = I, we have C (g) 2 = I and (28) replicates (22). The transformation C (g)
re-scales the step m(g+1) − m(g) within the coordinate system given by B (g) .
−1 −1 T
The single factors of the transformation C (g) 2
= B (g) D (g) B (g) can be explained
as follows (from right to left):
T
B (g) rotates the space such that the columns of B (g) , i.e. the principal axes of the
distribution N(0, C (g) ), rotate into the coordinate axes. Elements of the resulting
vector relate to projections onto the corresponding eigenvectors.
−1
D (g) applies a (re-)scaling such that all axes become equally sized.
(g)
B rotates the result back into the original coordinate system. This last transforma-
tion ensures that the principal axes of the distribution are not rotated by the overall
transformation and directions of consecutive steps are comparable.
16 Random (g+1)
selection means that the index i : λ (compare (6)) is independent of the value of xi:λ for all i =
1, . . . , λ, e.g. i : λ = i.

19
 −1 (g+1)
Consequently, the transformation C (g) 2 makes the expected length of pσ independent of
(g)
its direction, and for any sequence of realized covariance matrices C g=0,1,2,... we have under
(g+1) (0)
random selection pσ ∼ N (0, I), given pσ ∼ N (0, I) [7].
(g+1)
To update σ (g) , we “compare” kpσ k with its expected length EkN (0, I) k, that is
cσ  
ln σ (g+1) = ln σ (g) + kp(g+1)
σ k − EkN (0, I) k , (29)
dσ EkN (0, I) k

where
dσ ≈ 1, damping parameter, scales the change magnitude of ln σ (g) . The factor cσ /dσ /EkN (0, I) k
is based on in-depth investigations of the algorithm [7].
√ √
EkN (0, I) k = 2 Γ( n+1 n
2 )/Γ( 2 ) ≈ n + O(1/n), expectation of the Euclidean norm of a
N (0, I) distributed random vector.

(g+1)
For kpσ k = EkN (0, I) k the second summand in (29) is zero, and σ (g) is unchanged,
(g+1) (g+1)
while σ (g) is increased for kpσ k > EkN (0, I) k, and σ (g) is decreased for kpσ k<
EkN (0, I) k.
(g+1)
Alternatively, we might use the squared norm kpσ k2 in (29) and compare with its
expected value n [2]. In this case (29) would read
cσ  (g+1) 2 
ln σ (g+1) = ln σ (g) + kpσ k −n
dσ 2n
!
(g+1) 2
(g) cσ kpσ k
= ln σ + −1 . (30)
2dσ n

This update will presumable lead to faster step-size increments and slower step-size decre-
ments.
  
The step-size change is unbiased on the log scale, because E ln σ (g+1) σ (g) = ln σ (g)
(g+1)
for pσ ∼ N (0, I). The role of unbiasedness is discussed in Sect. 5. Equations (28)
−1
and (29) cause successive steps of the distribution mean m(g) to be approximately C (g) -
conjugate.

−1
In order to show that successive steps are approximately C (g) -conjugate first we re-
(g+1)
mark that (28) and (29) adapt σ such that the length of pσ equals approximately
(g+1) 2 (g+1) T (g+1)
EkN (0, I) k. Starting from (EkN (0, I) k)2 ≈ kpσ k = pσ pσ =
1
(g) − 2
RHST RHS of (28) and assuming that the expected squared length of C (m(g+1) −
m(g) ) is unchanged by selection (unlike its direction) we get
T −1
p(g)
σ C (g) 2
(m(g+1) − m(g) ) ≈ 0 , (31)
and  T
1 −1
 
C (g) 2 p(g)
σ C (g) m(g+1) − m(g) ≈ 0 . (32)

20
 (g−1) T −1
Given 1/(c1 + cµ )  1 and (31) we assume also pσ C (g) 2 (m(g+1) − m(g) ) ≈ 0
and derive T −1
  
m(g) − m(g−1) C (g) m(g+1) − m(g) ≈ 0 . (33)
−1
That is, the steps taken by the distribution mean become approximately C (g) -conjugate.

Because σ (g) > 0, (29) is equivalent to

!!
(g+1)
(g+1) (g) cσ kpσ k
σ =σ exp −1 (34)
dσ EkN (0, I) k

The length of the evolution path is an intuitive and empirically well validated goodness
measure for the overall step length. For µeff > 1 it is the best measure to our knowledge.
Nevertheless, it fails to adapt nearly optimal step-sizes on very noisy objective functions [4].

5 Discussion
The CMA-ES is an attractive option for non-linear optimization, if “classical” search meth-
ods, e.g. quasi-Newton methods (BFGS) and/or conjugate gradient methods, fail due to a
non-convex or rugged search landscape (e.g. sharp bends, discontinuities, outliers, noise, and
local optima). Learning the covariance matrix in the CMA-ES is analogous to learning the
inverse Hessian matrix in a quasi-Newton method. In the end, any convex-quadratic (ellip-
soid) objective function is transformed into the spherical function fsphere . This can improve
the performance on ill-conditioned and/or non-separable problems by orders of magnitude.
The CMA-ES overcomes typical problems that are often associated with evolutionary al-
gorithms.
1. Poor performance on badly scaled and/or highly non-separable objective functions.
Equation (27) adapts the search distribution to badly scaled and non-separable prob-
lems.

2. The inherent need to use large population sizes. A typical, however intricate to diagnose
reason for the failure of population based search algorithms is the degeneration of the
population into a subspace. This is usually prevented by non-adaptive components in
the algorithm and/or by a large population size (considerably larger than the problem
dimension). In the CMA-ES, the population size can be freely chosen, because the
learning rates c1 and cµ in (27) prevent the degeneration even for small population sizes,
e.g. λ = 9. Small population sizes usually lead to faster convergence, large population
sizes help to avoid local optima.
3. Premature convergence of the population. Step-size control in (34) prevents the pop-
ulation to converge prematurely. It does not prevent the search to end up in a local
optimum.

21
 Therefore, the CMA-ES is highly competitive on a considerable number of test functions
[7, 11, 13, 15, 16] and was successfully applied to many real world problems.17 Finally we
discuss a few basic design principles that were applied in the previous sections.

Change rates We refer to a change rate as the expected parameter change per sampled
search point, given a certain selection situation. To achieve competitive performance on a
wide range of objective functions, the possible change rates of the adaptive parameters need
to be adjusted carefully. The CMA-ES separately controls change rates for the mean value of
the distribution, m, the covariance matrix, C, and the step-size, σ.

• The change rate for the mean value m, given a fixed sample distribution, is determined
by the parent number and the recombination weights. The larger µeff , the smaller is the
possible change rate of m. Similar holds for most evolutionary algorithms.

• The change rate of the covariance matrix C is explicitly controlled by the learning rates
c1 and cµ and therefore detached from parent number and population size. The learning
rate reflects the model complexity. In evolutionary algorithms, the explicit control of
change rates of covariances, independent from the population size, is a rare feature.
• The change rate of the step-size σ is explicitly controlled by the damping parameter dσ
and is in particular independent from the change rate of C. The time constant 1/cσ ≤ n
ensures a sufficiently fast change of the overall step length in particular with small
population sizes.

Invariance Invariance properties of a search algorithm denote identical behavior on a set, or

a class of objective functions. Invariance is an important property of the CMA-ES.18 Trans-
lation invariance should be taken for granted in continuous domain optimization. Translation
invariance means that the search behavior on the function x 7→ f (x + a), x(0) = b − a,
is independent of a ∈ Rn . Further invariances, e.g. to certain linear transformations of the
search space, are highly desirable: they imply uniform performance on classes of functions
and therefore allow for generalization of empirical results. The CMA-ES exhibits the follow-
ing invariances.
• Invariance to order preserving (i.e. strictly monotonic) transformations of the objective
function value. The algorithm only depends on the ranking of function values.
• Invariance to angle preserving (rigid) transformations of the search space (rotation, re-
flection, and translation) if the initial search point is transformed accordingly.
• Scale invariance if the initial scaling, e.g. σ (0) , and the initial search point, m(0) , are
chosen accordingly.
• Invariance to a scaling of variables if the initial diagonal covariance matrix C (0) , and
the initial search point, m(0) , are chosen accordingly.
17 For a list of published applications see http://www.lri.fr/˜hansen/cmaapplications.pdf
18 Special acknowledgments to Iván Santibán̄ez-Koref for pointing this out to me.

22
 • Invariance to any invertible linear transformation of the search space, A, if the initial

T
covariance matrix C (0) = A−1 A−1 , and if the initial search point, m(0) , are trans-
formed accordingly.
Invariance should be a fundamental design criterion for any search algorithm. Together with
the ability to efficiently adapt the invariance governing parameters, invariance is a key to
competitive performance.

Stationarity An important design criterion for a stochastic search procedure is unbiasedness

of variations of object and strategy parameters [5, 16]. Consider random selection, e.g. the
objective function f (x) = rand to be independent of x. The population
 mean is unbiased
 if its
expected value remains unchanged in the next generation, that is E m(g+1)  m(g) = m(g) .
For the population mean, stationarity under random selection is a rather intuitive concept.
In the CMA-ES, stationarity is respected for all parameters in the basic equation (5). The
distribution mean m, the covariance matrix C, and ln σ are unbiased. Unbiasedness
  of ln σ
does not imply that σ is unbiased. Actually, under random selection, E σ (g+1)  σ (g) > σ (g) ,
compare (29).19
For distribution variances (or step-sizes) a bias toward increase or decrease entails the
danger of divergence or premature convergence, respectively, whenever the selection pressure
is low. Nevertheless, on noisy problems a properly controlled bias towards increase might be
appropriate even on the log scale.

Acknowledgments
The author wishes to gratefully thank Anne Auger, Christian Igel, Stefan Kern, and Fabrice
Marchal for the many valuable comments on the manuscript.

References
[1] Auger A, Hansen N. A restart CMA evolution strategy with increasing population size.
In Proceedings of the IEEE Congress on Evolutionary Computation, 2005.
[2] Arnold DV, Beyer HG. Performance analysis of evolutionary optimization with cumu-
lative step length adaptation. IEEE Transactions on Automatic Control, 49(4):617–622,
2004.
[3] Beyer HG. The Theory of Evolution Strategies. Springer, Berlin, 2001.
[4] Beyer HG, Arnold DV. Qualms regarding the optimality of cumulative path length con-
trol in CSA/CMA-evolution strategies. Evolutionary Computation, 11(1):19–28, 2003.

[5] Beyer HG, Deb K. On self-adaptive features in real-parameter evolutionary algorithms.

IEEE Transactions on Evolutionary Computation, 5(3):250–270, 2001.
19 (g+1) ,
 Alternatively,
  (34) would be designed to be unbiased for σ
if this would imply that
E ln σ (g+1)  σ (g) < ln σ (g) , to our opinion a less desirable variant.

23
 [6] G. Collange, N. Delattre, N. Hansen, I. Quinquis, and M. Schoenauer. Multidisciplinary
optimisation in the design of future space launchers. In P. Breitkopf and R. F. Coelho,
editors, Multidisciplinary Design Optimization in Computational Mechanics, chapter 12,
pages 487–496. Wiley, 2010.
[7] Hansen N. Verallgemeinerte individuelle Schrittweitenregelung in der Evolutionsstrate-
gie. Mensch und Buch Verlag, Berlin, 1998.
[8] Hansen N. Invariance, self-adaptation and correlated mutations in evolution strategies.
In Schoenauer M, Deb K, Rudolph G, Yao X, Lutton E, Merelo JJ, Schwefel HP, editors,
Parallel Problem Solving from Nature - PPSN VI, pages 355–364. Springer, 2000.
[9] Hansen N. The CMA evolution strategy: a comparing review. In Lozano JA, Larranaga P,
Inza I, and Bengoetxea E, editors, Towards a new evolutionary computation. Advances
on estimation of distribution algorithms, pages 75–102. Springer, 2006.
[10] Hansen N. Benchmarking a BI-Population CMA-ES on the BBOB-2009 Function
Testbed. In the workshop Proceedings of the GECCO Genetic and Evolutionary Com-
putation Conference, pages 2389–2395. ACM, 2009.
[11] Hansen N, Kern S. Evaluating the CMA evolution strategy on multimodal test functions.
In Xin Yao et al., editors, Parallel Problem Solving from Nature - PPSN VIII, pages
282–291. Springer, 2004.
[12] Hansen N, Niederberger SPN, Guzzella L, Koumoutsakos P. A method for handling un-
certainty in evolutionary optimization with an application to feedback control of com-
bustion. IEEE Transactions on Evolutionary Computation, 2009.
[13] Hansen N, Müller SD, Koumoutsakos P. Reducing the time complexity of the deran-
domized evolution strategy with covariance matrix adaptation (CMA-ES). Evolutionary
Computation, 11(1):1–18, 2003.
[14] Hansen N, Ostermeier A. Adapting arbitrary normal mutation distributions in evolution
strategies: The covariance matrix adaptation. In Proceedings of the 1996 IEEE Confer-
ence on Evolutionary Computation (ICEC ’96), pages 312–317, 1996.
[15] Hansen N, Ostermeier A. Convergence properties of evolution strategies with the deran-
domized covariance matrix adaptation: The (µ/µI , λ)-CMA-ES. In Proceedings of the
5th European Congresson Intelligent Techniques and Soft Computing, pages 650–654,
1997.
[16] Hansen N, Ostermeier A. Completely derandomized self-adaptation in evolution strate-
gies. Evolutionary Computation, 9(2):159–195, 2001.
[17] Kern S, Müller SD, Hansen N, Büche D, Ocenasek J, Koumoutsakos P. Learning proba-
bility distributions in continuous evolutionary algorithms – a comparative review. Natu-
ral Computing, 3:77–112, 2004.
[18] Larrañaga P. A review on estimation of distribution algorithms. In P. Larrañaga and J. A.
Lozano, editors, Estimation of Distribution Algorithms, pages 80–90. Kluwer Academic
Publishers, 2002.

24
[19] Larrañaga P, Lozano JA, Bengoetxea E. Estimation of distribution algorithms based
on multivariate normal and Gaussian networks. Technical report, Dept. of Computer
Science and Artificial Intelligence, University of the Basque Country, 2001. KZAA-IK-
1-01.
[20] Rechenberg I. Evolutionsstrategie ’94. Frommann-Holzboog, Stuttgart, Germany, 1994.
[21] Rubenstein RY, Kroese DP. The Cross-Entropy Method: a unified approach to combina-
torial optimization, Monte-Carlo simulation, and machine learning. Springer, 2004.

A Algorithm Summary: The (µ/µW , λ)-CMA-ES

Figure 6 outlines the complete algorithm, summarizing (5), (6), (22), (27), (28), and (34).
Used symbols are:
y k ∼ N (0, C), for k = 1, . . . , λ, are realizations from a multivariate normal distribution
with zero mean and covariance matrix C.
B, D result from an eigendecomposition of the covariance matrix C with C = BD 2 B T =
BDDB T (cf. Sect. 0.1). Columns of B are an orthonormal basis of eigenvectors. Di-
agonal elements of the diagonal matrix D are square roots of the corresponding positive
eigenvalues. While (36) can certainly be implemented using a Cholesky decomposition
of C the eigendecomposition is needed to compute (40) correctly.
xk ∈ Rn , for k = 1, . . . , λ. Sample of λ search points.
Pµ
hyiw = i=1 wi y i:λ , step of the distribution mean disregarding step-size σ.
y i:λ = (xi:λ − m)/σ, see xi:λ below.
xi:λ ∈ Rn , i-th best point out of x1 , . . . , xλ from (37). The index i : λ denotes the index of
the i-th ranked point, that is f (x1:λ ) ≤ f (x2:λ ) ≤ · · · ≤ f (xλ:λ ).
2 −1
Pµ

µeff = i=1 wi is the variance effective selection mass, see (8). It holds 1 ≤ µeff ≤ µ.
1 def
C − 2 = BD −1 B T , see B, D above. The matrix D can be inverted by inverting its diag-
− 12
Pµ elements. From the definitions we find that C hyiw = B hziw with hziw =
onal
i=1 wi z i:λ .
√ √
EkN (0, I) k = 2 Γ( n+1 n 1 1


2 )/Γ( 2 ) ≈ n 1 − 4n + 21n 2 .

kpσ k
(
2
1 if √ < (1.4 + n+1 )EkN (0, I) k
hσ = 1−(1−cσ )2(g+1) , where g is the generation
0 otherwise
number. The Heaviside function hσ stalls the update of pc in (42) if kpσ k is large.
This prevents a too fast increase of axes of C in a linear surrounding, i.e. when the
step-size is far too small. This is useful when the initial step-size is chosen far too small
or when the objective function changes in time.
δ(hσ ) = (1 − hσ )cc (2 − cc ) ≤ 1 is of minor relevance. In the (unusual) case of hσ = 0, it
substitutes for the second summand from (42) in (43).

25
Set parameters

Set parameters λ, µ, wi=1...µ , cσ , dσ , cc , c1 , and cµ to their default values according to

Table 1.
Initialization
Set evolution paths pσ = 0, pc = 0, covariance matrix C = I, and g = 0.

Choose distribution mean m ∈ Rn and step-size σ ∈ R+ problem dependent.1

Until termination criterion met, g ← g + 1
Sample new population of search points, for k = 1, . . . , λ

zk ∼ N (0, I) (35)
yk = BDz k ∼ N (0, C) (36)
2


xk = m + σy k ∼ N m, σ C (37)

Selection and recombination

µ
X µ
X
hyiw = wi y i:λ where wi = 1, wi > 0 (38)
i=1 i=1
µ
X
m ← m + σ hyiw = wi xi:λ (39)
i=1

Step-size control
1
← (1 − cσ )pσ + cσ (2 − cσ )µeff C − 2 hyiw
p
pσ (40)
  
cσ kpσ k
σ ← σ × exp −1 (41)
dσ EkN (0, I) k

Covariance matrix adaptation

p
pc ← (1 − cc )pc + hσ cc (2 − cc )µeff hyiw (42)
µ
X
← (1 − c1 − cµ ) C + c1 pc pT wi y i:λ y T


C c + δ(hσ )C + cµ i:λ (43)
i=1

1
The optimum should presumably be within the initial cube m ± 3σ(1, . . . , 1)T . If the optimum is ex-
pected to be in the initial search interval [a, b]n we may choose the initial search point, m, uniformly randomly
in [a, b]n , and σ = 0.3(b − a). Different search intervals ∆si for different variables can be reflected by a
different initialization of C, in that the diagonal elements of C obey cii = (∆si )2 . Remark that the ∆si
should not disagree by several orders of magnitude. Otherwise a scaling of the variables should be applied.

Figure 6: The (µ/µW , λ)-CMA Evolution Strategy. Symbols: see text

26
 Pµ
Table 1: Default Strategy Parameters, where µeff = Pµ 1 2 ≥ 1 and wi = 1, taken
i=1 wi
i=1
from [10]
Selection and Recombination:
λ
λ = 4 + b3 ln nc, µ = bµ0 c, µ0 = (44)
2
w0
w i = Pµ i 0 , wi0 = ln(µ0 + 0.5) − ln i for i = 1, . . . , µ, (45)
j=1 wj

Step-size control:
r !
µeff + 2 µeff − 1
cσ = , dσ = 1 + 2 max 0, − 1 + cσ (46)
n + µeff + 5 n+1

Covariance matrix adaptation:

4 + µeff /n
cc = (47)
n + 4 + 2µeff /n
2
c1 = (48)
(n + 1.3)2 + µeff
 
µeff − 2 + 1/µeff
cµ = min 1 − c1 , αµ with αµ = 2 (49)
(n + 2)2 + αµ µeff /2

Default Parameters The (external) strategy parameters are λ, µ, wi=1...µ , cσ , dσ , cc , c1 ,

and cµ . Default strategy parameters values are given in Table 1. An in-depth discussion of
most parameters is given in [16]. The default parameters of (46)-(49) are in particular chosen
to be a robust setting and therefore, to our experience, applicable to a wide range of functions
to be optimized. We do not recomment to change this setting. In contrast, the population
size λ in (44) can be increased.20 If the λ-dependent default values for µ and wi are used, the
population size λ has a significant influence on the global search performance [11]. Increasing
λ usually improves the global search capability and the robustness of the CMA-ES, at the
price of a reduced convergence speed. The convergence speed decrease at most linearly with
λ. Independent restarts with increasing population size [1], automated or manually conducted,
are an appropriate policy to perform well on most problems.

B Implementational Concerns
We discuss a few implementational concerns.
20 Decreasing λ is not recommended. Too small values regularly have strong adverse effects on the performance.

27
B.1 Multivariate normal distribution
Let the vector z ∼ N (0, I) have independent, (0, 1)-normally distributed components that
can easily be sampled on a computer. To generate a random vector y ∼ N(0, C) for (36),
we set y = BDz (see above symbol descriptions of B and D and Sects. 0.1 and 0.2, and
1
compare lines 52–53 and 83–84 in thePsource code below). Given y k = BDz k and C − 2 =
1 µ
BD −1 B T we have C − 2 hyiw = B i=1 wi z i:λ (compare (40) and lines 61 and 64 in the
source code below).

B.2 Strategy internal numerical effort

In practice, the re-calculation of B and D needs to be done not until max(1, b1/(10n(c1 +
cµ ))c) generations. For reasonable c1 + cµ values, this reduces the numerical effort due to the
eigendecomposition from O(n3 ) to O(n2 ) per generated search point, that is the effort of a
matrix vector multiplication.
On a Pentium 4, 2.5 GHz processor the overall strategy internal time consumption is
roughly 3 × (n + 5)2 × 10−8 seconds per function evaluation [17].
Remark that it is not sufficient to compute a Cholesky decomposition of C, because then
(40) cannot be computed correctly.

B.3 Termination criteria

In general, the algorithm should be stopped whenever it becomes a waste of CPU-time to con-
tinue, and it would be better to restart (eventually with increased population size [1]) or to
reconsidering the encoding and/or objective function formulation. We recommend the follow-
ing termination criteria [1, 10] that are mostly related to numerical stability:

• NoEffectAxis: stop if adding a 0.1-standard deviation vector in any principal axis

direction of C does not change m.21
• NoEffectCoord: stop if adding 0.2-standard deviations in any single coordinate does
not change m (i.e. mi equals mi + 0.2 σci,i for any i).

• EqualFunValues: stop if the range of the best objective function values of the last
10 + d30n/λe generations is zero.
• Stagnation: we track a history of the best and the median fitness in each iteration
over the last 20% but at least 120 + 30n/λ and no more than 20 000 iterations. We stop,
if in both histories the median of the last (most recent) 30% values is not better than the
median of the first 30%.
• ConditionCov: stop if the condition number of the covariance matrix exceeds 1014 .
• TolXUp: stop if σ ×max(diag(D)) increased by more than 104 . This usually indicates
a far too small initial σ, or divergent behavior.
21 More formally, we terminate if m equals to m + 0.1 σd b , where i = (g mod n) + 1, and d2 and b are
ii i ii i
respectively the i-th eigenvalue and eigenvector of C, with kbi k = 1.

28
 Two other useful termination criteria should be considered problem dependent:

• TolFun: stop if the range of the best objective function values of the last 10+d30n/λe
generations and all function values of the recent generation is below TolFun. Choosing
TolFun depends on the problem, while 10−12 is a conservative first guess.
• TolX: stop if the standard deviation of the normal distribution is smaller than in all
coordinates and σpc is smaller than TolX in all components. By default we set TolX
to 10−12 times the initial σ.

B.4 Flat fitness

In the case of equal function values for several individuals in the population, it is feasible
to increase the step-size (see lines 92–96 in the source code below). This method can inter-
fere with the termination criterion TolFun. In practice, observation of a flat fitness should
be rather a termination criterion and consequently lead to a reconsideration of the objective
function formulation.

B.5 Boundaries and Constraints.

The handling of boundaries and constraints is to a certain extend problem dependent. We
discuss a few useful approaches.

Best solution is stricly inside If the optimal solution is not too close to the infeasible do-
main (or at least not in a “corner”), a simple and sufficient way to handle any type of
boundaries and constraints are

1. set the fitness as

ffitness (x) = fmax + kx − xfeasible k , (50)
where fmax is larger than the worst fitness in the feasible population or in the
feasible domain (in case of minization) and xfeasible is a constant feasible point,
preferably in the middle of the feasible domain.
2. re-sampling any infeasible solution x until it become feasible.

Repairement available as for example with box-constraints.

Simple repairment It is possible to simply repair infeasible individuals before the up-
date equations are applied. This is not recommended, because the CMA-ES makes
implicite assumptions on the distribution of solution points, which can be heavily
violated by a repairment. The main resulting problem might be divergence or too
fast convergence of the step-size (however a solution will soon be available). Note
also, that this might be intricate to implement, in particular if y or z are used for
implementing the update equations in the original code.
Penalization We evaluate the objective function on a repaired search point, xrepaired ,
and add a penalty depending on the distance to the repaired solution.

ffitness (x) = f (xrepaired ) + α kx − xrepaired k2 . (51)

29
 The repaired solution is disregarded afterwards.
In case of box-boundaries, xrepaired is set to the feasible solution with the smallest
distance kx − xrepaired k. In other words, components that are infeasible in x are
set to the (closest) boundary value in xrepaired . A similar boundary handling with
a component-wise adaptive α is described in [12].

No repair mechanism available The fitness of the infeasible search point x might similarly
compute to X
ffitness (x) = foffset + α 11ci >0 × ci (x)2 (52)
i

where, w.l.o.g., the (non-linear) constraints ci : Rn → R, x 7→ ci (x) are satisfied

for ci (x) ≤ 0 , and the indicator function 11ci >0 equals to one for ci (x) > 0, zero
otherwise, and foffset = mediank f (xk ) equals, for example, to the median or 25%-tile
function value of the feasible points in the same generation. If no other information
is available, ci (x) might be computed as the squared distance of x to the best feasible
solution in the population. This approach has not yet been experimentally evaluated by
the author. A different, slightly more involved approach is given in [6].

In either case of (51) and (52), α should be chosen such that the differences in f and the
differences in the second summand have a similar magnitude.

30
C MATLAB Source Code
1 function xmin=purecmaes
2 % CMA-ES: Evolution Strategy with Covariance Matrix Adaptation for
3 % nonlinear function minimization.
4 %
5 % This code is an excerpt from cmaes.m and implements the key parts
6 % of the algorithm. It is intendend to be used for READING and
7 % UNDERSTANDING the basic flow and all details of the CMA *algorithm*.
8 % Computational efficiency is sometimes disregarded.
9
10 % -------------------- Initialization --------------------------------
11
12 % User defined input parameters (need to be edited)
13 strfitnessfct = ’felli’; % name of objective/fitness function
14 N = 10; % number of objective variables/problem dimension
15 xmean = rand(N,1); % objective variables initial point
16 sigma = 0.5; % coordinate wise standard deviation (step-size)
17 stopfitness = 1e-10; % stop if fitness < stopfitness (minimization)
18 stopeval = 1e3*Nˆ2; % stop after stopeval number of function evaluations
19
20 % Strategy parameter setting: Selection
21 lambda = 4+floor(3*log(N)); % population size, offspring number
22 mu = lambda/2; % lambda=12; mu=3; weights = ones(mu,1); would be (3_I,12)-ES
23 weights = log(mu+1/2)-log(1:mu)’; % muXone recombination weights
24 mu = floor(mu); % number of parents/points for recombination
25 weights = weights/sum(weights); % normalize recombination weights array
26 mueff=sum(weights)ˆ2/sum(weights.ˆ2); % variance-effective size of mu
27
28 % Strategy parameter setting: Adaptation
29 cc = (4+mueff/N) / (N+4 + 2*mueff/N); % time constant for cumulation for C
30 cs = (mueff+2)/(N+mueff+5); % t-const for cumulation for sigma control
31 c1 = 2 / ((N+1.3)ˆ2+mueff); % learning rate for rank-one update of C
32 cmu = 2 * (mueff-2+1/mueff) / ((N+2)ˆ2+2*mueff/2); % and for rank-mu update
33 damps = 1 + 2*max(0, sqrt((mueff-1)/(N+1))-1) + cs; % damping for sigma
34
36 % Initialize dynamic (internal) strategy parameters and constants
37 pc = zeros(N,1); ps = zeros(N,1); % evolution paths for C and sigma
38 B = eye(N); % B defines the coordinate system
39 D = eye(N); % diagonal matrix D defines the scaling
40 C = B*D*(B*D)’; % covariance matrix
41 eigeneval = 0; % B and D updated at counteval == 0
42 chiN=Nˆ0.5*(1-1/(4*N)+1/(21*Nˆ2)); % expectation of
43 % ||N(0,I)|| == norm(randn(N,1))
44
45 % -------------------- Generation Loop --------------------------------
46
47 counteval = 0; % the next 40 lines contain the 20 lines of interesting code
48 while counteval < stopeval
49
50 % Generate and evaluate lambda offspring
51 for k=1:lambda,
52 arz(:,k) = randn(N,1); % standard normally distributed vector
53 arx(:,k) = xmean + sigma * (B*D * arz(:,k)); % add mutation % Eq. 37
54 arfitness(k) = feval(strfitnessfct, arx(:,k)); % objective function call
55 counteval = counteval+1;
56 end
57
58 % Sort by fitness and compute weighted mean into xmean
59 [arfitness, arindex] = sort(arfitness); % minimization
60 xmean = arx(:,arindex(1:mu))*weights; % recombination % Eq. 39
61 zmean = arz(:,arindex(1:mu))*weights; % == Dˆ-1*B’*(xmean-xold)/sigma
62
63 % Cumulation: Update evolution paths
64 ps = (1-cs)*ps + (sqrt(cs*(2-cs)*mueff)) * (B * zmean); % Eq. 40
65 hsig = norm(ps)/sqrt(1-(1-cs)ˆ(2*counteval/lambda))/chiN < 1.4+2/(N+1);
66 pc = (1-cc)*pc + hsig * sqrt(cc*(2-cc)*mueff) * (B*D*zmean); % Eq. 42
67

31
 68 % Adapt covariance matrix C
69 C = (1-c1-cmu) * C ... % regard old matrix % Eq. 43
70 + c1 * (pc*pc’ ... % plus rank one update
71 + (1-hsig) * cc*(2-cc) * C) ... % minor correction
72 + cmu ... % plus rank mu update
73 * (B*D*arz(:,arindex(1:mu))) ...
74 * diag(weights) * (B*D*arz(:,arindex(1:mu)))’;
75
76 % Adapt step-size sigma
77 sigma = sigma * exp((cs/damps)*(norm(ps)/chiN - 1)); % Eq. 41
78
79 % Update B and D from C
80 if counteval - eigeneval > lambda/(cone+cmu)/N/10 % to achieve O(Nˆ2)
81 eigeneval = counteval;
82 C=triu(C)+triu(C,1)’; % enforce symmetry
83 [B,D] = eig(C); % eigen decomposition, B==normalized eigenvectors
84 D = diag(sqrt(diag(D))); % D contains standard deviations now
85 end
86
87 % Break, if fitness is good enough
88 if arfitness(1) <= stopfitness
89 break;
90 end
91
92 % Escape flat fitness, or better terminate?
93 if arfitness(1) == arfitness(ceil(0.7*lambda))
94 sigma = sigma * exp(0.2+cs/damps);
95 disp(’warning: flat fitness, consider reformulating the objective’);
96 end
97
98 disp([num2str(counteval) ’: ’ num2str(arfitness(1))]);
99
100 end % while, end generation loop
101
102 % -------------------- Final Message ---------------------------------
103
104 disp([num2str(counteval) ’: ’ num2str(arfitness(1))]);
105 xmin = arx(:, arindex(1)); % Return best point of last generation.
106 % Notice that xmean is expected to be even
107 % better.
108
109 % ---------------------------------------------------------------
110 function f=felli(x)
111 N = size(x,1); if N < 2 error(’dimension must be greater one’); end
112 f=1e6.ˆ((0:N-1)/(N-1)) * x.ˆ2; % condition number 1e6

32
D Reformulation of Learning Parameter ccov
For sake of consistency and clarity, we have reformulated the learning coefficients in (43) and
replaced
ccov
with c1 (53)
µcov
 
1
ccov 1 − with cµ and (54)
µcov
1 − ccov with 1 − c1 − cµ , (55)

and choosen (in replacing (49))

2
c1 = (56)
(n + 1.3)2 + µcov
1
!
µcov − 2 + µcov
cµ = min 2 , 1 − c1 , (57)
(n + 2)2 + µcov

The resulting coefficients are quite similar to the previous. In contrast to the previous formu-
lation, c1 becomes monotonic in µ−1 eff and c1 + cµ becomes virtually monotonic in µeff .
Another alternative, depending only on the degrees of freedom in the covariance matrix
and additionally correcting for very small λ, reads

min(1, λ/6)
c1 = √ (58)
m + 2 m + µneff
αµ0 + µeff − 2 + µ1eff
!
cµ = min 1 − c1 , √ (59)
m + 4 m + µ2eff
αµ0 = 0.3 , (60)
2
where m = n 2+n is the degrees of freedom in the covariance matrix. For µeff = 1, the
coefficient cµ is now chosen to be larger than zero, as αµ0 > 0. Figure 7 compares the new
learning rates with the old ones.

33
Figure 7: Learning rates c1 , cµ (solid) and ccov (dash-dotted) versus µeff . Above: Equations
(56) etc. for n = 3; 10. Below: Equations (58) etc. for n = 2; 40. Black: c1 √ + cµ and
ccov ; blue: c1 and ccov /µcov ; green: cµ and (1 − 1/µcov )ccov ; cyan: 2/(n2 + 2); red:
(c1 + cµ )/ccov , above divided by ten. For µcov ≈ 2 the difference is maximal, because c1
decreases much slower with increasing µcov and ccov is non-monotonic in µcov (a main reason
for the new formulation).

34