Adaptive Generation-Based Evolution Control for Gaussian Process Surrogate Models

Abstract

The interest in accelerating black-box optimizers has resulted in several surrogate model-assisted version of the Covariance Matrix Adaptation Evolution Strategy, a state-of-the-art continuous black-box optimizer. The version called Surrogate CMA-ES uses Gaussian processes or random forests surrogate models with a generation-based evolution control. This paper presents an adaptive improvement for S-CMA-ES, in which the number of generations using the surrogate model before retraining is adjusted depending on the performance of the last instance of the surrogate. Three algorithms that differ in the measure of the surrogate model's performance are evaluated on the COCO/BBOB framework. The results show a minor improvement on S-CMA-ES with constant model lifelengths, especially when larger lifelengths are considered.

1 Introduction

The problem of optimization of real-valued functions without a known mathematical expression, arising in many engineering tasks, is referred to as continuous black-box optimization. Evolutionary strategies, a class of randomized population-based algorithms inspired by natural evolution, are a popular choice for continuous black-box optimization. Especially the Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [ 6 ] is considered the stateof-the-art continuous black-box optimizer of the several past decades. Since values of a black-box function can only be obtained empirically and at considerable costs in practice, the number of function evaluations needed to obtain a desired function value is a key criterion for evaluating black-box optimizers.

The technique of surrogate modelling aims at saving function evaluations by building a surrogate model of the fitness and using that for a portion of function evaluations conducted in the course of the evolutionary search. Several surrogate model-assisted versions of the CMA-ES have been developed (see [ 11 ] for a recent comparison of some of the most notable algorithms). Surrogate CMA-ES (S-CMA-ES) [ 2 ] utilizes random forests- or Gaussian processes-based surrogate models, which possess an inherent capability to quantify uncertainty of the prediction.

In order to control surrogate model’s error, S-CMA-ES uses the surrogate model for a given number of generations gm before a new instance of the model is trained on a population evaluated with the fitness, which is a strategy called generation-based evolution control [ 9 ]. In [ 2 ], two model lifelengths, in particular gm ∈ {1, 5}, have been benchmarked on the COCO/BBOB framework. In many cases, the higher lifelength outperformed the lower one in earlier phases of the optimization, but the reverse order was observed towards later phases of the optimization. In this paper, we propose an adaptive version of the S-CMA-ES(A-S-CMA-ES), in which an adequate number of model-evaluated generations is estimated as a function of previous model’s error. We restrict our attention to Gaussian processes, since they outperformed random forest-based surrogates [ 2 ]. We consider three different surrogate error measures and compare all of them on the COCO/BBOB framework.

The remainder of this paper is organized as follows. Section 2 outlines basic concepts of S-CMA-ES. The proposed adaptive improvement is described in Section 3. Experimental setup is given in Section 4. Experimental results are reported in Section 6. Section 7 concludes the paper. 2

Surrogate CMA-ES

The CMA-ES operates on a population of λ candidate solutions sampled from a multivariate normal distribution: xk ∼ N (m, σ 2C) k = 1, . . . , λ , (1) where N is the normal distribution function; m and C are the mean and the covariance matrix of the estimated search distribution, respectively; and the σ is the overall search step size. The candidate solutions are ranked according to their fitness values: Upon a (weighted) selection of μ < λ highest ranked points, the mean and the covariance matrix of the multivariate normal distribution are adapted according to a procedure that takes as input, among other variables, a cumulation of the past search steps [ 5 ]. The S-CMA-ES modifies the CMA-ES by replacing its sampling (1) and fitnessevaluation (2) steps with a procedure depicted in Algorithm 1.

4: 5: Algorithm 1 Surrogate part of S-CMA-ES Input: g (generation) gm (number of model-evaluated generations) σ , λ , m, C (CMA-ES internal variables) r (maximal distance between m and a training point) nreq (minimal number of points for training) nmax (maximal number of points for training)

A (archive), fM (model), f (fitness) 1: xk ∼ N (m, σ 2C) k = 1, . . . , λ {sampling} 2: if g is original-fitness-evaluatedthen 3: yk ← f (xk) k = 1, . . . , λ {fitness evaluation} A ← A ∪ {(xk, yk)}kλ=1 (Xtr, ytr) ← choose ntr training points within the Mahalanobis distance r from A , assuring that nreq ≤ ntr ≤ nmax 6: fM ← train_model(Xtr, ytr) 7: mark (g + 1) as model-evaluated 8: else 9: yˆk ← fM (xk) {model evaluation} 10: if gm model generations have passed then 11: mark (g + 1) as original-fitness-evaluated 12: end if 13: end if

