Objective function evaluation in continuous optimization tasks is often the operation that dominates the algorithm's cost. In particular in the case of black-box functions, i.e. when no analytical description is available, and the function is evaluated empirically. In such a situation, utilizing information from a surrogate model of the objective function is a well known technique to accelerate the search. In this paper, we review two traditional approaches to surrogate modelling based on Gaussian processes that we have newly reimplemented in MATLAB: Metamodel Assisted Evolution Strategy using probability of improvement and Gaussian Process Optimization Procedure. In the research reported in this paper, both approaches have been for the first time evaluated on BlackBox Optimization Benchmarking framework (BBOB), a comprehensive benchmark for continuous optimizers.
An analytical definition of the objective function in real-world optimization tasks is sometimes hard to obtain. Therefore, neither information about the function’s smoothness nor its derivatives are available. Moreover, evaluation of the function is usually expensive as it can only be done empirically, e.g. by measurement, testing, or running a computer simulation. Such functions are called black-box.
One class of optimization algorithms that are
successfully applied to black-box optimization are
evolution strategies. An evolution strategy is an optimization
method that works on a population of candidate solutions
using evolutionary operators of selection, mutation and
recombination [
Still the whole population must be evaluated which might make running the algorithm for a sufficient number of generations infeasible due to expensive evaluation. To address this issue, techniques that involve a surrogate regression model of the fitness have been proposed.
Two major requirements of incorporating a surrogate model (also called metamodel in the literature) into evolutionary strategies are model management and model selection.
Model management is the task to control the surrogate model’s impact on algorithm’s convergence by using the original fitness alongside its surrogate model in the course of the search.
In evolution control, a certain fraction of individuals or
generations is controlled, i.e. evaluated with the fitness
function, while the remainder is evaluated with the
surrogate model [
For example, Metamodel-Assisted Evolution Strategy
(MAES) uses a surrogate model to pre-select the most
promising individuals before they enter a selection
procedure of a standard ES [
In contrast to evolution control, surrogate approach [
Regarding the model selection, Gaussian processes (GPs) are a non-parameterized regression model that is appealing for the task as it gives its prediction in terms of a Gaussian distribution. The variance of this prediction can be utilized as a confidence measure that promotes exploration of insufficiently modelled areas.
This paper reviews two traditional algorithms
interconnecting Gaussian process-based surrogate models with the
CMA-ES: Metamodel-Assisted Evolution strategy with
improved pre-selection criterion by Ulmer, Strechert and
Zell [
The former is a GP-based MAES with probability of improvement (POI) as a pre-selection criterion.
The latter represents the surrogate approach – in each iteration, a local GP model is built and four functions designed to balance predicted value and variance are optimized.
While both algorithms are GP-based, they differ both in the model-management approach as well as in utilizing GP’s confidence.
The framework COCO/BBOB (Comparing Continuous
Optimizers / Black-Box Optimization Benchmarking) [
Both tested methods had been proposed before the BBOB framework originated. In the research reported in this paper, we evaluated both methods on the noiseless part of the BBOB framework. For that purpose a new implementation1 was required as the original source codes were not available to us.
In the following, we first briefly introduce Gaussian processes as a suitable surrogate fitness model in Section 2. An exposition of the tested methods is given in Section 3 and experimental setup in Section 4. Section 5 presents experimental results and finally Section 6 concludes the paper. 2
Both algorithms under review feature the choice of Gaussian processes as a surrogate fitness function model.
GPs are a probabilistic model with several properties that make it well suited for fitness function modelling: its hyperparameters are comprehensible and limited in number and it provides a confidence measure given by standard deviation of predicted value at new data points.
In the following, we define a GP using notation and
equations from Büche [
Consider f : RD → R an unknown real-parameter function to be approximated. GP model is specified by a set XN = xi | xi ∈ RD iN=1 of N training data points with known function values tN = {ti | f (xi) = ti}iN=1. The data are assumed to be a sample of zero-mean multivariate Gaussian distribution with joint probability density p(tN | XN ) = exp(− 12 tTN CN−1tN ) p(2π)N det(CN ) (1) where the covariance matrix CN is defined by means of a covariance function C(xi, x j, Θ) i, j ∈ {1, . . . , N} with a fixed set of hyperparameters Θ.
1The sources are available at ❤tt♣s✿✴✴❣✐t❤✉❜✳❝♦♠✴r❡♣❥❛❦✴ s✉rr♦❣❛t❡✲❝♠❛❡s✴ C (xp, xq) = θ exp
For a set tN+1 that includes a new observation tN+1 = f (xN+1), we obtain
Using Bayesian rule for conditional probabilities, the prediction at a new data point has the density function p(tN+1 | XN+1, tN ) = p(tN+1 | XN+1) .
p(tN | XN ) The covariance matrix CN+1 can be written with the use of the covariance matrix CN as
CN+1 =
CN kT k κ where k = (C(xi, xN+1))1N is a vector of covariances between the new point and XN and κ = C(xN+1, xN+1) is the new point’s variance.
Using (1) and (2) together with the fact that the inverse
CN−+11 can also be expressed by the means of CN−1, (3) can
be simplified to a univariate Gaussian [
p(tN+1 | XN+1, tN ) ∝ exp with mean and variance given by 1 (tN+1 − tˆN+1)2 ! − 2 σt2N+1 (5) tˆN+1 = kT CN−1tN , σt2N+1 = κ − kT CN−1k.
A comprehensive exposition of Gaussian processes can
be found in [
The covariance function plays an important role as it expresses prior assumptions about the shape of the modelled function. The vector Θ of model’s hyperparameters is usually fitted with the maximum likelihood method.
One of the most commonly used covariance functions is squared exponential covariance function:
C (xp, xq) = θ exp
1 (xp − xq)T (xp − xq) − 2 r2 where the parameter θ scales the covariance between two points and radius r is the characteristic length scale of the process.
When two points are close to each other compared to the characteristic length scale r, the covariance is close to one while it exponentially decays to zero with their distance growing.
The squared exponential covariance function can be improved with automatic relevance determination (ARD): where the radii ri scale the impact of two points distance on their correlation separately in each dimension.
To address the need to balance the exploitation of a fitted model prediction and the exploration of regions where model’s confidence is low, the standard deviation of GP’s prediction can be utilized.
One possibility are merit functions proposed in [
Another option is the probability of improvement (POI). Let us assume Gaussian process prediction to be a random variable Y (x) with mean tˆ(x) and standard deviation σ (x). For a chosen threshold T less than or equal to the so-far best obtained fitness value fmin, the probability of improvement in point x is defined as:
POIT (x) = p(Y (x) ≤ T ) = Φ T − tˆ(x) σ (x) (9) where Φ is the cumulative distribution function of the distribution N (0, 1).
3.1
GP Model Assisted Evolution Control A standard ES operates on λ offsprings generated from μ parents by evolutionary operators of recombination and mutation. After the fitness of the offsprings is evaluated, a population of μ best individuals is selected to reproduce to a new offspring in the next iteration. In (μ, λ ) ES, μ best individuals are selected from the λ offsprings, whereas in (μ + λ ) ES, μ best individuals are selected from the union of the offprings and their parents.
MAES [
The GP model is trained in every generation on a set of Ntr most recently evaluated individuals (step 8).
In this paper we consider two pre-selection criteria in
accordance with [
Algorithm 1 GP Model-Assisted Evolution Strategy Input: f – fitness function μ – number of parents λ – population size λPre – size of pre-selected population Ntr – size of training dataset χM – the preselection criterion that depends on a GP model M, e.g. mean model prediction or POI 1: Pop ← generate and evaluate λ initial samples 2: M ← a GP model trained on points from Pop 3: while termination criteria not reached do 4: Offspring ← reproduce Pop into λPre new points 5: Offspring ← mutate Offspring 6: Offspring ← select best λ points according to the pre-selection criterion χM 7: evaluate Offspring with f 8: M ← update model M on Ntr points most recently evaluated with f 9: Pop ← select μ points best according to f 10: end while 3.2
Gaussian Process Optimization Procedure
Gaussian Process Optimization Procedure is due to Büche,
Schraudolph and P. Koumoutsakos [
CMA-ES is used as the evolution strategy, with the number of parents μ set to 2 and the population size λ set to 10.
The pseudocode is given in Figure 2. After generating an initial population, a local GP model is built and utilized in an iterative procedure. Considering possibly low accuracy and computational infeasibility of global models, the training dataset is restricted to NC points closest to the current best known solution xbest (step 6) and NR most recently evaluated points (step 7).
If the GP model M has been successfully trained then the CMA-ES optimizes four merit functions fM,α (8) for each α ∈ {0, 1, 2, 4}. Areas that might be approximated inaccurately are avoided by bounding the ES to the hypercube spanning the set of NC points selected from xbest’s neighborhood (steps 10 and 12). The points that are optima of the considered merit functions are evaluated by the original fitness and added to the dataset of known points.
In the case that no new solution is found, a random perturbation (step 23) is evaluated and added to the dataset. Unfortunately, authors don’t specify the parameter m that occurs in 23. We set it to the value m = 1.
The authors used the following covariance function in
their GP model:
C (xp, xq) = θ1exp
where ri, θ1, θ2, θ3 > 0 and δpq is the Kronecker delta.
The function is the sum of the squared exponential
covariance function with ARD (7), constant shift θ2 and a white
noise scaled with θ3.
Benchmarking) [
Each function is defined on [
Since source codes were available for neither of the tested methods, we implemented them in MATLAB.
For Gaussian processes, we chose MATLAB’s default implementation, ❢✐tr❣♣, a part of Statistics and Machine Learning Toolbox. The GP hyperparameters fitting was done with a quasi-newton optimizer, which is the default in ❢✐tr❣♣.
Parameters of the benchmarked algorithms were set as
follows:
CMA-ES. A multi-start version with the population size
doubled after each restart was used in the tests (MATLAB
code v. 3.62.beta). Number of restarts was set to 4,
while other parameters settings were left default: λ =
4 + ⌊3log10D⌋, σstart = 38 , IncPopSize = ⌊λ /2⌋.
MAES. We implemented GP Model Assisted Evolution
Strategy on top of a framework developed for S-CMA-ES
algorithm, which employs a GP model in conjunction with
a generation evolution control [
The number of parents and population size were set to correspond with the CMA-ES settings.
Two pre-selection criteria were tested: MMP and the POI with threshold equal to the so-far best sampled fitness value fmin. In both cases λPre was set to 3λ . The training set was comprised from 2λ most recently evaluated points.
We used the same covariance function as in the GPOP case, i.e. (10).
GPOP. We adhered to [
The termination criteria were chosen as follows: • number of consecutive iterations with no new solution found is larger than 2 while the tolerance on the two points euclidean distance for them to be considered equal is 10−8 23: xR ← 24: yR ← f (x) 25: A ← A ∪ (xR, yR) 26: end if 27: end while 28: return xbest 4
The framework COCO/BBOB (Comparing Continuous Optimizers / Black-Box Optimization • the overall change of fitness values during the last 10 iterations is lower than 10−9 • the target fitness value is reached
Training set size parameters NC and NR were chosen in
accordance with empirical results reported in [
Although the performance is measured in terms of the number of fitness evaluations, other operations such as optimizing a surrogate model may also be costly in benchmarking scenarios. If we consider functions in 10 D, a run of GPOP on one core of a computational grid took approximately 27.8 real hours per function on average. For those reasons, we limited the maximum number of fitness evaluations to 100 ∗ D for all tested methods. 5
Results from experiments on all the 24 noiseless BBOB benchmark functions are presented in Figures 1–3.
The expected running time (ERT) depicted in Figure 1
depends on a given target value, ft = fopt + Δ f , i.e. the
true optimum fopt of the respective benchmark function
raised by a small value Δ f . The ERT is computed over
all relevant trials as the number of the original function
evaluations (FE) executed during each trial until the
target function value ft reached, summed over all trials and
divided by the number of successful trials, i.e. trials that
actually reached ft: [
ERT( ft) = #FE( fbest ≥ ft) #succ (11)
In the graphs for functions 1, 2, 5, 8, 9, 12, 14, GPOP achieved significantly better results compared to all other algorithms for some dimensions. In contrast, the differences between MAES and CMA-ES are rather small, regardless of the pre-selection criterion.
The graphs in Figure 2 summarize the performance over subgroups of the benchmark functions for the highest tested dimension 10. The graphs show the proportion of algorithm runs that reached a target value ft ∈ 10[−1..2] (see the figure caption for further details).
The GPOP speedup is most eminent on the group of
separable functions (functions 1–5), that is functions,
optimization of which can be reduced to D one-dimensional
problems [
Similar results may be observed on Figure 3 that for each function shows the dependence of the relative best fitness value on the number of evaluations in 10 D. As can be seen on the graphs for functions 13 and 24, GPOP in some cases outperforms all other algorithms in early stages, but then gets trapped in a local minimum.
On the graph for functions 21 on Figure 3, MAES with MMP as the pre-selection criterion visibly outperforms all other algorithms. 6
In this paper, we compared the CMA-ES and our implementation of two traditional methods which improve upon the CMA-ES with Gaussian process-based surrogate models: MAES, which extends the CMA-ES with a preselection step, and GPOP, which iteratively optimizes the GP model.
The benchmarks on the BBOB framework did not show any significant speedup of MAES compared to the CMA-ES. On the other hand, GPOP in many cases outperforms all the other methods, especially in early optimization stages. On some functions, though, it tends to get trapped in local minima. This might be explained with the fact that GPOP requires considerably fewer function evaluations per iteration than other methods. However, the model is built locally and might not be sufficient for exploration of the search space in later phases.
Our MAES implementation relies on a modification of
the sampling step in CMA-ES, thereby changing the
distribution of the sampled population. This might mislead
CMA-ES a bit and requires further research, in particular
considering the proposal in [
Another area for further research is the exploration of various confidence measures across tested methods, especially in connection with GPOP. 7
Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme "Projects of Large Research, Development, and Innovations Infrastructures" (CESNET LM2015042), is greatly appreciated.
The research performed by Lukáš Bajer has been supported by the Czech Health Research Council project NV15-33250A.
The research performed by Jakub Repický was supported by SVV project number 260 224. 9 Rosenbrock rotated -1
0 1 2 log10 of (# f-evals / dimension) 5 Linear slope