Revisiting Population Models in Differential
Evolution on a Limited Budget of Evaluations?
Ryoji Tanabe
Yokohama National University, Yokohama, Japanrt.ryoji.tanabe@gmail.com

Abstract. No previous study has reported that differential evolution(DE) is competitive with state-of-the-art black-box optimizers on a limited budget of evaluations (i.e., the expensive optimization scenario).
This is true even for surrogate-assisted DEs. The basic framework of DE
should be reconsidered to improve its performance substantially. In thiscontext, this paper revisits population models in DE on a limited budget of evaluations. This paper analyzes the performance of DE with five
population models on the BBOB function set. Results demonstrate thatthe traditional synchronous model is unsuitable for DE in most cases. In
contrast, the performance of DE can be significantly improved by using
the plus-selection model and the worst improvement model. Results also
demonstrate that DE with a suitable population model is competitive
with covariance matrix adaptation evolution strategy depending on thenumber of evaluations and the dimensionality of a problem.

1

Introduction

Single-objective black-box numerical optimization involves finding a solution x =(x1 , ..., xn )> that minimizes a given objective function f : X → R. Here, X is then-dimensional solution space. Any explicit knowledge of f is unavailable.
This paper considers black-box optimization on a limited budget of evaluations. Optimization with a small number of function evaluations (e.g., 100 × n
evaluations) is generally called expensive optimization. In contrast, this paper denotes optimization with a relatively large number of function evaluations (e.g.,
10 000 × n evaluations) as cheap optimization. Some real-world problems requirea long computation time to evaluate a solution x by expensive computer simulations [3,24] (e.g., CFD [6]). Also, expensive optimization frequently appears in
the filed of hyperparameter optimization of machine learning models [10].
In the evolutionary computation community, surrogate-assisted evolutionary
algorithms are representative approaches for expensive optimization [3,24]. In
general, surrogate-assisted evolutionary algorithms replace an expensive objective function with a cheap surrogate model. Then, the objective value of a new
solution is predicted by the surrogate model based on past solutions found during?

The final authenticated publication is available online at https://doi.org/10.
1007/978-3-030-58112-1_18.

the search process. It is expected that surrogate-assisted evolutionary algorithms
can effectively reduce the number of evaluations by an actual objective function.
Covariance matrix adaptation evolution strategy (CMA-ES) is a variant of
evolution strategies (ES) for numerical optimization [13,19]. Surrogate-assisted
approaches have been intensively studied in the field of ES, especially CMA-ES[4,14,26,33]. Surrogate-assisted CMA-ES has shown state-of-the-art performancefor expensive optimization. Even a non-surrogate-assisted CMA-ES performs
well for expensive optimization. Although Bayesian optimizers (e.g., EGO [25])generally show good performance within a very small number of evaluations, anon-surrogate-assisted CMA-ES performs well after that in most cases [32,36].
Similar to CMA-ES, differential evolution (DE) is a numerical optimizer [42].
The results in the annual IEEE CEC competitions have demonstrated that some
DEs are competitive with more complex optimizers for cheap optimization (e.g.,[46]). However, the number of previous studies on DE algorithms for expensive
optimization is much smaller than that for cheap optimization. In contrast to
the ES community, as pointed out in [7], surrogate-assisted approaches have notreceived much attention in the DE community. Also, a surrogate-assisted DE has
not been competitive even with a non-surrogate-assisted CMA-ES. Most previous studies “avoided” to compare DE with CMA-ES for expensive optimization.
Nevertheless, the simplicity of DE is still attractive in practice. The easy-touse property of DE is valuable for a non-expert user who just uses an evolutionary
algorithm as a black-box optimization tool. In this context, we want to substantially improve the performance of DE for expensive optimization. To achieve thisgoal while keeping the simplicity of DE, this paper revisits population models in
DE for expensive optimization. A population model determines how to update
the population for each iteration. The original DE [42] uses the synchronous
model, which updates all individuals in the population simultaneously by using
one-to-one survivor selection. In addition to the synchronous model, this paper
analyzes the performance of DE with the following four population models: theasynchronous [8,53,54], (µ + λ) [37,49], worst improvement [1], and subset-tosubset [12] models. Although some previous studies (e.g., [8]) investigated theeffect of population models in DE, such an analysis is only for cheap optimization, not for expensive optimization. Thus, little is known about the influence ofthe population model on the performance of DE for expensive optimization.
Our contributions in this paper are at least threefold:
1. We demonstrate that the performance of DE can be significantly improvedby replacing the traditional synchronous model with the (µ + λ) or worst
improvement models. This means that a more efficient surrogate-assisted DE
for expensive optimization could be designed by using the two models.
2. We compare DE with the two models to CMA-ES and other optimizers. We
show that DE with a suitable population model performs better than orsimilar to CMA-ES depending on the number of function evaluations andthe dimensionality of a problem. This is the first study to report such apromising performance of DE for expensive optimization.

Algorithm 1: Synchronous model
Initialize P = {x1 , ..., xµ } randomly;
while The termination criteria are not met do
3
for i ∈ {1, ..., µ} do
4ui ← generateTrialVector(P );
1
2

5
6

for i ∈ {1, ..., µ} do
if f (ui ) ≤ f (xi ) then xi ← ui ;3. DE with the two models can be viewed as a base-line. Any surrogate-assisted
DE should outperform the two non-surrogate-assisted DEs. Some surrogateassisted DEs have been proposed (e.g., [30,34,51,56]). However, benchmarking a surrogate-assisted DE has not been standardized in the DE community.
It is also difficult to accurately reproduce experimental results of most existing surrogate-assisted DEs since their source code is not available throughthe Internet. Consequently, the progress is unclear. Our base-line addresses
this issue, facilitating a constructive development of a surrogate-assisted DE.
The rest of this paper is organized as follows. Section 2 explains the five
population models in DE. Section 3 describes the setting of our computational
experiments. Section 4 shows analysis results. Section 5 concludes this paper.