λ Output: fM , A , (yk)k=1

Depending on the generation number g, the procedure evaluates all candidate solutions either with the real fitness or with the model. In each case, the sampling of the estimated multivariate normal distribution is unchanged (step 1).

If the population is original-fitness-evaluated (step 3), the new evaluations are saved in an archive of known solutions (step 4). Afterwards, a new model is trained on a set of points within the Mahalanobis distance r from the current CMA-ES distribution N (m, σ C) (step 5).

In model-evaluated generations, the fitness values of the whole population of candidate solutions are estimated by the model (step 9). 2.1

Gaussian Processes A Gaussian process (GP) is a collection of random variables ( f (x))x∈RD , such that any finite subcollection f = ( f (x1), . . . , f (xN )) has an N-dimensional normal distribution. A Gaussian process is defined by a mean function μ(x) (often assumed to be zero) and a covariance function k(x, x; θ ), where θ is a vector of parameters of k, hence hyperparameters of the Gaussian process. Given a set of training data X = {x1, . . . , xN }, the covariance matrix of a GP prior is KN + σn2IN , where KN is a N × N matrix given by {KN }i, j = k(xi, x j; θ ) for all i, j = 1, . . . , N; σn2 is the variance of an additive, i. i. d. noise and IN is a N × N identity matrix. Given a new point x∗ ∈/ X , Gaussian process regression is derived by conditioning the joint normal distribution of ( f (x1), . . . , f (xN ), f (x∗)) on the prior, which yields a univariate Gaussian (see [ 12 ] for more details). The hyperparameters θ of a GP regression model are estimated using the maximum likelihood estimation method. 3

Adaptive Evolution Control for Surrogate CMA-ES

The generation-based evolution strategy optimizes the fitness function and the surrogate model thereof in certain proportion. On problem areas that can be approximated well, a surrogate-assisted optimization might benefit from frequent utilization of the model, while on areas that are hard for the surrogate to approximate, frequent utilization of the model might degenerate the performance due to the model’s inaccuracy.

The proposed evolution control adapts the lifelength of a new model depending on the previous model’s error measured on the current generation.

Consider a generation g that is marked as originalfitness-evaluated, and a newly-trained surrogate model fM . If fM is the first surrogate trained so far, put gm = 1. Otherwise, an error ε of a previous surrogate model f last is estimated on the newly evaluated pop

M ulation (x(1g+1), . . . , x(λg+1)) (Algorithm 2). The error ε is then mapped into a number of consecutive generations gm, gm ∈ [0, gmmax], for which the surrogate fM will be used (Algorithm 3).

We consider three different model error measures. As the CMA-ES depends primarily on the ranking of candidate solutions, the first two proposed error measures, Kendall correlation coefficient and Rank difference are based on ranking. The third one, Kullback-Leibler divergence a. k. a. information gain, measures a difference between a multivariate normal distribution estimated from the fitness valuesy and a multivariate normal distribution estimated for the predicted values yˆ.

