Sample Reuse in the Covariance Matrix Adaptation
Evolution Strategy Based on Importance Sampling
Shinichi Shirakawa

University of Tsukuba
1-1-1 Tennodai, Tsukuba,
Ibaraki, Japan

Youhei Akimoto

Shinshu University
4-17-1 Wakasato, Nagano
City, Nagano, Japanshirakawa@iit.tsukuba.ac.jp y_akimoto@shinshuu.ac.jp
Kouzou Ohara

Kazuki Ouchi

Aoyama Gakuin University
5-10-1 Fuchinobe, Chuo-ku,
Sagamihara-shi, Kanagawa,Japan

kzk0520.cb@gmail.com

Aoyama Gakuin University
5-10-1 Fuchinobe, Chuo-ku,
Sagamihara-shi, Kanagawa,Japan

ohara@it.aoyama.ac.jp

ABSTRACT
Recent studies reveal that the covariance matrix adaptation evolution strategy (CMA-ES) updates the parameters based on the natural gradient. The rank-based weight is considered the result of thequantile-based transformation of the objective value and the parameters are adjusted in the direction of the natural gradient estimatedby Monte-Carlo with the samples drawn from the current distribution. In this paper, we propose a sample reuse mechanism for the
CMA-ES. On the basis of the importance sampling, the past samples are reused to reduce the estimation variance of the quantile andthe natural gradient. We derive the formula for the rank-µ updateof the covariance matrix and the mean vector update using the pastsamples, then incorporate it into the CMA-ES without the step-size
adaptation. From the numerical experiments, we observe that the
proposed approach helps to reduce the number of function evaluations on many benchmark functions, especially when the numberof samples at each iteration is relatively small.

Categories and Subject Descriptors
G.1.6 [Numerical Analysis]: Optimization—Global optimization

Keywords
Importance Sampling; Natural Gradient; Covariance Matrix Adaptation Evolution Strategy; Information Geometric Optimization

1. INTRODUCTION

The covariance matrix adaptation evolution strategy (CMA-ES)[6, 7] is a comparison-based stochastic search algorithm for continuous optimization. The CMA-ES controls the mean vector andPermission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others thanACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission
and/or a fee. Request permissions from permissions@acm.org.

GECCO ’15, July 11 - 15, 2015, Madrid, Spain© 2015 ACM. ISBN 978-1-4503-3472-3/15/07. . . $15.00
DOI: http://dx.doi.org/10.1145/2739480.2754704

the covariance matrix of the multivariate normal distribution, fromwhich candidate solutions are sampled. At each iteration it samples candidate solutions, then they are evaluated on the problem
and their ranking is computed. The parameter of the distribution
N (m, σ 2 C) are updated by using the candidate solutions and theirranking, so that the candidate solutions tend to converge towards
the optimum of the problem at a geometric rate. Since it does
not explicitly use the objective value, the CMA-ES is invariant
to the order-preserving transformation of the objective function,
f !→ g ◦ f , where f : Rd → R is an objective function andg : R → R is any increasing function. Invariances to the orderpreserving transformation of the objective function and to the linear
transformation of the search space enable the CMA-ES to eﬃciently
solve non-convex, non-separable and ill-conditioned problems [8].
Recent studies [2,5,10] reveal that some components of the CMAES, namely m-update and the rank-µ C-update, are regarded as thestochastic natural gradient ascent [4]. Ollivier et al. [10] haveproposed a generalized framework of the stochastic search method
for black-box optimization problems on an arbitrary search space,called information geometric optimization (IGO). The IGO framework includes not only the rank-µ update CMA-ES but also a binary optimization algorithm called population-based incremental
learning (PBIL) as instantiations with diﬀerent probability models. The IGO framework provides a mathematical interpretation of
the weighted recombination, based on which we develop a novel
technique in this paper.
In the IGO framework, the objective function (to be minimized,
w.l.o.g.) is transformed into the utility function on the basis of
the quantile of the objective value under the sampling distribution.
The objective of the parameter update is considered to maximize
the expected utility function over the sampling distribution, and theparameters are updated in its steepest ascent direction. With respect
to the KL-divergence, the steepest ascent direction is known to begiven by the so-called natural gradient [4], which is computed by
the product of the inverse Fisher information matrix and the gradient. The natural gradient is estimated by Monte-Carlo using thecandidate solutions sampled from the current distribution. Larger
sample size typically results in higher accuracy (lower variance) of
the natural gradient estimate, requires less iterations, but it costs
more function evaluations.

To accurate the natural gradient estimate without any additional
function evaluations and accelerate the search, we consider reusing
the candidate solutions of past iterations. In the CMA-ES and the
IGO setting, the samples drawn at past iterations are discarded. If wejust reuse the samples from past iterations without the modificationof the update, it will introduce a bias to the Monte-Carlo estimateof the natural gradient since the past distributions are diﬀerent fromthe current distribution. However, if the past distributions are close
to the current distribution, the past samples can be useful to estimate
the natural gradient. Reusing the past samples is natural and widelyused in the evolutionary computation, for instance, the (µ + λ)-ES,
though their purpose diﬀers from ours.
In this paper, we propose a sample reuse mechanism for the CMAES based on the importance sampling. We derive the formula toestimate the utility values (recombination weights) and the natural
gradient by using the samples from the mixture of the current andpast distributions. Theoretically it does not introduce any bias in theestimation and is expected to reduce the estimation variance. We
derive the formula of the rank-µ update of the covariance matrix
and the mean vector update with the sample reuse, then incorporate
rank-one update into this mechanism.
The rest of the paper is organized as follows. In Section 2 weoverview the CMA-ES, the IGO, and an existing sample reuse mechanism called importance mixing [13, 14]. In Section 3, we introduce the importance sampling in the IGO framework and derive
the update formula for the rank-µ update CMA-ES. In Section 4,
we evaluate the performance of the CMA-ES with and without the
proposed sample reuse mechanism. We also compare the proposed
strategies with the importance mixing. In Section 5, we summarize
the contribution and discuss the future work.