2

Five population models in DE

Here, we do not consider any method that aims to maintain the diversity in the
population, such as crowding DE [47], island/distributed models [8,52], and thecellular topology [8,35]. These methods aim to prevent the premature convergence of DE to find a good solution with a large number of function evaluations.
Thus, these methods are not suitable for expensive optimization, where DE needs
to quickly find a good solution with a small number of function evaluations.
We carefully surveyed the literature on population models that aim to accelerate the convergence speed of DE. As a result, we found the following four
population models: the asynchronous model [8,53,54], the (µ + λ) model [37,49],
the worst improvement model [1], and the subset-to-subset (STS) model [12].
Algorithms 1–5 show the overall procedure of the basic DE with the traditional synchronous model and the four population models. We extracted only the
population models from their corresponding original algorithms. For example,
the worst improvement model is derived from DE with generalized differential(DEGD) [1], which consists of multiple components. In this study, we want tofocus only on the population model. For this reason, we generalized it so that
we can examine its effectiveness in an isolated manner. For the same reason, we
do not use any surrogate model and any parameter adaptation method for the
scale factor F and the crossover rate C [45]. Below, we explain the five models.
• Synchronous model (Algorithm 1)

The synchronous model is used in the original DE [42] and most recent DE
variants, including jDE [5], JADE [57], CoDE [50], and SHADE [43]. Here, we
explain the synchronous model and some basic operations in DE.
At the beginning of the search, the population P = {x1 , ..., xµ } is initialized(line 1 in Algorithm 1), where µ is the population size. For each i ∈ {1, ..., µ}, xi
is the i-th individual in the population P . Here, xi is an n-dimensional solutionof a problem. For each j ∈ {1, ..., n}, xi,j is the j-th element of xi . According to
the DE terminology, we use the terms “individual” and “vector” synonymously.
After the initialization of P , the following steps (lines 2–6 in Algorithm 1)
are repeatedly performed until a termination condition is satisfied. For eachi ∈ {1, ..., µ}, a trial vector (child) ui is generated (lines 3–4 in Algorithm 1). In
this operation, first, a mutant vector v i is generated by applying a differential
mutation to some individuals in P . Table 1 shows seven representative mutation
strategies in DE. The scale factor F controls the magnitude of the mutation.
Parent indices r1 , r2 , ... are randomly selected from {1, ..., µ} \ {i} such that theydiffer from each other. In Table 1, xbest is the best individual in P . For eachi ∈ {1, ..., µ}, xpbest is randomly selected from the top max(bp µc, 2) individuals
in P , where p ∈ [0, 1] controls the greediness of current-to-pbest/1 [57] and randto-pbest/1 [55]. Also, x̃r2 and x̃r3 in current-to-pbest/1 and rand-to-pbest/1 are
randomly selected from the union of P and an external archive A.
After the mutant vector v i has been generated, a trial vector ui is generated
by applying crossover to xi and v i . In this study, we use binomial crossover [42]:(vi,j if qj ≤ C or j = jrandui,j :=,(1)xi,j otherwise
where qj is randomly selected from [0, 1], and jrand is randomly selected from
{1, ..., n}. The crossover rate C in (1) controls the number of inherited elementsfrom the target vector xi to the trial vector ui .
After the µ trial vectors have been generated, all individuals are updated
simultaneously (lines 5–6 in Algorithm 1). If f (ui ) ≤ f (xi ) for each i ∈ {1, ..., µ},xi is replaced with ui . The individuals that were worse than the trial vectors are
stored in the external archive A. When |A| exceeds a pre-defined size, randomly
selected individuals are deleted to keep the archive size constant.
• Asynchronous model (Algorithm 2)
In contrast to the synchronous model, individuals in the population are updated in an asynchronous manner. Thus, immediately after the trial vector u
has been generated, xi can be replaced with u (lines 4–5 in Algorithm 2). It isdifficult to find out the first study that proposed the asynchronous model in DE
since such an idea is general. In fact, some previous studies did not explicitly
describe that they dealt with the asynchronous model (e.g., [53]). While DE withthe asynchronous model can immediately use a new superior individual for the
search, DE with the synchronous model needs to “wait” by the next iteration.
For this reason, the asynchronous model is generally faster than the synchronous
model in terms of the convergence speed of the population [7,8].
• (µ + λ) model (Algorithm 3)

Table 1: Seven representative mutation strategies for DE.
Strategies

Definitions

rand/1v i := xr1 + F (xr2 − xr3 )rand/2v i := xr1 + Fi (xr2 − xr3 ) + Fi (xr4 − xr5 )best/1v i := xbest + Fi (xr1 − xr2 )best/2v i := xbest + Fi (xr1 − xr2 ) + Fi (xr3 − xr4 )current-to-best/1 v i := xi + Fi (xbest − xi ) + Fi (xr1 − xr2 )current-to-pbest/1 v i := xi + Fi (xpbest − xi ) + Fi (xr1 − x̃r2 )rand-to-pbest/1

v i := xr1 + Fi (xpbest − xr1 ) + Fi (xr2 − x̃r3 )Algorithm 2: Asynchronous model
Initialize P = {x1 , ..., xµ } randomly;
while The termination criteria are not met do
3
for i ∈ {1, ..., µ} do4u ← generateTrialVector(P );5
if f (u) ≤ f (xi ) then xi ← u;
1
2

The elitist (µ + λ) model is general in the field of evolutionary algorithms,
including genetic algorithm and ES. For each iteration, a set of λ trial vectors
Q = {u1 , ..., uλ } are generated (lines 3–5 in Algorithm 3). For each u, the target
vector is randomly selected from the population. Then, the best µ individuals
in P ∪ Q survive to the next iteration (line 6 in Algorithm 3). Unlike other
evolutionary algorithms, the (µ + λ) model has not received much attention in
the DE community. Only a few previous studies (e.g., [37,39,49]) considered the(µ + λ) model. As pointed out in [37], the synchronous model may discard a
trial vector that performs worse than its parent but performs better than otherindividuals in the population. The (µ + λ) model addresses such an issue.
• Worst improvement model (Algorithm 4)
In the worst improvement model [1], only λ worst individuals can generate
trial vectors (lines 3–5 in Algorithm 4). Then, the λ individuals and their λ trial
vectors are compared as in the synchronous model (lines 6–7 in Algorithm 4).
Ali [1] demonstrated that a better individual is rarely replaced with its trial
vector. In other words, the better the individual xi (i ∈ {1, ..., µ}) is, the moredifficult it is to generate ui such that f (ui ) ≤ f (xi ). In contrast, a worse individual is frequently replaced with its trial vector. Based on this observation, Ali
proposed the worst improvement model that allows only the λ worst individuals to generate their λ trial vectors. The number of function evaluations could
possibly be reduced by not generating the remaining µ − λ trial vectors that are
unlikely to outperform their µ − λ parent individuals.
• Subset-to-subset (STS) model (Algorithm 5)

Algorithm 3: (µ + λ) model
Initialize P = {x1 , ..., xµ } randomly;
while The termination criteria are not met do3Q ← ∅;
4for i ∈ {1, ..., λ} do
5u ← generateTrialVector(P ), Q ← Q ∪ {u};1
2

6

P ← Select the µ best individuals from the union P ∪ Q;

Algorithm 4: Worst improvement model
Initialize P = {x1 , ..., xµ } randomly;
while The termination criteria are not met do3
K ← Select indices of the λ worst individuals in P ;4for i ∈ K do
5ui ← generateTrialVector(P );
1
2

6
7

for i ∈ K do
if f (ui ) ≤ f (xi ) then xi ← ui ;The aim of the STS model [12] is the same as that of the (µ + λ) model. The
STS model uses the index-based ring topology. In the STS model, µ individuals
and µ trial vectors are grouped based on their indices and the subset size s ≥ 2.
After µ trial vectors have been generated (lines 4–5 in Algorithm 5), an index i
is randomly selected from {1, ..., µ} (line 6 in Algorithm 5), where i determines
the start position on the index-based ring topology. First, l is set to s or the sizeof the remaining individuals (line 8 in Algorithm 5). Then, a set of l individuals
and l trial vectors are stored into X based on the index-based ring topology(line 9 in Algorithm 5). In Algorithm 5, the function “modulo(a, b)” returns theremainder of the division of a by b. Finally, the l best individuals in X survive
to the next iteration (lines 10–11 in Algorithm 5).
For example, suppose that µ = 5, s = 2, and i = 3. In this case, 5+5 individuals are grouped as follows: {x3 , x4 , u3 , u4 }, {x5 , x1 , u5 , u1 }, and {x2 , u2 }. Whenf (x3 ) = 0.3, f (x4 ) = 0.1, f (u3 ) = 0.4, and f (u4 ) = 0.2, x4 and u4 survive tothe next iteration as follows: x3 := x4 and x4 := u4 . Notice that u4 cannot
survive to the next iteration in the synchronous model since f (u4 ) > f (x4 ).

3Experimental setup

We performed all experiments using the COCO software [16]. The source codeused in our experiments is available at https://github.com/ryojitanabe/de_
expensiveopt. We used the 24 BBOB noiseless functions f1 , ..., f24 [18], which
are grouped into the following five categories: separable functions (f1 , ..., f5 ),functions with low or moderate conditioning (f6 , ..., f9 ), functions with highAlgorithm 5: STS model
Initialize P = {x1 , ..., xµ } randomly;h ← µ/s ;// If modulo(µ, s) 6= 0, h ← µ/s + 1
3 while The termination criteria are not met do
4for i ∈ {1, ..., µ} do
5ui ← generateTrialVector(P );
1
2

6
7
8
9
10
11
12

i ← Randomly select an index from {1, ..., µ};for j ∈ {1, ..., h} dol ← min{µ − ((j − 1) × s), s};K ← {modulo(i, µ), ..., modulo(i + l − 1, µ)},X ← {xk |k ∈ K} ∪ {uk |k ∈ K};
for k ∈ K do
xk ← arg min{f (x)}, X ← X \ {xk };x∈Xi ← modulo(i + s, µ);conditioning and unimodal (f10 , ..., f14 ), multimodal functions with adequate
global structure (f15 , ..., f19 ), and multimodal functions with weak global structure (f20 , ..., f24 ). The dimensionality n of the functions was set to 2, 3, 5, 10, 20,
and 40. For each function, 15 runs were performed. These settings adhere to
the procedure in COCO. According to the expensive optimization scenario in
COCO, the maximum number of function evaluations was set to 100 × n.
In general, the best parameter setting in EAs depends on the allowed budgetof evaluations [9,38,41,44]. Although parameter studies on DE have been well
performed for cheap optimization (e.g., [5,11,58]), little is known about a suitableparameter setting for expensive optimization. Thus, it is unclear how to set
control parameters in DE. We try to address this issue by “hand-tuning” and“automated algorithm configuration”.
We performed “hand-tuning” to select µ, the mutation strategy, and someparameters (i.e., λ and s)Pin the five population models. We mainly used then2Sphere function f (x) =i=1 xi for benchmarking. As a result, we set µ tobαµ ln nc, and αµ = 13. The population was initialized using Latin hypercube
sampling. We used rand-to-pbest/1 in Table 1. As in [57], we set the parametersof rand-to-pbest/1 as follows: p = 0.05 and |A| = µ. As in the standard settingin DE for cheap optimization, we set F and C to 0.5 and 0.9, respectively. Thesesettings were suitable for most population models. In the (µ + λ) and worst
improvement models, we set λ to 1. We also set s in the STS model to 2.
We performed “automatic algorithm configuration” using SMAC [22], which
is a surrogate-model based configurator. We used the latest version of SMAC
(version 2.10.03) downloaded from the authors’ website. We set the cost functionin SMAC (i.e., the estimated performance of a configuration) to the error value|f (xbsf ) − f (x∗ )|, where xbsf is the best-so-far solution found by DE, and x∗ isthe optimal solution of a training problem. We used the 28 CEC2013 functions[28] with n ∈ {2, 5, 10, 20, 40} as the training problems. For each run of DE, theTable 2: Parameters tuned by SMAC, where µ = max{bαµ ln nc, 6}, |A| = bαarc µc,
λ = max{bαλ µc, 1}, and s = max{bαs µc, 2}.
αµStrategyRange[5, 20]

See Table 1

Default

10

rand/1

p

αarc

F

C

αλ

αs

[0, 1] [0, 3] [0, 1] [0, 1] [0, 1] [0, 0.2]
0.5

0

0

Synchronous

9.27 rand-to-pbest/1 0.17 0.58 0.51 0.52--

Asynchronous

9.09 rand-to-pbest/1 0.29 1.81 0.50 0.62--(µ + λ)9.50 rand-to-pbest/1 0.34 1.96 0.53 0.65 0.64-

Worst improvement 7.33 rand-to-pbest/1 0.89 1.62 0.58 0.75 0.22-

STS

0.05 1.0

0.5

5.44 rand-to-pbest/1 0.36 1.95 0.61 0.61

-

0.18

maximum number of function evaluations was set to 100 × n. For each run of
SMAC, the maximum number of configuration evaluations was set to 5 000. For
each population model, five independent SMAC runs were performed. Then, we
evaluated the performance of DE with the five configurations found by SMAC
on the CEC2013 set with n ∈ {2, 5, 10, 20, 40}. Finally, we selected the best one
from the five configurations based on their average rankings by the Friedman
test. Table 2 shows the range of each parameter, the default configuration, andthe best configuration for each model. Interestingly, all configurations include
rand-to-pbest/1. For all configurations, F and C values are relatively similar.

4

Results

This section analyzes the performance of DE with the five population models.
We mainly discuss our results based on the anytime performance of optimizers,
rather than the end-of-the-run results exactly at 100 × n evaluations. For thesake of simplicity, we refer to “a DE with a population model” as “a population
model”. We also use the following abbreviations in Figures 1, 2, and 3: the synchronous (Syn), asynchronous (Asy), (µ + λ) (Plus), worst improvement (WI),and subset-to-subset (STS) models. Section 4.1 demonstrates the performanceof the five population models with the hand-tuned parameters. Section 4.2 examines the effect of automatic algorithm configuration in the population models.
Section 4.3 compares two population models to CMA-ES and other optimizers.
4.1

Comparison of the five population models

Figure 1 shows results of the five population models with the hand-tuned parameters on all 24 BBOB functions with n ∈ {2, 3, 5, 10, 20, 40}. In Figure 1, “best2009” is a virtual algorithm portfolio that consists of the performance data of
31 algorithms participating in the GECCO BBOB 2009 workshop [17].
Figure 1 shows the bootstrapped empirical cumulative distribution (ECDF)
of the number of function evaluations (FEvals) divided by n (FEvals/n) for
31 targets for all 24 BBOB functions. We used the COCO software with the

WI

0.6

Plus

0.4

Asy

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

STS

0.2

0

best 2009

0.8

WI

0.6

Plus

0.4

STS

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.0

2

Asy

0.2

0

log10(# f-evals / dimension)

best 2009

0.8

Plus

0.6

WI

0.4

Asy

0.2

STS

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instancesSynv2.3.2

2

Fraction of function,target pairs

Fraction of function,target pairs

best 2009

0.8

WI

0.6

Plus

0.4

Asy

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.0

2

STS

0.2

0

2
log10(# f-evals / dimension)(c) n = 5

1.0 bbob f1-f24, 20-D

best 2009

0.8

WI

0.6

Plus

0.4

Asy

0.2

STS

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.0

Syn

v2.3.20Synv2.3.2(b) n = 3

1.0 bbob f1-f24, 10-D0

1.0 bbob f1-f24, 5-D

log10(# f-evals / dimension)(a) n = 2

0.0

Syn

v2.3.2

2

Fraction of function,target pairs

0.0

Syn

v2.3.2

1.0 bbob f1-f24, 3-D

Fraction of function,target pairs

best 2009

0.8

Fraction of function,target pairs

Fraction of function,target pairs

1.0 bbob f1-f24, 2-D

1.0 bbob f1-f24, 40-D

best 2009

0.8

WI

0.6

Asy

0.4

Plus

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

STS

0.2
0.0

Syn

v2.3.202log10(# f-evals / dimension)

log10(# f-evals / dimension)

log10(# f-evals / dimension)(d) n = 10(e) n = 20(f) n = 40

Fig. 1: Results of the five population models with the hand-tuned parameters. For eachabbreviation (Syn, Asy, Plus, WI, and STS), see the beginning of Section 4.

“--expensive” option to generate all ECDF figures in this paper. Although
the target error values are usually in {102 , ..., 10−8 }, they are adjusted based on“best 2009”. In ECDF figures, the vertical axis indicates the proportion of target
error values reached by an optimizer within specified function evaluations. For
example, in Figure 1 (e), the synchronous model reaches about 20 percent of
all 31 target error values within 10 × n evaluations on all 24 BBOB functions
with n = 20 in all 15 runs. For more details of ECDF, see [15]. In addition to
ECDF, we analyzed the performance of the population models based on averagerun-time with the rank-sum test (p = 0.05), but the results are consistent with
the ECDF figures in most cases. For this reason, we show only ECDF figures.
Figure 1 shows that the (µ + λ) and worst improvement models show agood performance on any dimensional problems. Although the (µ + λ) model
performs slightly worse than the worst improvement model exactly at 100 × n
evaluations, it performs better than the worst improvement model at the earlystage, especially for n ∈ {20, 40}. The asynchronous model performs better than
the synchronous model. This is consistent with the results in previous studies(e.g., [7,8]). The asynchronous model is also competitive with the (µ + λ) modelfor n = 40. Compared to the STS model, the asynchronous model has theadvantage that it does not require any parameter such as the subset size s.
Overall, we can say that the choice of a population model significantly influences the performance of DE for expensive optimization. The (µ + λ) andworst improvement models are the best when using the hand-tuned parameters.

WI
Plus
T-WI
T-Plus

0.6

Asy
STS

0.4

Syn
T-Syn

0.2

T-STS

0

Plus
STS

0.6

T-Plus
Asy

0.4

Syn
T-STS

0.2

T-Syn

0

1.0 bbob f1-f24, 10-D

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

Plus

0.8

WI
T-STS
T-Plus

0.6

Asy
STS

0.4

T-WI
Syn

0.2

T-Asy
T-Syn

v2.3.2

2

Plus
T-WI
Asy
T-STS
STS

0.4

T-Plus
Syn

0.2

T-Asy

0

2

(c) n = 5best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

T-Plus

0.8

WI
T-STS
Plus

0.6

T-WI
Asy

0.4

T-Asy
STS

0.2

Syn
T-Syn

v2.3.20

T-Syn

v2.3.2log10(# f-evals / dimension)

1.0 bbob f1-f24, 20-D

0.0

WI

0.6

0.0

2

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.8(b) n = 3
Fraction of function,target pairs

Fraction of function,target pairs

1.0 bbob f1-f24, 5-D

log10(# f-evals / dimension)(a) n = 2

0

T-Asy

v2.3.2log10(# f-evals / dimension)

0.0

WI
T-WI

0.0

2

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.8

2

Fraction of function,target pairs

0.0

T-Asy

v2.3.2

1.0 bbob f1-f24, 3-D

Fraction of function,target pairs

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.8

Fraction of function,target pairs

Fraction of function,target pairs

1.0 bbob f1-f24, 2-D

1.0 bbob f1-f24, 40-D

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

T-Plus

0.8

T-STS
T-WI
T-Asy

0.6

WI
T-Syn

0.4

Asy
Plus

0.2

STS
0.0

Syn

v2.3.202log10(# f-evals / dimension)

log10(# f-evals / dimension)

log10(# f-evals / dimension)(d) n = 10(e) n = 20(f) n = 40

Fig. 2: Results of the five population models with automatically-tuned parameters.
The prefix “T-” means that the model uses the automatically-tuned parameters.

In contrast, the traditional synchronous model performs the worst in the five
population models for all dimensions. Based on these observations, we do notrecommend the use of the synchronous model for expensive optimization.

4.2

Effect of automatic algorithm configuration

This section investigates the effect of automatic algorithm configuration in the
population models. Figure 2 shows results of the five population models withthe automatically-tuned parameters. For details of the parameters, see Table 2.
Except for the STS model, the hand-tuned parameters are more suitable
than the automatically-tuned parameters for n ∈ {2, 3, 5, 10}. In contrast, all
the five models with the automatically-tuned parameters outperform those withthe hand-tuned parameters for n = 40. The reason is discussed in Section 4.3.
Although the (µ + λ) model with the automatically-tuned parameters performs
the best exactly at 100×n evaluations for n ∈ {20, 40}, that with the hand-tuned
parameters shows a good performance at the early stage.
In summary, we can obtain almost the same conclusion in Section 4.1 evenwhen using the automatically-tuned parameters. For example, the synchronous
model performs poorly even when using the automatically-tuned parameters.

4.3

Comparison to other optimizers

The results in Section 4.1 and 4.2 show that the worst improvement model withthe hand-tuned parameters (WI) and the (µ + λ) model with the automaticallytuned parameters (T-Plus) perform the best in the five population models forn ≤ 10 and n ≥ 20, respectively. Here, we compare the two best population
models with the following six optimizers. We downloaded the performance dataof the six optimizers from the COCO data archive.
• SMAC-BBOB [23] is a Bayesian optimizer. Although the original version of SMAC
[22] is an algorithm configurator as explained in Section 3, this version of SMAC
(SMAC-BBOB) is a numerical optimizer. SMAC-BBOB is similar to EGO [25].
• lmm-CMA [2] is a surrogate-assisted CMA-ES with local meta-models [4].
• CMAES Hutter [23] is a CMA-ES with the default parameters.
• texp liao [29] is a CMA-ES with automatically-tuned parameters for expensive optimization as in Section 4.2. The irace tool [31] was used for tuning.
• R-SHADE-10e2 [44] is a SHADE [43] with automatically-tuned parameters for
expensive optimization as in Section 4.2. SHADE is one of the state-of-the-art
adaptive DE algorithms. R-SHADE-10e2 uses the synchronous model.
• DE-scipy [48] is DE from the Python SciPy library (https://www.scipy.
org/). DE-scipy can be viewed as a DE used by a non-expert user. DE-scipy
uses the asynchronous model.
Figure 3 shows results of the two population models and the six optimizers.
For n = 40, only the data of texp liao, CMAES Hutter, and DE-scipy wereavailable. Unsurprisingly, SMAC-BBOB and lmm-CMA perform the best at around
10 × n and 100 × n evaluations, respectively. We wanted to know how poorly thetwo population models perform compared to SMAC-BBOB and lmm-CMA.
As seen from Figure 3, the (µ + λ) and worst improvement models perform
significantly better than R-SHADE-10e2 and DE-scipy for any n. This result
provides important information to design an efficient DE, i.e., the basic DE witha suitable population model can outperform even SHADE. Our results indicatethat the default version of DE-scipy is unsuitable for expensive optimization.
The (µ + λ) and worst improvement models perform significantly better than
CMAES Hutter and texp liao for n ∈ {2, 3, 5}. The worst improvement model
performs better than CMAES Hutter and texp liao until 100×n evaluations andperforms similar to CMAES Hutter and texp liao exactly at 100 × n evaluations.
Although the two models perform significantly worse than texp liao for n ∈{20, 40}, they have a better performance than CMAES Hutter at the early stage.
The (µ + λ) model is also competitive with CMAES Hutter at 100 × n evaluations.
The (µ+λ) model in this section and texp liao use the automatically-tuned
parameters for expensive optimization. Interestingly, both optimizers perform
well for high-dimensional functions (n ≥ 20), but they perform poorly for lowdimensional functions (n ≤ 10). This unintuitive observation may come from thedifficulty in finding a parameter set with a good scalability to n. As reported in[40], automatically-tuned parameters are likely to fit for hard-to-solve training
problems. Here, parameters in the (µ + λ) model were tuned on the 28 CEC2013
functions with n ∈ {2, 5, 10, 20, 40}, and parameters in texp liao were tuned on

WI
T-Plus

0.6
0.4
0.2

0

WI
T-Plus0.6

CMAES Hut

CMAES Hut 0.4

texp liao

R-SHADE-1

SMAC-BBOB

DE-scipy-

0.2
0.0

2

R-SHADE-1

0

Fraction of function,target pairs

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

CMAES Hut
WI

0.6

texp liao
0.4

T-Plus

CMAES Hut
0.4

texp liao
R-SHADE-1

0.2

SMAC-BBOB

0.2

SMAC-BBOB
DE-scipy-

v2.3.22

DE-scipy-

v2.3.202(c) n = 5best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

1.0 bbob f1-f24, 40-D

best 2009

0.8

texp liao

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

lmm-CMA-E

0.8

texp liao

CMAES Hut 0.6

0.6

T-Plus
0.4

WI

R-SHADE-1

0

T-Pluslog10(# f-evals / dimension)

1.0 bbob f1-f24, 20-D

lmm-CMA-E

0.8

0.0

WI

(b) n = 3

1.0 bbob f1-f24, 10-D

lmm-CMA-E

0.6

0.0

2

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.8

log10(# f-evals / dimension)(a) n = 2
Fraction of function,target pairs

DE-scipy-

v2.3.2log10(# f-evals / dimension)

1.0 bbob f1-f24, 5-D

lmm-CMA-E

SMAC-BBOB

texp liao
v2.3.2

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

0.8

Fraction of function,target pairs

0.0

1.0 bbob f1-f24, 3-D

lmm-CMA-E

0.8

Fraction of function,target pairs

best 2009

31 targets RLs/dim: 0.5..50
from refalgs/best2009-bbob.tar.gz
15 instances

Fraction of function,target pairs

Fraction of function,target pairs

1.0 bbob f1-f24, 2-D

CMAES Hut
T-Plus

0.4

R-SHADE-1
0.2
0.0

SMAC-BBOB
DE-scipy-

v2.3.2

0

2

WI

0.2
0.0

DE-scipy-

v2.3.202log10(# f-evals / dimension)

log10(# f-evals / dimension)

log10(# f-evals / dimension)(d) n = 10(e) n = 20(f) n = 40

Fig. 3: Results of the two population models (WI and T-Plus) and the six optimizers.

the 19 SOCO functions [20] with n ∈ {5, 10, 20, 40}. It seems that parameters in
the two optimizers fit only for high-dimensional functions. Although addressing
this issue is beyond the scope of this paper, an in-depth analysis is needed.

5

Conclusion

We analyzed the performance of the five population models in DE for expensive
optimization on the 24 BBOB functions. The traditional synchronous model
performs the worst in most cases. In contrast, the worst improvement and (µ +λ) models are suitable for DE on a limited budget of evaluations. The worst
improvement model with the hand-tuned parameters and the (µ+λ) model withthe automatically-tuned parameters perform significantly better than CMA-ES
for n ∈ {2, 3, 5} and are competitive with CMA-ES for n = 10. In summary,our results demonstrated that DE with a suitable population model performs
better than or similar to CMA-ES depending on the number of evaluations anddimensionality of a problem. This is the first study to report such a promising
performance of DE for expensive optimization. DE with the two models should
be a base-line for benchmarking new surrogate-assisted DE algorithms.
We believe that the poor performance of DE for n ∈ {20, 40} compared to
CMA-ES can be addressed by using an efficient parameter adaptation method
for F and C [45]. It is promising to design a surrogate-assisted DE with the worst
improvement and (µ + λ) models. It is also promising to design an algorithm
portfolio [21,27,36] that consists of DE, CMA-ES, and Bayesian optimizers.

Acknowledgments. This work was supported by Leading Initiative for Excellent Young Researchers, MEXT, Japan.

References
1. M. M. Ali. Differential evolution with generalized differentials. J. Comput. Appl.
Math., 235(8):2205–2216, 2011.
2. A. Auger, D. Brockhoff, and N. Hansen. Benchmarking the local metamodel CMAES on the noiseless BBOB’2013 test bed. In GECCO (Companion), pages 1225–
1232, 2013.
3. T. Bartz-Beielstein and M. Zaefferer. Model-based methods for continuous anddiscrete global optimization. Appl. Soft Comput., 55:154–167, 2017.
4. Z. Bouzarkouna, A. Auger, and D. Y. Ding. Local-meta-model CMA-ES for partially separable functions. In GECCO, pages 869–876, 2011.
5. J. Brest, S. Greiner, B. Bošković, M. Mernik, and V. Žumer. Self-Adapting Control
Parameters in Differential Evolution: A Comparative Study on Numerical Benchmark Problems. IEEE TEVC, 10(6):646–657, 2006.
6. S. J. Daniels, A. A. M. Rahat, R. M. Everson, G. R. Tabor, and J. E. Fieldsend. A Suite of Computationally Expensive Shape Optimisation Problems Using
Computational Fluid Dynamics. In PPSN, pages 296–307, 2018.
7. S. Das, S. S. Mullick, and P. N. Suganthan. Recent advances in differential evolution- an updated survey. Swarm and Evol. Comput., 27:1–30, 2016.
8. B. Dorronsoro and P. Bouvry. Improving Classical and Decentralized Differential
Evolution With New Mutation Operator and Population Topologies. IEEE TEVC,
15(1):67–98, 2011.
9. A. S. D. Dymond, A. P. Engelbrecht, and P. S. Heyns. The sensitivity of single
objective optimization algorithm control parameter values under different computational constraints. In IEEE CEC, pages 1412–1419. IEEE, 2011.
10. M. Feurer and F. Hutter. Hyperparameter Optimization. In Automated Machine
Learning - Methods, Systems, Challenges, pages 3–33. Springer, 2019.
11. R. Gämperle, S. D. Müller, and P. Koumoutsakos. A parameter study for differential evolution. In Int. Conf. on Adv. in Intelligent Systems, Fuzzy Systems, Evol.
Comput., pages 293–298, 2002.
12. J. Guo, Z. Li, and S. Yang. Accelerating differential evolution based on a subsetto-subset survivor selection operator. Soft Comput., 23(12):4113–4130, 2019.
13. N. Hansen. The CMA Evolution Strategy: A Tutorial. CoRR, abs/1604.00772,
2016.
14. N. Hansen. A global surrogate assisted CMA-ES. In A. Auger and T. Stützle,
editors, GECCO, pages 664–672, 2019.
15. N. Hansen, A. Auger, D. Brockhoff, D. Tušar, and T. Tušar. COCO: Performance
Assessment. CoRR, abs/1605.03560, 2016.
16. N. Hansen, A. Auger, O. Mersmann, T. Tušar, and D. Brockhoff. COCO: A
Platform for Comparing Continuous Optimizers in a Black-Box Setting. CoRR,
2016.
17. N. Hansen, A. Auger, R. Ros, S. Finck, and P. Posı́k. Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009. In GECCO(Companion), pages 1689–1696, 2010.
18. N. Hansen, S. Finck, R. Ros, and A. Auger. Real-parameter black-box optimization
benchmarking 2009: Noiseless functions definitions. Technical report, INRIA, 2009.

19. N. Hansen and A. Ostermeier. Completely Derandomized Self-Adaptation in Evolution Strategies. Evol. Comput., 9(2):159–195, 2001.
20. F. Herrera, M. Lozano, and D. Molina. Test suite for the spec. iss. of Soft Computing on scalability of evolutionary algorithms and other metaheuristics for large
scale continuous optimization problems. Technical report, Univ. of Granada, 2010.
21. M. D. Hoffman, E. Brochu, and N. de Freitas. Portfolio Allocation for Bayesian
Optimization. In UAI, pages 327–336, 2011.
22. F. Hutter, F. H. Hoos, and K. Leyton-Brown. Sequential Model-Based Optimization for General Algorithm Configuration. In LION, pages 507–523, 2011.
23. F. Hutter, H. H. Hoos, and K. Leyton-Brown. An evaluation of sequential modelbased optimization for expensive blackbox functions. In GECCO (Companion),
pages 1209–1216, 2013.
24. Y. Jin. Surrogate-assisted evolutionary computation: Recent advances and futurechallenges. Swarm and Evol. Comput., 1(2):61–70, 2011.
25. D. R. Jones, M. Schonlau, and W. J. Welch. Efficient global optimization of
expensive black-box functions. J. Global Optimiz., 13(4):455–492, 1998.
26. A. Kayhani and D. V. Arnold. Design of a Surrogate Model Assisted (1 + 1)-ES.
In PPSN, pages 16–28, 2018.
27. P. Kerschke and H. Trautmann. Automated Algorithm Selection on ContinuousBlack-Box Problems by Combining Exploratory Landscape Analysis and Machine
Learning. Evol. Comput., 27(1):99–127, 2019.
28. J. J. Liang, B. Y. Qu, P. N. Suganthan, and A. G. Hernández-Dı́az. Problem
definitions and evaluation criteria for the cec 2013 special session on real-parameter
optimization. Technical report, Nanyang Technological Univ., 2013.
29. T. Liao and T. Stützle. Expensive optimization scenario: IPOP-CMA-ES with apopulation bound mechanism for noiseless function testbed. In GECCO (Companion), pages 1185–1192, 2013.
30. B. Liu, Q. Zhang, and G. G. E. Gielen. A Gaussian Process Surrogate Model
Assisted Evolutionary Algorithm for Medium Scale Expensive Optimization Problems. IEEE TEVC, 18(2):180–192, 2014.
31. M. López-Ibáñez, J. Dubois-Lacoste, L. P. Cáceres, M. Birattari, and T. Stützle.
The irace package: Iterated racing for automatic algorithm configuration. Oper.
Res. Perspect., 3:43–58, 2016.
32. I. Loshchilov and F. Hutter. CMA-ES for Hyperparameter Optimization of Deep
Neural Networks. CoRR, abs/1604.07269, 2016.
33. I. Loshchilov, M. Schoenauer, and M. Sebag. Self-adaptive surrogate-assisted covariance matrix adaptation evolution strategy. In GECCO, pages 321–328, 2012.
34. X. Lu, K. Tang, B. Sendhoff, and X. Yao. A new self-adaptation scheme fordifferential evolution. Neurocomputing, 146:2–16, 2014.
35. N. Lynn, M. Z. Ali, and P. N. Suganthan. Population topologies for particle swarm
optimization and differential evolution. Swarm and Evol. Comput., 39:24–35, 2018.
36. H. Mohammadi, R. L. Riche, and E. Touboul. Making EGO and CMA-ES Complementary for Global Optimization. In LION, 2015.
37. N. Noman and H. Iba. A new generation alternation model for differential evolution. In GECCO, pages 1265–1272, 2006.
38. A. P. Piotrowski. Review of Differential Evolution population size. Swarm andEvol. Comput., 32:1–24, 2017.
39. S. Rahnamayan, H. R. Tizhoosh, and M. M. A. Salama. Opposition-Based Differential Evolution. IEEE TEVC, 12(1):64–79, 2008.
40. S. K. Smit and A. E. Eiben. Parameter Tuning of Evolutionary Algorithms: Generalist vs. Specialist. In EvoApplicatons, volume 6024, pages 542–551, 2010.

41. K. Socha. The Influence of Run-Time Limits on Choosing Ant System Parameters.
In GECCO, pages 49–60, 2003.
42. R. Storn and K. Price. Differential Evolution - A Simple and Efficient Heuristic for
Global Optimization over Continuous Spaces. J. Glo. Opt., 11(4):341–359, 1997.
43. R. Tanabe and A. Fukunaga. Success-History Based Parameter Adaptation for
Differential Evolution. In IEEE CEC, pages 71–78, 2013.
44. R. Tanabe and A. Fukunaga. Tuning differential evolution for cheap, medium, andexpensive computational budgets. In IEEE CEC, pages 2018–2025, 2015.
45. R. Tanabe and A. Fukunaga. Reviewing and Benchmarking Parameter Control
Methods in Differential Evolution. IEEE Trans. Cyber., 50(3):1170–1184, 2020.
46. R. Tanabe and A. S. Fukunaga. Improving the search performance of SHADE
using linear population size reduction. In IEEE CEC, pages 1658–1665, 2014.
47. R. Thomsen. Multimodal optimization using crowding-based differential evolution.
In IEEE CEC, pages 1382–1389, 2004.
48. K. Varelas and M. Dahito. Benchmarking multivariate solvers of SciPy on thenoiseless testbed. In GECCO (Companion), pages 1946–1954, 2019.
49. Y. Wang and Z. Cai. Constrained Evolutionary Optimization by Means of (µ +λ)-Differential Evolution and Improved Adaptive Trade-Off Model. Evol. Comput.,19(2):249–285, 2011.
50. Y. Wang, Z. Cai, and Q. Zhang. Differential Evolution With Composite Trial
Vector Generation Strategies and Control Parameters. IEEE TEVC, 15(1):55–66,
2011.
51. Y. Wang, D. Yin, S. Yang, and G. Sun. Global and Local Surrogate-Assisted
Differential Evolution for Expensive Constrained Optimization Problems With Inequality Constraints. IEEE Trans. Cyber., 49(5):1642–1656, 2019.
52. M. Weber, F. Neri, and V. Tirronen. Distributed differential evolution with
explorative–exploitative population families. GPEM, 10(4):343–371, 2009.
53. M. Wormington, C. Panaccione, K. M. Matney, and D. K. Bowen. Characterizationof structures from X-ray scattering data using genetic algorithms. Phil. Trans. R.
Soc. Lond. A, 357(1761):2827–2848, 1999.
54. M. Zhabitsky and E. Zhabitskaya. Asynchronous Differential Evolution with Adaptive Correlation Matrix. In GECCO, pages 455–462, 2013.
55. J. Zhang and A. Sanderson. Adaptive differential evolution: a robust approach tomultimodal problem optimization, volume 1. Springer, 2009.
56. J. Zhang and A. C. Sanderson. DE-AEC: A differential evolution algorithm basedon adaptive evolution control. In IEEE CEC, pages 3824–3830. IEEE, 2007.
57. J. Zhang and A. C. Sanderson. JADE: Adaptive Differential Evolution With Optional External Archive. IEEE TEVC, 13(5):945–958, 2009.
58. K. Zielinski, P. Weitkemper, R. Laur, and K. D. Kammeyer. Parameter study
for differential evolution using a power allocation problem including interference
cancellation. In IEEE CEC, pages 1857–1864, 2006.