Kendall rank correlation coefficient Kendall rank correlation coefficientτ measures similarity between two different orderings of the same set. Let y = ( f (x1) , . . . , f (xλ )) and yˆ= f Mlast(x1), . . . , f Mlast(xλ ) be the sequences of the fitness values and the predicted values of a population x1, . . . , xλ , respectively. A pair of indices (i, j), such that i 6= j, i, j ∈ {1, . . . , λ }, is said to be concordant, if both yi < y j and yˆi < yˆj or if both yi > y j and yˆi > yˆj. A discordant pair (i, j), i 6= j, i, j ∈ {1, . . . , λ } is one fulfilling that both yi < y j and yˆi > yˆj or both yi > y j and yˆi < yˆj. Let nc and dc denote the number of concordant and discordant 9: 10: 11: 12: 13: Algorithm 2 Model error estimation Input: error_type (one of {“Kendall”, “RankDifference”, “Kullback-Leibler”}) g (CMA-ES generation number) x(1g+1), . . . , x(λg+1) (a newly sampled population) y, yˆ(fitness values and model predictions in generation g) ccma = (cc, c1, cμ , cσ , dσ ) (CMA-ES constants) v(cgm)a = (m(g), C(g), p(σg), p(cg), σ (g)) (CMA-ES variables at generation g) εmax (maximal error so far) 1: if error_type = “Kendall” then 2: τ ← Kendall rank correlation coefficient betweeny and yˆ 3: ε ← 12 (1 − τ) 4: else if error_type = “Rank-Difference” then 5: ε ← εRμD(yˆ, y) 6: else if error_type = “Kullback-Leibler” then 7: (m(g+1), C(g+1), σ (g+1)) ← cma_update((x(1g+1), . . . , x(λg+1)), y, ccma, v(cgm)a) (m(Mg+1), C(Mg+1), σ (g+1)) ← cma_update((x(1g+1), . . . , x(λg+1)), yˆ, ccma, v(cgm)a) ε ← DKL(N (m(Mg+1), σ M(g+1)C(Mg+1))k N (m(g+1), σ (g+1)C(g+1))) if ε > εmax then

εmax ← ε end if ε ← εmεax {normalize in proportion to the historical maximum} 14: end if Output: ε ∈ [ 0, 1 ] pairs of indices from {1, . . . , λ }, respectively. The Kendall correlation coefficientτ between vectors y and yˆis defined as: τ =

2 λ (λ − 1) (nc − nd ).

In the corresponding branch of Algorithm 2, the value τ is decreasingly scaled into interval [ 0, 1 ].

Ranking difference error The ranking difference error is a normalized version of a measure used in [ 10 ]. Given r1(i) the rank of the i-th element of yˆand r2(i) the rank of the i-th element of y, the ranking difference error is the sum of element-wise differences between r1 and r2 taking into account only the μ best-ranked points from yˆ: εRμD(yˆ, y) =

∑i:r1(i)≤μ |r2(i) − r1(i)| maxπ ∈Sλ ∑i:π(i)≤μ | i − π(i)| where Sλ is the group of all permutations of set {1, . . . , λ }. Kullback-Leibler divergence Kullback-Leibler divergence from a continuous random variable Q with probability density function q to a continuous random variable P with probability density function p is defined as:

DKL(PkQ) =

p(x) log Z ∞ −∞ p(x) q(x) dx.

For two multivariate normal distributions N1(μ1, Σ1) and N2(μ2, Σ2) with the same dimension k, the KullbackLeibler divergence from N2 to N1 is:

The algorithm of model error estimation (Algorithm 2) in generation g computes Kullback-Leibler divergence from a CMA-estimated multivariate normal distribution N (m(g+1), C(g+1)) w. r. t. fitness values y to a CMA-estimated multivariate normal distribution N (m(Mg+1), C(Mg+1)) w. r. t. predicted values yˆ. Procedure cma_update in steps 7 and 8 refers to one iteration of the CMA-ES from the point when a new population has been sampled. The result is normalized by the historical maximum (step 13).

Adjusting the number of model generations The number of consecutive model generations gm is updated by Algorithm 3. The history of surrogate model errors ε is exponentially smoothed with a rate ru (step 1). The error is truncated at a threshold εT so that resulting gm = gmax for m all values ε ≥ εT (step 3). We consider two different transfer functions T1, T2 : [ 0, 1 ] → [ 0, 1 ] (plotted in Figure 1) that scale the error into the admissible interval [0, gmmax]:

T1(x) = x T2(x; k) =

1 x − 2 Algorithm 3 Updating the number of model generations Input: ε (estimation of surrogate model error, ε ∈ [ 0, 1 ]) εT ∈ [ 0, 1 ] (a threshold at which the error is truncated to 1) γ : [ 0, 1 ] → [ 0, 1 ] (transfer function) ru (error update rate) εlast (model error from the previous iteration) gmax (upper bound for admissible number of model m generations) 1: ε ← (1 − ru)εlast + ruε {exponential smoothing} 2: εlast ← ε

1 min {ε, εT } {truncation to 1} 3: ε ← εT 4: gm ← round(γ(1 − ε)gmmax) {scaling into the admissible interval} Output: gm – updated number of model-evaluated generations Both functions are defined on [ 0, 1 ], moreover, Ti(0) = 0 and Ti(1) = 1 for i = 1, 2. Transfer function T2 is a simple sigmoid function defined to be slightly less sensitive near the edges than in the middle. More control can thus be achieved in the regions of low and high error values. The parameter k determines the steepness of the sigmoid curve. 4