2. RELATED ALGORITHMS
The covariance matrix adaptation evolution strategy (CMA-ES)[6, 7] maintains the Gaussian distribution N (m, σ 2 C), from which
the candidate solutions are sampled. It repeats sampling, evaluation
and update of the parameters of the distribution. At each iteration
the CMA-ES generates the λ candidate solutions x i for i = 1, · · · , λ.
Their objective values are evaluated and the ranking of each solution is computed. The parameters of the Gaussian distribution is
updated by using the current candidate solutions and their ranking.
The Gaussian distribution in the CMA-ES is parameterized with amean vector m ∈ Rd and a covariance matrix σ 2 C, where σ ∈ R >is called global step-size and C ∈ Rd×d is a positive-definite symmetric matrix. These parameters are updated by the diﬀerent update
rules, the weighted recombination for m, the rank-µ and rank-one
update for C, and the cumulative step-size adaptation (CSA) for
σ. Here, we only explain the update rules of m and C since weconsider the CMA-ES without CSA to investigate the basic eﬀect ofthe proposed approach in this paper. Therefore, w.l.o.g., we assume
σ = 1.
Let rk(x it ) be the ranking of x it among the λ candidate solutions
with respect to f . The weight is assigned to each solution according
to its ranking. Each weight defined in [7] ismax(0, ln( λ+12 ) − ln(i))λ max(0, ln( λ+1 ) − ln( j))2
j=1

and the weight wrk(x t ) is assigned to each solution x it .
im t+1 = m t + cm

(1)

λ
"
i=1

