Minimizing the number of function evaluations became a very challenging problem in the field of blackbox optimization, when one evaluation of the objective function may be very expensive or time-consuming. Gaussian processes (GPs) are one of the approaches suggested to this end, already nearly 20 years ago, in the area of general global optimization. So far, however, they received only little attention in the area of evolutionary black-box optimization. This work investigates the performance of GPs in the context of black-box continuous optimization, using multimodal functions from the CEC 2013 competition. It shows the performance of two methods based on GPs, Model Guided Sampling Optimization (MGSO) and GPs as a surrogate model for CMA-ES. The paper compares the speed-up of both methods with respect to the number of function evaluations using different settings to CMAES with no surrogate model.
Evolutionary computation became very successful during the past few decades in continuous black-box optimization. In such optimization, we have no mathematical definition of the optimized function available, thus we can calculate analytically neither that function itself, nor its derivatives. In such cases, there is no option but to empirically evaluate the objective function through measurements, tests or simulations.
In various real-world optimization problems, the
evaluation of the objective function is very expensive or
time-consuming, e.g., protein’s folding stability
optimization [
In this paper, we employ two approaches addressing that
task. The first, called Model Guided Sampling
Optimization (MGSO) [
This work investigates the performance of both methods
on the set of niching functions from the CEC 2013
competition [
The following section describes the theoretical fundamentals of GPs and introduces the MGSO method. It also explains surrogate modeling and using GPs as a surrogate model for CMA-ES. Section 3 presents results of experimental evaluation of the considered methods. Section 5 summarizes the results and concludes the paper. 2
GP is a random process such that any finite sequence X1, . . . , Xk of the involved random variables has a multivariate Gaussian distribution. GP is defined by its mean value and covariance matrix described by a function with relatively small number of hyper-parameters, which are usually fitted by the maximum likelihood method. Firstly, GP is trained with N data points from the input space X, D N
XN = {xi|xi ∈ R }i=1 with known input-output values (xN , yN ), then it is used for predicting the (N + 1)-st point. The conditional density of the extended vector yN+1 = (y1, . . . , yN , yN+1), conditioned on XN+1 = XN ∪ {xN+1} is p(yN+1|XN+1) = exp(− 12 yTN+1CN−+11yN+1) p(2π)N+1det(CN+1) where CN+1 is the covariance matrix of a (N + 1)dimensional Gaussian distribution. The covariance matrix can be expressed as
CN+1 =
CN
kT
κk ,
where κ is the variance of the new point itself, k is is the
vector of covariances between the new point and training
data and CN is the covariance of the Gaussian distribution
corresponding to the N training data points [
Covariance functions provide prior information about the objective function and express the covariance between the function values of each two data points xi, x j as cov( f (xi), f (x j)) = k(xi, x j). Because the matrix x1, . . . , xN serves as a covariance matrix, it has to positive semidefinite. 2.1
MGSO has the ability to use a regression model for
prediction and error estimation in order to get the
probability of obtaining a better solution. It was inspired by two
previously proposed methods in the field of black-box
optimization. The first method, Estimation of Distribution
Algorithms [
MGSO was proposed as an alternative method for
Jones’ Efficient Global Optimization (EGO) [
The crucial step in the MGSO algorithm is the sampling of PoI, which is determined by the predicted mean fˆ(x) = yˆand the standard deviation sˆ(x) = sy of the GP model fˆ in any point x of the input space
PoIT (x) = Φ
T − fˆ(x) sˆ(x) = P( fˆ(x ≤ T )), which corresponds to the value of cumulative distribution function of the Gaussian for the value of threshold T . Although all the variables are sampled from Gaussian distribution, PoI(x) is not Gaussian-shaped as it depends on the threshold T and the modelled function f . 2.2
Surrogate modeling is a technique used in optimization in order to decrease the number of expensive function evaluations. A surrogate model, which is a regression model of suitable kind (in our case a GP), is constructed by training with known values of the objective function for some inputs first, and then it is used by the employed evolutionary optimization algorithm instead of the original objective function (in evolutionary optimization usually called fitness) during the search for the global optimum. Although, creating and training a model increases time complexity of the optimization algorithm, using a model instead of the original fitness function decreases the number of its evaluations, which is a crucial objective in expensive optimization.
Every regression model approximates the original
fitness function with some error. To prevent the
optimization from being mislead from such an erroneous
approximation, it is necessary to use the original fitness function
for some subset of evaluations. That subset is determined
by the evolution control (EC) strategy [
An individual-based EC strategy consists in
determining the subset of individuals evaluated by the original
fitness function in each generation. The following
description is illustrated by Figure 1. Denote λ to be the size
of the population provided by CMA-ES. First, λ 00 = α λ
points are sampled from N(m, σ 2C), where α ∈ [
A generation-based EC strategy determines the
number of model-evaluated generations between two
generations evaluated by the original function. After a generation
is evaluated by the original fitness function, the model is
trained using the obtained values. The number of
consecutive model-evaluated generations can be determined also
dynamically, as introduced in so-called adaptive EC
strategy [
Determining the most suitable EC parameters, however, is an open problem, which depends on the properties of the fitness function and the current performance of the surrogate model. Moreover, the most suitable parameters change during the optimization process. 3
Previous investigations compared the performance of
MGSO [
MGSO performance was examined using two covariance functions, isotropic squared exponential (KiSsEo) and squared exponential with automatic relevance determination (KaSrEd), with parameters shown in Table 1. The results in Tables 2 and 3 show the speed-up of MGSO with respect to CMA-ES. As can be seen, the KiSsEo covariance function performed better among these two in more than the half of cases. Table 1 shows used parameter settings in our evaluations. 3.2
The speed-up results are shown in Tables 2 and 3. In performed evaluations, four covariance functions in the GP surrogate model were used, two types of the squared exponential covariance function, the isotropic version KiSsEo(xi, x j) = σ 2f exp
1 − 2`2 (xi − x j)>(xi − x j) , (1) covariance functions starting values
of (σ 2f, `) starting values of σn2 covariance functions starting values
of (σ 2f, `) starting values of σn2
S-CMA-ES
cov ∈ {KνM=at52érn, Kexp, KiSsEo, KaSrEd}
(0.1, 10) for KiSsEo (0.05 × J1,D, 0.1) for KaSrEd (0.5, 2) otherwise
0.01 MGSO cov ∈ {KiSsEo, KaSrEd} (0.1, 10) for KiSsEo (0.05 × (J1,D), 0.1) for KaSrEd 0.01 Kexp, KνM=at52érn, denote, respectively, the isotropic squared exponential, squared exponential with automatic relevance determination, exponential and Matérn with parameter ν = 52 covariance functions. J1,D denotes the vector of ones of length equal to the dimension D of the input space. and the version using automatic relevance determination KaSrEd(xi, x j) = σ 2f exp
1 − 2 (xi − x j)>λ −2(xi − x j) , (2) where λ stands for the characteristic length scale which measures the distance for being uncorrelated along xi. The covariance matrices 1 and 2 differ only when λ is a diagonal matrix instead of a scalar. Two types of the Matérn covariance function were used, (3) ν= 21 KMate´rn(r) = exp
r − ` , which is better known as exponential covariance function (Kexp), and
ν= 52 KMate´rn(r) = σ 2f 1 + √5r ` + 5r2 ! 3`2 exp
√5r !
− `
, (4)
where r = (xi − x j) and the parameter ` is the
characteristic length-scale with which the distance of two considered
data points is compared and σ 2f is the signal variance. The
description of the listed covariance functions can be found
in [
In the performed experiments, different variants of the chosen EC strategies, described in Section 2.2, were examined, generation-based and individual-based. The result are discussed in the following sections. 3.3
Apart from covariance function selection, generationbased EC strategy was determined by two other parameters, the number of model-evaluated generations and the multiplication factor of CMA-ES’ step size σ , which is used in the original-evaluated generations in order to provide points for model training from a broader region of the input space. In the implementation, the first parameter was varied among the values 1, 2, 4 and 8 consecutive modelevaluated generations and the second parameter was varied among the values 1 and 2. 3.4
Apart from covariance function selection, three other
parameters, described in Section 2.2, were examined in the
case of individual-based EC strategy. The first parameter
α ∈ [
In performed evaluations, the parameter α was varied among the values 0, 0.0625, 0.125, and 0.25, the parameter β was varied among the values 5 and 10 and γ was varied among the values 0, 0.1 and 0.2. 4 4.1
“+” mean that, unlike the employed method, no CMAES run was able to reach the target. Signs “*” mean that neither the considered method nor CMA-ES were able to reach the target. Speed-ups written in bold mark cases where the S-CMA-ES’ or MGSO’s median of the ERT is significantly lower than the median of the CMA-ES according to the one-sided Wilcoxon’s test on the significance level α = 0.05. 4.2
The MGSO method brought the highest speed-up in the case of f 2, f 4, 3D version of f 7, 2D version of f 11 and 5D version of f 12. The worst results were observed in the case of f 1, 2D version of f 10 and 3D and 10D versions of f 12. The best results were achieved using KiSsEo covariance function, however, in the case of f 4 and 2D version of f 7 KaSrEd covariance function brought much better results.
In the case of the generation-based EC, the overall best settings with respect to the median values are (KiSsEo, 8, 1) – 8 consecutive model-evaluated generations with unmodified step size in combination with KiSsEo covariance function. The overall best generation-based EC settings showed to be also the best generation-based EC settings of the respective functions, except for 3D version of f 12, where S-CMA-ES performed better using larger step size. Using different covariance functions didn’t bring much better results than the overall best covariance function.
The individual-based EC strategy achieved the best results with the overall settings (KaSrEd, 0, 5, 0.1) – squared exponential covariance matrix with automatic relevance determination, no pre-sampling before training the model, 5 as the multiplicator determining the size of extended population and 0.1 as a multiplicator determining the amount of best points chosen from the extended population. The best results using described parameters were achieved in the case of functions f 3 and 2D and 3D version of f 6. However, the overall performance of the individual-based EC strategy lags far behind the generation-based EC strategy, MGSO and even CMA-ES itself. 4.3
Figure 2 shows examples of the best-fitness progress with the best observed settings (see Table 2 for details). Medians and the first and third quartiles of the best fitness reached are shown; medians and quartiles measured for MGSO and S-CMA-ES on 15 and 10 independent runs (for both EC strategies), respectively.
The optimization progress of the individual-based EC strategy showed to be the slowest in comparison to other methods. MGSO outperformed CMA-ES in most cases and the highest speed-up was achieved in the later phase of the optimization process. The generation-based EC strategy achieved the highest speed-up in the middle phase of the optimization process. However, the generation-based (KaSrEd, 0, 5, 0.1) (KvM=a52térn, 0, 10, 0.2) KiSsEo
KaSrEd param: (KiSsEo, 8, 1) (KaSrEd, 0, 5, 0.1) (KvM=a52térn, 0, 10, 0)
KiSsEo param: (KiSsEo, 8, 1) (KvM=a25térn, 8, 1) (KaSrEd, 0, 5, 0.1) (Kexp, 2−4, 5, 0)
KiSsEo
KaSrEd Table 3: Speed-up of S-CMA-ES using individual- and generation-based strategies and MGSO, compared to CMA-ES without a surrogate model – functions f 11 − f 20 (see Section 4.1 for details). Empty columns signify that the best observed settings for the respective function-dimension combination are identical with the overall best observed settings. EC strategy generally brought the best results, it was outperformed by MGSO in the case of 2D version of f 8 and 5D version of f 12 in the later phase of the optimization process.
In this paper, two optimization approaches based on
Gaussian processes were tested on the set of multimodal fitness
functions from the CEC 2013 competition [
In the case of S-CMA-ES, two evolution control strategies were used, the individual- and generation-based. Although S-CMA-ES using generation-based EC strategy outperformed MGSO, both methods showed the performance improvement in most cases. On the other hand, the individual-based EC strategy brought the worst results of all considered methods. We also observed, that SCMA-ES performs better using generation-based EC setting with more consecutive model-evaluated generations. Isotropic squared exponential covariance function showed to be the most suitable for the optimization from all tested covariance functions.
This is a work in progress that is a part of a broader ongoing research. Terefore, it would be premature to draw deeper conclusions at this stage. We hope to be able to draw such conclusions after further investigations will be performed in the future.
The research reported in this paper has been supported by the Czech Science Foundation (GA CˇR) grant 13-17187S.