Experimental Setup

The proposed adaptive generation-based evolution control for the S-CMA-ES with three different surrogate model error measures is evaluated on the noiseless testbed of the COCO/BBOB (Comparing Continuous Optimizers / Black-Box Optimization Benchmarking) framework [ 7, 8 ] and compared with the S-CMA-ES and CMA-ES.

Each function is defined everywhere onRD and has its optimum in [ −5, 5 ]D for all dimensionalities D ≥ 2. For every function and every dimensionality, 15 trials of the optimizer are run on independent instances, which differ in linear transformations of the x-space or shifts of the f -space. In our experiments, instances recommended for BBOB 2015 workshop, i. e., {1, . . . , 5, 41, . . . 50}, were used. Each trial is terminated when the fopt is reached within a small tolerance Δ ft = 10−8 or when a given budget of function evaluations, 250D in our case, is used up. Experiments were run for dimensionalities 2, 3, 5, 10 and 20. The algorithms’ settings are summarized in the following subsections. 4.1

CMA-ES The CMA-ES results in BBOB format were downloaded from the BBOB 2010 workshop archive 1. The CMA-ES used in those experiments was in version 3.40.beta and utilized a restart strategy (known as IPOP-CMA-ES), where the population size is increased by factor IncPopSize after 1http://coco.gforge.inria.fr/data-archive/bbob/ 2010/ each restart [ 1 ]. The default parameter values employed in the CMA-ES are λ = 4 + b3logDc, μ = b λ2 c, σstart = 83 , IncPopSize = 2. The S-CMA-ES was tested with two numbers of modelevaluated generations, gm = 1 (further denoted as “GP-1”) and gm = 5 (“GP-5”). All other S-CMA-ES settings were left as described in [ 2 ]. In particular, the Mahalanobis distance was r = 8, the starting values (θ , l) of the Matérn covariance function kMν=at5é/rn2 were (0.5, 2) and the starting value of the GP noise parameter σn2 was 0.01. If not mentioned otherwise, the corresponding settings of adaptive versions of the S-CMA-ES are as just stated.

In order to find the most promising settings for each considered surrogate error measure, a full factorial experiment was conducted on one half of the noiseless testbed, namely on functions fi for i ∈ {2, 3, 6, 8, 12, 13, 15, 17, 18, 21, 23, 24}. The discretization of continuous parameters (γ, εT , gmmax, ru) is reported in Table 1. All possible combinations of the parameters were ranked on the 12 selected functions according to the lowest achieved Δ f med (see Section 6) for different numbers of function evaluations #FEs/D = 25, 50, 125, 250. The best settings were chosen according to the highest sum of 1-st rank counts. Ties were resolved according to the lowest sum of ranks. All of the best settings included maximum model-evaluated generations gmax = 5. The remaining of m the winning values are summarized in the following paragraphs.

Kendall correlation coefficient (ADA-Kendall) Transfer function γ = T2, error threshold εT = 0.5 and update rate ru = 0.2.

Ranking difference error (ADA-RD) The same, except transfer function was γ = T1.

Kullback-Leibler divergence (ADA-KL) Transfer function γ = T2, error threshold εT = 0.9 and update rate ru = 0.5. 5

CPU Timing

In order to assess computational costs other than the number of function evaluations, we calculate CPU timing per function evaluation for each algorithm and each dimensionality. Each experiment was divided into jobs by dimensionalities, functions and instances. All jobs were run in a single thread on the Czech national grid MetaCentrum. The average time per function evaluation for each algorithm and each tested dimensionality is summarized in Table 2. 6

Results

We test the difference in algorithms’ convergence for significance on the whole noiseless testbed with the non-parametric Friedman test [ 3 ]. The algorithms are ranked on each BBOB function with respect to medians of log-scaled minimal distance Δ f from the function optimum, denoted as Δ f med, at a fixed budget of function evaluations.