wrk(x t ) (x it − m t ),i(2)where cm is the learning rate for m and usually cm = 1. The update
rule of the covariance matrix C consists of two components: rankone update and rank-µ update. The rank-one update accelerates the
update of the covariance matrix using the evolution path which isthe cumulation of the consecutive steps. The evolution path for the
rank-one update, pc , is updated with the following formulaλ"#pct+1 = (1 − cc )pct + cc (2 − cc ) µeﬀ
wrk(x t ) (x it − m t ), (3)ii=1where c!
c is the cumulation parameter for the evolution path andλ w 2 ) −1 . Then, the update rule of the covariance matrixµeﬀ = ( i=1iC with the rank-one and rank-µ update is
C t+1 = C t + c1 (pct+1 (pct+1 ) T − C t )$!!!!!!!!!!!!!!!!!!!!!!!%&!!!!!!!!!!!!!!!!!!!!!!!'
rank-one update

λ
"

(
)
wrk(x i ) (x it − m t )(x it − m t ) T − C t ,i=1$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!%&!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'

+ cµ

(4)

rank-µ update

where c1 and c µ are the learning rates of the rank-one and rankµ updates for C, respectively. The rank-one update enlarges the
eigenvalue of C corresponding to the direction of the evolution path,and the rank-µ update is considered the natural gradient ascent of
the expectation of the transformed objective function as describedin Section 2.2. The CMA-ES described above with c1 = 0 is called
the pure rank-µ update CMA-ES.

2.2

2.1 The CMA-ES without CSA

wi = !

The mean vector m is updated byInformation Geometric Optimization

Given a probability model pθ (x) on X with the parameter θ ∈Θ, the IGO [10] transforms the original minimization problem off : X → R to the maximization of Jθ t : Θ → R at each iterationand performs the so-called natural gradient ascent. The functionJθ t depends on the parameter θ t at each iteration. It is the expectedfvalue of the utility Wθ t (x) that is a nonlinear and non-increasingtransformation of the objective function f .
f
The utility function1, Wθ (x), is defined by using the quantileof f (x) under x ∼ pθ . Let w : [0, 1] → R be a non-increasing
function, qθ! (x) = Pθ [y : f (y) ! f (x)] represent the probabilityfof sampling a better point than x from pθ , and Wθ (x) = w(qθ! (x)).
Then, Jθ t (θ) is defined as the expectation of the utility over the
distribution pθ , namely,*fJθ t (θ) =Wθ t (x)pθ (x)dx.

˜ θ t (θ), is given by the product of the
The natural gradient, ∇J
inverse of the Fisher information matrix F (θ) and the vanilla gradient ∇Jθ t (θ), i.e., the vector of the partial derivative w.r.t. eachparameter, which can be written as*f
˜ θ t (θ) =˜ (θ; x)pθ (x)dx.
∇JWθ t (x) ∇l(5)
1 This definition assumes the probability of getting the same objective value among diﬀerent samples is 0. This assumption holds only
if the Lebesgue measure of the level set of the objective function is0, and it is no problem in most practical cases. The rigorous formulation that allows the same objective value is explained in [10].

˜ (θ; x) = F −1 (θ)∇l (θ; x) is the natural gradient of the logHere, ∇llikelihood l (θ; x) = ln pθ (x).
In practice, the integral in (5) cannot be computed analyticallysince the objective function is black-box, and have to be estimatedby Monte-Carlo using the samples drawn from pθ t (x). Replacing
the expectation in (5) with the arithmetic average over the samples
x i ∼ pθ t , we haveprobability distribution
pk . Our objective here is to estimate theexpectation g(x)p0 (x)dx by using λ(K + 1) samples x ki . In the
CMA-ES, p0 and pk corresponds to the current distribution and the
distribution of previous kth iteration. The simplest way is to apply
the Monte-Carlo estimate by using only the samples x 0i from p0 ,λ1"
g(x 0i )λλ"f˜ θ t (θ t ) ≈ 1
˜ (θ t ; x i ) .
∇JWθ t (x i ) ∇l
λ
i=1

f

We also need to estimate the utility Wθ t (x i ) = w(qθ!t (x i )) for eachx i . Since qθ!t (x i ) is the probability of sampling a better point
than x i from pθ t , it can be approximated by using the number ofbetter point than x i among the λ samples drawn from pθ t , namelyf(rk(x i ) − 1/2)/λ. Then, the utility Wθ t (x i ) is approximated by+,ŵ(x i ) = w (rk(x i ) − 1/2)/λ .
With this approximation, we
˜ θ t (θ t ) ≈ λ −1 ! λ ŵi ∇l˜ (θ t ; x i ).
have an estimate of (5) at θ t as ∇J
i=1Then, we obtain the parameter updateθ t+1 = θ t + η

λ
"
ŵ
i=1

i ˜

λ

∇l (θ t ; x i ) ,

(7)

i=1

implying that we discard all the other samples x ki (k " 1) despitethat they may help to estimate the quantity more accurately. To
utilize all the possible samples, we employ the idea of importancesampling.
!K
Let p̄ be the mixture p̄(x) = (K + 1) −1 k=0
pk (x). One canrewrite the expectation as**
p0 (x)
g(x)p0 (x)dx =g(x)p̄(x)dx(8)p̄(x)
Considering the set of points x ki to be sampled from p̄, we canapproximate the RHS of the above equality by the average

(6)

K "
λ
"
p0 (x ki )1
g(x ki ).
λ(K + 1)p̄(x k )k=0 i=1(9)iwhere η denotes the learning rate and ŵi = ŵ(x i ) for simplicity.
The pure rank-µ update CMA-ES (c1 = 0 in (4)) is considered an
instantiation of the IGO algorithm. Given the multivariate Gaussian
distribution N (m, C) parameterized by θ = (m, C), the parameter
update (6) in the IGO algorithm is equivalent to (2) and (4) if
cm = c µ , c1 = 0, and any function w such that w((i−1/2)/λ) = λwi
for wi defined in (1).

It is an unbiased estimator of (8) shown in [11,15], i.e., the expectedvalue of (9) equals to (8). In this way, we can utilize all the possible
samples without introducing any bias and expect that this leads tomore accurate estimation than the simple Monte-Carlo using only
x 0i . This is the key idea of our work. Note that the case K = 0
recovers (7).2

2.3 Importance Mixing

3.2

An idea of reusing previously generated solutions has been introduced in [13, 14], called importance mixing. It creates the setof points that can be considered as taken from the current distribution pθ t . Provided a population of λ points drawn frompθ t −1 , it accepts each point x in the population with probabilitymin{1, (1 − α)pθ t (x)/pθ t −1 (x)} into new population. Let λ ′ ! λbe the number of accepted points. Then, sample a point x frompθ t and accept it with probability max{α, 1 − pθ t −1 (x)/pθ t (x)} into
the new population and repeat this procedure until λ − λ ′ points
are accepted. Then, the new population is considered as being distributed as pθ t . Since this idea is introduced to reduce the numberof points for which the function value is evaluated at each iteration,
the estimation variance of the natural gradient does not necessarily
lessen. In Section 4, we compare the proposed strategies with theimportance mixing. The parameter α called the minimal refresh
rate is set to α = 0 to reduce the function evaluations as much as
possible in the experiments.

3. SAMPLE REUSE IN THE CMA-ES

In this section, we propose a sample reuse mechanism for the
CMA-ES based on the importance sampling that is a technique toestimate the expectation over a probability distribution by using
samples drawn from a diﬀerent probability distribution. We mix
the current and past samples as drawn from the mixture of thecurrent and past distributions, and estimate the natural gradient by
Monte-Carlo with the samples from the mixture.

3.1 Importance Sampling

Suppose that we have K + 1 probability distributions pk (k =0, 1, . . . , K) and λ independent samples x ki (i = 1, . . . , λ) from each

Sample Reuse in the IGO

We introduce (9) in the IGO framework to estimate the natural˜ θ t (θ t ) using the samples x t−k (i = 1, . . . , λ) from thegradient ∇Jicurrent and!past K distributions, pt−k (k = 0, 1, . . . , K). Let p̄(x) =K pt−k (x).
(K + 1) −1 k=0
fFirst, we need to estimate Wθ t (x it−k ) = w(qθ!t (x it−k )). By
definition, we can write qθ!t (x it−k ) in the form of (8) with g(x) =I{ f (x) ! f (x it−k )}, where I{·} indicates the indicator function.
Then, we can apply the formula (9) to estimate qθ!t (x it−k ) and we
obtainK "λ I{ f (x t−l ) ! f (x t−k )} pt (x t−l )"jij
l=0 j=1

λ(K + 1)p̄(x t−lj )=: q̄(x it−k ) .

(10)

f

Given w, the estimate ŵit−k of the weight value Wθ t (x it−k ) for eachx it−k is computed as()
ŵit−k = w q̄(x it−k ) .
(11)
Next we estimate the natural gradient (5) using (11). The naturalf˜ (θ t ; x).
gradient (5) at θ t is of the form (8) with g(x) = Wθ t (x) ∇l

2 One can create a diﬀerent unbiased estimator that uses all thepossible samples. Instead of considering the mixture p̄, one com!λputes K +1 unbiased estimator λ −1 i=1g(x ki )p0 (x ki )/pk (x ki ) andaverages them, resulting in!K !λλ −1 (K + 1) −1 k=0
g(x ki )p0 (x ki )/pk (x ki ).
i=1

This is also an unbiased estimator of (8). We can prove that (9)has a smaller estimation variance than this estimator, which is also
claimed by [11, 12]. Therefore, we employ (9).

fWith an approximation of Wθ t (x it−k ) for each i and k definedin (11), we can apply the formula (9) and obtain the estimate of˜ θ t (θ t ) as
∇J
K

λ

""
1
˜ (θ t ; x t−k ) ,r̂ (x it−k ) ∇liλ(K + 1)(12)

k=0 i=1

where r̂ (x ki ) = ŵik pθ t (x ki )/ p̄(x ki ). We remark that the proposed
estimator (12) as well as the original estimator (the second termon the RHS of (6)) is biased; its expectation diﬀers from (5) sincef˜ θ t , but they are
the same samples are used to estimate Wθ t and ∇J
consistent; the bias vanishes and both estimates converge to (5)w.p.1 as λ → ∞. See Remark 2 of [1] for the impact of the bias.

3.3 Sample Reuse in the CMA-ES

We introduce the importance sampling to the CMA-ES describedin Section 2.1. More precisely, we replace the mean vector update(2) and the rank-µ update of the covariance matrix in (4) with novel
ones using the samples from the current and past K iteration.
As discussed in Section 2.2, the pure rank-µ update CMA-ES
is considered as an instantiation of the IGO algorithm with themultivariate Gaussian distribution. Therefore, we can simply replace the second term on the RHS of (2) and the rank-µ update of
the covariance matrix in (4) with the novel natural gradient estimate (12). All we need is to choose the function w. We choosew(s) = −2 ln(2s)I{s ! 1/2}, since wi defined in (1) is considered
as an approximation of w((i − 1/2)/λ) for large λ. The weightvalue ŵit−k in (12) is given by (11) with this w. Then, we obtainthe update rules for the mean vector and the covariance matrix asm t+1 = m t + cm
C t+1 = C t + c µK "λ r̂ (x t−k )"iλ(K + 1)k=0 i=1K "λ r̂ (x t−k )"i(x it−k − m t ),λ(K + 1)+ t−k,
· (x i − m t )(x it−k − m t ) T − C t .

(13)
(14)

k=0 i=1

Note that if K = 0, (13) and (14) is equivalent to (2) and (4) with
c1 = 0, except that ŵit /λ is diﬀerent from wi defined in (1). In the
preliminarily experiments we have observed that the sum of ŵit /λtends to be smaller than one when K = 0, whereas the sum of wi isalways one. Its impact is observed in the next section.
We combine the rank-one update to (14). We simply add thesecond term on the RHS of (4) to (14). To update the evolution path,
we solely uses the current samples and use wrk(x t ) as the weightifor each x it rather than ŵit . It means, the evolution path is updated
as it is done in (3). Since the evolution path itself accumulates
the past information, the rank-one update is considered utilizing thepast samples. Therefore, the proposed algorithm exploits the past
information in diﬀerent ways to update the covariance matrix, andwe expect the synergy.
To implement the proposed method, one needs to keep the mean
vectors and the covariance matrices of the past K distributions tocompute the likelihood ratio pθ t (x it−k )/ p̄(x it−k ), which requiresO(K d 2 ) additional memory space. For eﬃcient and numericallystable computation, we keep the log-likelihood of each point, l ik,l =ln pθ t −k (x it−l ), for each i = 1, . . . , λ and k, l = 0, . . . , K, which requires O(λK 2 ) additional space. The likelihood ratio is then com!K
puted as pθ t (x it−k )/ p̄(x it−k ) = (K + 1)( l=0
exp(l il,k − l i0,k )) −1 .
Then the computational complexity at each iteration is O(λK d 2 ).

4.

EXPERIMENTS AND RESULTS

In this section, we evaluate the proposed approach on a standard
benchmark functions. We consider the following four variants of
the CMA-ES using the past solutions:(A) Reuse-m, C The parameters are updated using (13) and (14).
(B) Reuse-C The covariance matrix is updated using (14), while
the mean vector is updated by (2), only current samples areused and the rankings are computed among the current samples as in the standard CMA-ES.
(C) Reuse-m, C + rank-one The rank-one update is incorporated
into Algorithm (A), as is described in Section 3.3.
(D) Reuse-C + rank-one The rank-one update is incorporated into
Algorithm (B).
For these four variants, we observe the eﬀect of K and the population size λ. Moreover, we compare them with the pure rank-µ
update CMA-ES with and without importance mixing and the hybrid of the rank-one and rank-µ update.

4.1

Experimental Procedure

Our benchmark set consists of Sphere, Rosenbrock, Ellipsoid,
Cigar, Ackley, Bohachevsky, Schaﬀer, and Rastrigin functions3.
The global optimum is located at [1, · · · , 1]T for Rosenbrock function and [0, · · · , 0]T for the other functions, where the optimal value
is 0. The problem dimension d of each benchmark function is set
to 20 and 40.
The initial mean vector m0 of the Gaussian distribution is drawn
uniform randomly from [a, b]d for each run, and the covariancematrix is initialized by C 0 = σ 2 I, where σ = (b − a)/2. The
range is set to [1, 5]d for Sphere, Ellipsoid, Cigar, and Rastriginfunctions, [−2, 2]d for Rosenbrock function, [1, 30]d for Ackley
function, [1, 15]d for Bohachevsky function, and [10, 100]d for
Schaﬀer functions, according to [3]. The evolution path pc is
initialized to the zero vector.
Each run is terminated and regarded as success when the best
objective value reaches smaller than 10−10 . On the other hand,
the run is terminated and regarded as failure after d × 106 function evaluations, or when the minimal eigenvalue of the covariancematrix reaches 10−30 , except for Schaﬀer function, where the minimal eigenvalue is set to 10−60 . For each setting we conduct 50
independent runs.
Diﬀerent population sizes (sample size) λ and K are tested. For
Sphere, Rosenbrock, Ellipsoid, and Cigar functions, λ ∈ {4 +⌊3 ln d⌋, 0.2d, 0.4d, 0.8d, d,. 2d, 4d, 8d, 16d}.
For the other multi/
modal functions, λ = 2 2i/2−1 λ base for i = 0, 1, · · · , 7 with
λ base = 2 ln d, d, 2d, 10d for Ackley, Bohachevsky, Schaﬀer, andRastrigin functions, respectively. In preliminary experiments wehave observed the failure with high probability when λ ! λ base .
The number of iterations for sample reuse is K = 0, 1, 3, 5, 7, 9. As
is mentioned in Section 3.3, K = 0 leads to parameter updates
slightly diﬀerent from the standard mean vector update and therank-µ update.
Following [7], we set the learning rates as follows:
4 + µeﬀ /d2
cc =, c1 =,
λ + 4 + 2µeﬀ /λ(d + 1.3) 2 + µeﬀ
01(15)2(µeﬀ − 2 + 1/µeﬀ )c µ = min 1.0 − c1,.
(d + 2) 2 + 2µeﬀ /2

cm = 1,

3 The definition of each benchmark function can be found in [3, 9].

Func. eval. count / success rate

Sphere (20 dim.)

Cigar (20 dim.)

Rosenbrock (20 dim.)

106

Rastrigin (20 dim.)Ackley (20 dim.)

106

105

105

105

105

105
0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

0

2

4

6

8

10

Num. of reuse iterations K
Pure rank-µ CMA(B) Reuse-C

Pure rank-µ CMA + rank1(C) Reuse-m, C + rank1(A) Reuse-m, C(D) Reuse-C + rank1

Figure 1: Average number of function evaluations divided by the success probability.
/are shown on 20 dimensional problems. The default
population size is used for Sphere, Rosenbrock, and Cigar functions and λ = 2 2λ base for Ackley and Rastrigin functions.
Note that for Algorithms (A) and (B), the learning rate for the
rank-one update, c1 , is set to 0 and we do not update the evolution
path.

4.2 Results and Discussion

Figure 1 shows the average number of function evaluations overthe successful runs divided by the success probability on 20 dimensional functions for six algorithms: the pure rank-µ update CMA-ES(c1 = 0 in Section 2.1), the hybrid update CMA-ES (described in
Section 2.1), and Algorithms (A)–(D). The default population size
λ = 4 + ⌊3 ln d⌋. is used/ for Sphere, Rosenbrock, and Cigar functions and λ = 2 2λ base for Ackley and Rastrigin functions. The
result on Ellipsoid function is very similar to that on Cigar function,and the results on Bohachevsky and Schaﬀer functions are similar
to that on Ackley function. Moreover, we have observed a similar
trends on K on the 40 dimensional functions. E.g., with the default
population size, the pure rank-µ update CMA-ES and Algorithm(A) with K = 5 spend 7.5 × 105 and 3.8 × 105 function evaluations
to solve 40 dimensional Sphere function, respectively. Their results
are omitted due to the space limitation.
Figures 2 and 3 show the average number of function evaluations
over the successful runs divided by the success probability on 20
dimensional problems. The performance of the pure rank-µ update
CMA-ES, the pure rank-µ update CMA-ES with importance mixing,
and the pure rank-µ update CMA-ES with sample reuse (Algorithm(A): Reuse-m, C with K = 0, 1, 3, 5) for diﬀerent population sizes
are shown in Figure 2. The results of the hybrid update CMA-ES,
and the CMA-ES with sample reuse (Algorithm (D): Reuse-C with
rank-one with K = 0, 1, 3, 5) for diﬀerent population sizes are shownin Figure 3.
Figure 4 illustrates the transition of the sum of !
the productsof the weight and the likelihood ratio s kt = (1/λ) i r̂ (x ki ) forthe current (k = t) and past samples (k = t − 1, . . . , t − 3) in a
typical run. The results of Algorithm (A) and (B) with K = 3 andthe default population size on 20 dimensional Sphere function are
shown. Smaller values for t − k iteration indicate either that the
samples from t − k iteration are away from the current distribution
and the likelihood ratios are small, or that the samples from t − kiteration are low in the ranking and small or zero weight values are
assigned to them. Remark also that the sum of the coeﬃcients does
not amount to 1, as we have mentioned in Section 3.3.

Figure 5 displays the typical behaviors of Algorithm (A) withthe default population size on 20 dimensional Sphere, Rosenbrock,and Ellipsoid functions and Algorithm (D) with λ = 160 on 20
dimensional Rosenbrock function. The results of both K = 0 and
K = 5 are monitored.

4.2.1

Impact of K

We observe from Figure 1 that the performances of Algorithms(A)–(D) improve more or less monotonically as K increases for all
but Rastrigin function. The algorithms with sample reuse (A and B
or C and D) outperform the one without sample reuse (the pure rankµ update CMA-ES or the hybrid update CMA-ES, respectively)
already with K = 1, where the current and previous populations areused. We see the improvement over K up to 7 in Algorithm (A),whereas the improvement is observed for K up to 3 in Algorithms(B)–(D). The reason that the performance is saturated for large K is
because the likelihood ratios for old samples tend to be small since
past distributions are away from the current distribution and newer
samples tend to have better function values, resulting that very small
weights are assigned to very old samples. On Rastrigin function, alarge K value tends to lead to low success rate. This defect is lessvisible when λ becomes larger as it is seen in Figures 2 and 3.
From Figure 5, the speedup by the sample reuse (K = 5 overK = 0) is seen not only at the convergence stage but also at theadaptation stage.

4.2.2

Impact of λWhen the population size is relatively large, we observe from
Figures 2 and 3 that the sample reuse is not helpful very much to
reduce the number of function evaluations. For large λ, we expectthat the natural gradient estimate in (2) and (4) is relatively accurate
and the eﬀect of the sample reuse is less emphasized. Moreover,since the learning rate c µ defined in (15) becomes larger as λincreases, leading to a big change of the probability distribution,
the likelihood ratios for the past samples will be small. Therefore,
the sample reuse does not help when λ is relatively large. We canobserve the same eﬀect of large λ for Algorithms (B) and (C) as
well.
On Rosenbrock function, the sample reuse turns out to lead to
failure in the case of λ > 102 . As we see in Figure 5, the mean
vector needs to move from the origin towards the global optimum.

10

7

10

6

106

Function evaluation count / success rate

Cigar (20 dim.)

Ellipsoid (20 dim.)

Rosenbrock (20 dim.)Sphere (20 dim.)

106

106

105

105

105
105
101
106

101

102

Bohachevsky (20 dim.)

101

102

Ackley (20 dim.)

101

102

102

Rastrigin (20 dim.)

Schaffer (20 dim.)

107
106
106
105
101

102

105

101

102

105

102

105 2
10

103

Population size λPure rank-µ CMA(A) Reuse-m, C (K = 1)

Pure rank-µ CMA with IM(A) Reuse-m, C (K = 3)

(A) Reuse-m, C (K = 0)
(A) Reuse-m, C (K = 5)

Figure 2: Average number of function evaluations divided by the success probability versus population size λ of the pure rank-µ update
CMA-ES without and with importance mixing (IM), and the rank-µ update CMA-ES with sample reuse (Algorithm (A) for K = 0, 1, 3, 5) on
20 dimensional problems. Missing data implies the failure.
In the case of λ = 160 and K = 5, we observe that the eigenvaluesof the covariance matrix tends to be smaller compared to the caseof K = 0 during the movement of the mean vectors. Smaller
eigenvalues of the covariance matrix leads to producing shorter
steps, x − m, resulting in slow movement of the mean vectors.4

4.2.3 Sample reuse in the pure rank- µ update

Comparing Algorithms (A) and (B) in Figure 1, we find that
Algorithm (B) is better for K = 1 but Algorithm (A) reaches better
performance for larger K. The reason is observed in Figure 4,
where the sums of the coeﬃcients, r̂ (x ki ), which appears in (13)and (14), are more or less equal for k = t, . . . , t − 3 in Algorithm(A), whereas the sums for k = t and t − 1 amount to higher values
than the ones for k = t − 2 and t − 3 in Algorithm (B). This indicatesthat the past samples are used to estimate the natural gradient in
Algorithm (A) as well as the current samples, whereas mostly thecurrent and previous samples are used and the older samples hasless impact in Algorithm (B). In Algorithm (A), the mean vector
is updated by using the proposed update (13), where the natural
gradient is estimated with smaller variance than the one estimatedin the original update (2). Then, the mean vector less fluctuates
and the likelihood ratios are likely to be relatively large for the pastsamples.

4 Empirically, we observe that a larger population size itself leads
to smaller eigenvalues. Since more samples are used in the case of
K = 5, compared to the case of K = 0, the reason of adapting asmaller covariance matrix may be the same reason.

4.2.4

Sample reuse in the hybrid update

Each variant of the algorithms is improved by introducing therank-one update as seen in Figure 1. While the sample reuse for m
update results in better performance when the covariance matrix issolely updated by the rank-µ update, the standard m update leads tofaster convergence when the rank-one update is incorporated. As
discussed in Section 4.2.3, the proposed m update (13) tends tohave less fluctuation in a subspace that are less sensitive in function
value such as the first axis for Cigar function. Then the evolution
path will not be long in this subspace, while it should be in theexample of Cigar function so that the covariance matrix learns theprinciple axis. The standard m update does not disturb the rank-one
update and the speedup is achieved by the sample reuse in the rank-µ
update.
Comparing the performance on Rastrigin function shown in Figures 2 and 3, we notice that Algorithm (D) is less suﬀered from alarge value of K than Algorithm (A). This may be because the pastsamples have less eﬀect in Algorithm (D) for the reason partiallyobserved in Figure 4.
Figure 3 shows that smallest λ leads to the least performance
diﬀerence between diﬀerent K values in Algorithm (D) for theunimodal functions, while the diﬀerence between diﬀerent K values
becomes pronounced for smaller λ in Algorithm (A), shown in
Figure 2. The learning rate c1 for the rank-one update has a relatively
large value for a small population size compared with the learning
rate c µ for the rank-µ update, then the impact of the rank-one update
is dominative.

Sphere (20 dim.)

Function evaluation count / success rate

105

106

Cigar (20 dim.)

Ellipsoid (20 dim.)

Rosenbrock (20 dim.)

105

105

105
101
10

5

102

101

Bohachevsky (20 dim.)

101

102

Ackley (20 dim.)

101

102

102

Rastrigin (20 dim.)

Schaffer (20 dim.)

106

101

102

105

101

102

105

105 2
10

102

103

Population size λPure rank-µ CMA + rank1(D) Reuse-C + rank1 (K = 3)

(D) Reuse-C + rank1 (K = 0)
(D) Reuse-C + rank1 (K = 5)

(D) Reuse-C + rank1 (K = 1)

Figure 3: Average number of function evaluations divided by the success probability versus population size λ of the pure rank-µ update
CMA-ES with rank-one update, and the CMA-ES with sample reuse (Algorithm (D) for K = 0, 1, 3, 5) on 20 dimensional problems.

Sum of coefficients(A) Reuse-m, C
1.0

t−3
t−2
t−1
t

0.8
0.6
0.4
0.2
0.0
0.0

0.2

0.4

0.6
0.8
1.0
1.2
Evaluation count (×105 )

1.4

1.6

Sum of coefficients(B) Reuse-C
1.0

t−3
t−2
t−1
t

0.8
0.6
0.4
0.2
0.0
0.0

0.2

0.4

0.6
0.8
1.0
1.2
Evaluation count (×105 )

1.4

1.6

Figure 4: Transition of the sum of the coeﬃcients (the productsof the weight and the! likelihood ratio) for the current and pastsamples, s kt = 1/λ i r̂ (x ki ), k = t, . . . , t − 3, on 20 dimensional
Sphere function for Algorithm (A) and (B) with K = 3 and thedefault population size. The moving average with window size 11
is computed for smoothing.

4.2.5 Comparison with importance mixing
When introducing the importance mixing in the pure rank-µ update CMA-ES, the performance on the unimodal functions are improved. The improvement is emphasized when λ is relatively large

as is seen in Figure 2. The performance is better than Algorithm(A) for λ " 320. On the other hand, it is likely to be trappedby local minima on the multimodal functions. The motivation of
the importance mixing is to reduce the number of function evaluations by reusing the sample and not to estimate the natural gradientmore precisely. In fact, the estimate of the natural gradient will beless accurate since the number of samples drawn from the current
distribution is less than λ.
Combining the rank-one update and the importance mixing, itceases working properly when the refresh rate α = 0. By tuning
α, the unimodal functions can be solved. However, the α needs tobe tuned dependently of λ. Moreover, it does not work well on themultimodal functions.

5.

CONCLUSION

We have proposed the sample reuse technique based on the importance sampling to improve the accuracy of the estimate of thenatural gradient in the CMA-ES. The quantile based utility function
and then the natural gradient are estimated using the mixture of thecurrent and last K populations. This has been realized thanks to themathematical formulation of the rank-based selection proposed in
the IGO framework. From the experiments, we have found that the
sample reuse is eﬀective both for the mean vector and the covariance matrix update in the pure rank-µ update CMA-ES. On the otherhands, when the rank-one update is employed with rank-µ update(i.e., hybrid update), the sample reuse for the mean vector update
disturbs the evolution path. The best performance is achieved by
the variant that updates the mean vector without the sample reuseand combines the rank-one update and the rank-µ update using thesample reuse, where the reasonable choice for K seems to be K ! 3.

(A) Reuse-m, C (K = 0, λ = 12) on Sphere
105
103
101
10−1
10−3
10−5
10−7
10−9

5

0.5

1.0

1.5

2.0

Mean coordinate val.

0
2.5

−5
0.0

0.5

Evaluation count (×105 )

1.5

2.0

2.5

3.0

100
50
0
−50
0.0

Mean coordinate val.
0.5

2.0

Evaluation count (×105 )

2.5

10
5
0
−5
0.0

1.0

1.5

2.0

2.5

3.0

1.5

2.0

0.5

2.5

Evaluation count (×105 )

3.0

−50
0.0

0.5

1.0

1.5

2.0

Evaluation count (×105 )

1.5

2.0

2.5

3.0

(D) Reuse-C + rank1 (K = 5, λ = 160) on Rosenbrock

Mean coordinate val.

0

1.0

105
103
101
10−1
10−3
10−5
10−7
10−9

Best func. val.
Eigenvalues

50

Mean coordinate val.
1.0

Mean coordinate val.

Evaluation count (×105 )(A) Reuse-m, C (K = 5, λ = 12) on Ellipsoid
105
103
101
10−1
10−3
10−5
10−7
10−9

Best func. val.
Eigenvalues

0.5

2
1
0
−1
0.0

Best func. val.
Eigenvalues

Evaluation count (×105 )(A) Reuse-m, C (K = 5, λ = 12) on Rosenbrock

Mean coordinate val.
1.0

1.5

105
103
101
10−1
10−3
10−5
10−7
10−9

Best func. val.
Eigenvalues

0.5

1.0

(D) Reuse-C + rank1 (K = 0, λ = 160) on Rosenbrock
105
103
101
10−1
10−3
10−5
10−7
10−9

Best func. val.
Eigenvalues

Evaluation count (×105 )(A) Reuse-m, C (K = 5, λ = 12) on Sphere
105
103
101
10−1
10−3
10−5
10−7
10−9

4
2
0
−2
0.0

105
103
101
10−1
10−3
10−5
10−7
10−9

Best func. val.
Eigenvalues

5

Mean coordinate val.

0
−5
0.0

(A) Reuse-m, C (K = 0, λ = 12) on Ellipsoid(A) Reuse-m, C (K = 0, λ = 12) on Rosenbrock
105
103
101
10−1
10−3
10−5
10−7
10−9

Best func. val.
Eigenvalues

2.5

3.0

2
1
0
−1
0.0

Best func. val.
Eigenvalues

Mean coordinate val.
0.5

1.0

1.5

2.0

2.5

3.0

Evaluation count (×105 )Figure 5: The best objective function value, eigenvalues of C, and each coordinate of m with K = 0 (top) and K = 5 (bottom) are shown
for Algorithm (A) with the default population size on 20 dimensional Sphere (1st column), Rosenbrock (2nd column), and Ellipsoid (3rdcolumn) functions and Algorithm (D) with λ = 160 on 20 dimensional Rosenbrock (4th column) function.
Larger K values tend to lead to failure on Rastrigin function. We
have investigated the eﬀect of K on 20 and 40 dimensional problems
and observed a similar trend. More study on the eﬀect of K should
be done on higher dimensional problems.
The step-size adaptation is not introduced in the paper. In the
CMA-ES, the step-size adaptation improves the eﬃciency and therobustness of the algorithm drastically. For instance, it can solve 20
dimensional Ellipsoid function within 2 × 104 function evaluations
with the default population size, whereas the proposed algorithm
needs at least 7 × 104 function evaluations. However, since thisleads to a fast change of the distribution, the likelihood ratio for the
past populations will be very small and the past population cannot
contribute the natural gradient estimate. Therefore, we expect thatthe proposed sample reuse mechanism, as it is, is not eﬀective when
the step-size adaptation is employed. We may need a modification
to eﬀectively reuse the past samples.

6. REFERENCES[1] Y. Akimoto, A. Auger, and N. Hansen. Convergence of theContinuous Time Trajectories of Isotropic Evolution
Strategies on Monotonic C 2 -composite Functions. In
Problem Solving from Nature-PPSN XII, pages 42–51, 2012.
[2] Y. Akimoto, Y. Nagata, I. Ono, and S. Kobayashi.
Bidirectional Relation between CMA Evolution Strategies
and Natural Evolution Strategies. In Parallel Problem
Solving from Nature –PPSN XI, volume 6238 of LNCS, pages
154–163. Springer, 2010.
[3] Y. Akimoto, J. Sakuma, I. Ono, and S. Kobayashi.
Functionally Specialized CMA-ES : A Modification of
CMA-ES Based on the Specialization of the Functions of
Covariance Matrix Adaptation and Step Size Adaptation. In
Proceedings of the Genetic and Evolutionary Computation
Conference (GECCO ’08), pages 479–486, 2008.
[4] S. Amari. Natural Gradient Works Eﬃciently in Learning.
Neural Computation, 10(2):251–276, 1998.
[5] T. Glasmachers, T. Schaul, Y. Sun, D. Wierstra, and
J. Schmidhuber. Exponential Natural Evolution Strategies. In
Proceedings of the Genetic and Evolutionary Computation
Conference (GECCO ’10), pages 393–400, 2010.
[6] N. Hansen. The CMA Evolution Strategy: A Comparing
Review. In J. A. Lozano, P. Larrañaga, I. n. Inza, andE. Bengoetxea, editors, Towards a New Evolutionary

Computation, volume 192 of Studies in Fuzziness and Soft
Computing, pages 75–102. Springer, 2006.
[7] N. Hansen and A. Auger. Principled Design of Continuous
Stochastic Search: From Theory to Practice. In Y. Borensteinand A. Moraglio, editors, Theory and Principled Methods for
Designing Metaheustics, Natural Computing Series, pages
145–180. Springer, 2014.
[8] N. Hansen, A. Auger, R. Ros, S. Finck, and P. Pošík.
Comparing Results of 31 Algorithms from the Black-Box
Optimization Benchmarking BBOB-2009. In Proceedings ofthe Genetic and Evolutionary Computation Conference(GECCO ’10) Companion, pages 1689–1696, 2010.
[9] N. Hansen and A. Ostermeier. Completely Derandomized
Self-Adaptation in Evolution Strategies. Evolutionary
Computation, 9(2):159–195, 2001.
[10] Y. Ollivier, L. Arnold, A. Auger, and N. Hansen.
Information-Geometric Optimization Algorithms: A
Unifying Picture via Invariance Principles. eprint
arXiv:1106.3708v2, 2011.
[11] C. R. Shelton. Policy Improvement for POMDPs Using
Normalized Importance Sampling. Technical Report AI
Momo 2001-002, MIT AI Lab., 2001.
[12] C. R. Shelton. Policy Improvement for POMDPs Using
Normalized Importance Sampling. In Proceedings of theSeventeenth Conference on Uncertainty in Artificial
Intelligence (UAI ’01), pages 496–503, 2001.
[13] Y. Sun, D. Wierstra, T. Schaul, and J. Schmidhuber. Eﬃcient
Natural Evolution Strategies. In Proceedings of the Genetic
and Evolutionary Computation Conference (GECCO ’09),
pages 539–546, 2009.
[14] Y. Sun, D. Wierstra, T. Schaul, and J. Schmidhuber.
Stochastic Search Using the Natural Gradient. In
Proceedings of the 26th International Conference on
Machine Learning (ICML ’09), pages 1161–1168, 2009.
[15] E. Veach and L. J. Guibas. Optimally Combining Sampling
Techniques for Monte Carlo Rendering. In Proceedings ofthe 22nd Annual Conference on Computer Graphics andInteractive Techniques (SIGGRAPH ’95), pages 419–428,
1995.