To account for different optimization scenarios, the test is conducted separately for all considered dimensionalities of the input space and two function evaluation budgets, a higher and a lower one. Let #FEt be the smallest number of function evaluations at which at least one algorithm reached the target, i. e., satisfied Δ f med ≤ Δ ft , or #FEt = 250D if the target has not been reached. We set the higher budget for the tests to #FEt and the lower budget to #F3Et .

Mean ranks from the Friedman test are given in Table 3. The critical value for the Friedman test is 2.29.

The mean ranks differ significantly for all tested scenarios except for both tested numbers of function evaluations in 2D and the higher tested number of function evaluations in 3D. Starting from 5D upwards, the lowest mean rank is achieved either by ADA-Kendall or ADA-RD at both tested #FEs.

In order to show pairwise differences, we perform a pairwise N × N comparison of the algorithms’ average ranks by the post-hoc Friedman test with the BergmannHommel correction of the family-wise error [ 4 ] in cases when the null hypothesis of equal algorithms’ performance was rejected. To better illustrate algorithms differences, we also count the number of benchmark functions at which one algorithm achieved a higher rank than the other. The pairwise score and the statistical significance of the pairwise mean rank differences are reported in Table 4. In the post-hoc test, ADA-Kendall significantly outperforms both the CMA-ES and GP-5 in 10D and 20D. It also significantly outperforms GP-1 in 10D at the higher tested #FEs.

f8 Rosenbrock, original 20D 2

For illustration, the average control frequency given by the ratio of the number of total original-fitness-evaluated generations to the number of total model-evaluated generations within one trial, for data from 15 trials on f8 (Rosenbrock’s function) in 20D is given in Figure 2. The algorithm ADA-KL led to generally lower control frequencies than its competitors, which might explain its slightly inferior performance. Similar results were observed for the remaining functions and dimensionalities.

The cases when ADA-Kendall and ADA-RD are able to switch between more exploitation-oriented and more data-gathering-oriented behaviour can be studied on the results from COCO’s postprocessing. GP-5 outperforms both GP-1 and the CMA-ES on the lower and middle parts of the empirical distribution functions (ECDFs) basically for all dimensionalities (Figure 3). On the other hand, GP1 outperforms GP-5 especially in later phases of the search (Figure 3).

The ability of ADA-Kendall and ADA-RD to switch to a less-exploitation mode when appropriate is eminent on the ECDFs plots in 20D, especially on the moderate and the all-function groups (top right and bottom right on Figure 3), with exception of the well structured multimodal group (middle right), when they fail in the middle part and the weakly structured multimodal group (bottom left), when they fail towards the end of the search. 7

Conclusion

In this paper, we implemented several modifications of the Surrogate CMA-ES (S-CMA-ES), an algorithm using generation-based evolution control in connection with GPs. We considered three measures of surrogate model error according to which an adequate number of upcoming model-evaluated generations could be estimated online. Three resulting algorithms were compared on the COCO/BBOB framework with the S-CMA-ES parametrized by two different numbers of consecutive model-evaluated separable functions 1.0 s bbob - f1-f5, 20-D ira 31 target RLs/dim: 0.5..50 p from refalgs/best2009-bbob.tar.gz te0.8 15 instances g r a t moderate functions 1.0 s bbob - f6-f9, 20-D ira 31 target RLs/dim: 0.5..50 p from refalgs/best2009-bbob.tar.gz te0.8 15 instances g r a t 1 1 1 v2.1 v2.1 v2.1 2 2 2 bbeests2t0209009

AADDAA--KKenenda

AADDAA--RRDD

CCMMAA--EESS

GGPP11

ADA-KL ADA-KL GPP55 bbeests2t0209009

AADDAA--KKenenda

AADDAA--RRDD

GGPP11

CCMMAA--EESS

ADA-KL ADA-KL GPP55 bbeests2t0209009

GGPP11

GGPP55

AADDAA--KKLL

AADDAA--KKenenda

CMA-ES CMA-ES

ADDAA--RRDD n+0.6 o i t c n ff0.4 u o n o i tr0.2 o p o r P 0.0

0 n+0.6 o i t c n ff0.4 u o n o i tr0.2 o p o r P 0.0

0 n+0.6 o i t c n ff0.4 u o n o i tr0.2 o p o r P 0.0 0 n+0.6 o i t c n ff0.4 u o n o i tr0.2 o p o r P 0.0 BBOB 2009 for each selected target. algorithm from BBOB 2009 just not reached them within a given budget of k × DIM, with 31 different values of k chosen equidistant in logscale within the interval {0.5, . . . , 50}. The “best 2009” line corresponds to the best aRT observed during 1 1 1 v2.1 v2.1 v2.1 2 2 2 bbeests2t0209009

AADDAA--RRDD

AADDAA--KKenenda

GGPP11

CCMMAA--EESS

GP5 GP5 ADDAA--KKLL bbeests2t0209009

AADDAA--RRDD

AADDAA--KKenenda

CCMMAA--EESS

AADDAA--KKLL

GP1 GP1 GPP55 bbeests2t0209009

AADDAA--KKenenda

AADDAA--RRDD

GGPP11

CCMMAA--EESS

ADA-KL ADA-KL GPP55 FF 2D generations. Since the work on the adaptive extension is still in progress, the presented results summarize the performance of all compared algorithms on the whole BBOB framework or its function groups. We found two error measures, the Kendall rank correlation and the rank difference error, that significantly outperformed the S-CMA-ES used with a higher number of model-evaluated generations, especially in higher dimensionalities of the input space. However, both of these algorithms provided only a minor improvement of the S-CMA-ES used with a lower number of model-evaluated generations and in some tested scenarios fell behind both tested settings of the S-CMA-ES. An area for further research is the adjustment of other surrogate model parameters beside the control frequency, such as the number of the training points or the radius of the area from which they are selected. 8

Acknowledgments

The research reported in this paper has been supported by the Czech Science Foundation (GACˇ R) grant 17-01251.

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. 2D #FEs⁄#FEt 1 1⁄

3 1⁄

3

GP-1 1⁄ 3 GP-1

References

[1] A. Auger and N. Hansen . A restart CMA evolution strategy with increasing population size . In 2005 IEEE Congress on Evolutionary Computation. IEEE , 2005 .

[2] L. Bajer , Z. Pitra , and M. Holenˇa. Benchmarking Gaussian processes and random forests surrogate models on the BBOB noiseless testbed . In Proceedings of the Companion Publication of the 2015 on Genetic and Evolutionary Computation Conference - GECCO Companion '15 . Association for Computing Machinery (ACM) , 2015 .

[3] J. Demšar . Statistical comparisons of classifiers over multiple data sets . J. Mach. Learn. Res. , 7 : 1 - 30 , December 2006 .

[4] S. García and F. Herrera . An extension on "statistical comparisons of classifiers over multiple data sets" for all pairwise comparisons . J. Mach. Learn. Res. , 9 : 2677 - 2694 , 2008 .

[5] N. Hansen . The CMA evolution strategy: A tutorial . CoRR, abs/1604.00772 , 2016 .

[6] 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 of the 12th annual conference comp on Genetic and evolutionary computation - GECCO '10 , pages 1689 - 1696 , New York, NY, USA, 2010 . ACM.

[7] N. Hansen , S. Finck , R. Ros , and A. Auger . Real-parameter Black-Box Optimization Benchmarking 2009: Noiseless functions definitions . Technical report, INRIA , 2009 , updated 2010 .

[8] N. Hansen , S. Finck , R. Ros , and A. Auger . Real-parameter Black-Box Optimization Benchmarking 2012: Experimental setup . Technical report, INRIA , 2012 .

[9] Y. Jin and B. Sendhoff . Fitness approximation in evolutionary computation-A survey . In Proceedings of the Genetic and Evolutionary Computation Conference, GECCO '02 , pages 1105 - 1112 , San Francisco, CA, USA, 2002 . Morgan Kaufmann Publishers Inc.

[10] S. Kern , N. Hansen , and P. Koumoutsakos . Local Metamodels for Optimization Using Evolution Strategies , pages 939 - 948 . Springer Berlin Heidelberg, Berlin, Heidelberg, 2006 .

[11] Z. Pitra , L. Bajer , J. Repický , and M. Holenˇa. Overview of surrogate-model versions of covariance matrix adaptation evolution strategy . In Proceedings of the Genetic and Evolutionary Computation Conference 2017 , Berlin, Germany, July 15-19 , 2017 (GECCO '17) . ACM, July 2017 . To appear.

[12] C. E. Rassmusen and C. K. I. Williams . Gaussian processes for machine learning . Adaptive computation and machine learning series . MIT Press, 2006 .