Combining Gaussian Processes and Neural Networks in Surrogate Modeling for Covariance Matrix Adaptation Evolution Strategy

Abstract

This paper focuses on surrogate models for Covariance Matrix Adaptation Evolution Strategy (CMA-ES) in continuous black-box optimization. Surrogate modeling has proven to be able to decrease the number of evaluations of the objective function, which is an important requirement in some real-world applications where the evaluation can be costly or time-demanding. Surrogate models achieve this by providing an approximation instead of the evaluation of the true objective function. One of the stateof-the-art models for this task is the Gaussian process. We present an approach to combining Gaussian processes with artificial neural networks, which was previously successfully applied to other machine learning domains. The experimental part employs data recorded from previous CMA-ES runs, allowing us to assess different settings of surrogate models without running the whole CMA-ES algorithm. The data were collected using 24 noiseless benchmark functions of the platform for comparing continuous optimizers COCO in 5 different dimensions. Overall, we used data samples from over 2.8 million generations of CMA-ES runs. The results examine and statistically compare six covariance functions of Gaussian processes with the neural network extension. So far, the combined model did not show up to outperform the Gaussian process alone. Therefore, in conclusion, we discuss possible reasons for this and ideas for future research.

Keywords

• black-box optimization
• surrogate modeling
• artificial neural networks
• Gaussian processes
• covariance functions

-

In evolutionary black-box optimization, an increasing attention is paid to tasks with expensive evaluation of the black-box objective function. It is immaterial whether that expensiveness is due to time-consuming computation like in long simulations [ 17 ], or due to evaluation through costly experiments like in some areas of science [ 3 ]. To deal with such expensive evaluations, black-box optimization has in the late 1990s and early 2000s adopted an approach called surrogate modeling or metamodeling [ 7, 13, 15, 34, 42, 45, 48 ].

Basically, a surrogate model in continuous black-box optimization is any regression model that with a sufficient fidelity approximates the true black-box objective function and replaces it in some of its evaluations. Among the regression models most frequently used to this end, we are most interested in Gaussian processes (GPs) [ 47 ], due to the fact that they estimate not only the expected value of the true objective function, but the whole distribution of its values. Apart from regression, they are also encountered in classification, and play a key role in Bayesian optimization [ 22, 28 ], where the distribution of values provides selection criteria that are alternatives to the objective function value, such as the probability of improvement or the expected improvement.

The importance of GPs in machine learning incited attempts to integrate them with the leading learning paradigm of the last decades – neural learning, including deep learning. The attractiveness of this research direction is further supported by recent theoretical results concerning relationships of asymptotic properties of important kinds of artificial neural networks (ANNs) to properties of GPs [ 35,39,41 ]. The integration of GP with neural learning has been proposed on two different levels: (i) Proper integration of an ANN with a GP, in which the

GP forms the final output layer of the ANN [ 9, 52 ]. (ii) Only a transfer of the layered structure, which is a crucial feature, to the GP context, leading to the concept of deep GPs (DGPs) [ 8, 12, 22, 23 ].

The recalled research into the integration of GPs with neural learning has used regression data [ 8, 9, 12, 23, 52 ], mostly from the UCI Machine Learning Repository [ 50 ], but also data concerning locomotion of walking bipedal robots [ 9 ], and face patches and hadwritten digits [ 52 ]. In [ 12,23 ], also classification data were used. However, we are aware of only one application of such an integration, in particular of DGPs, to two very easy 1- and 2-dimensional optimization problems [ 22 ]. Hence, there is a gap between the importance of GPs in Bayesian optimization and missing investigations of the suitability of integrated GP-ANN models for surrogate modeling in black-box optimization. That gap motivated the research reported in this paper.

We have performed this novel kind of investigation of GP-ANN integration using data from more than 5000 runs of using GPs as Bayesian optimizers in black-box optimization, i.e. optimization of pointwise observable functions for which no analytic expression is available and the function values have to be obtained either empirically, e.g. through measuring or experiments, or through numerical simulations. Based on that data, our investigation addressed two research questions: 1. Does the integration of a GP with neural learning has an added value compared with employing a GP alone or an ANN alone? 2. What is the impact of each of the considered GP covariance functions on the result of GP-ANN integration?

The next section introduces the principles of surrogate modeling as well as of the evolutionary optimization method Covariance Matrix Adaptation Evolution Strategy (CMA-ES), in the context of which our research has been performed. In Section 3, the fundamentals of GPs are recalled and the method we have used for their integration with neural networks is explained. The core part of the paper is Section 4, which attempts to provide experimental answers to the above research questions. 2

Surrogate modeling in Black-Box

Optimization

Basically, the purpose of surrogate modeling – to approximate an unknown functional dependence – coincides with the purpose of response surface modeling in the design of experiments [ 24, 40 ]. Therefore, it is not surprising that typical response surface models, i.e., low order polynomials, belong also to the most traditional and most successful surrogate models [ 1, 2, 20, 29, 45 ]. Other frequently used kinds of surrogate models are artificial neural networks of the kind multilayer perceptron (MLP) or radial basis function network [ 18, 25–27, 42, 53 ], and Gaussian processes, in surrogate modeling also known as kriging [ 6, 7, 13–15, 33, 34, 37, 49, 51 ]. Occasinally encountered are support vector regression [ 10, 36 ] and random forests [ 5, 44 ]. 2.1

CMA-ES and Its Surrogate-Assisted Variant DTS-CMA-ES

The CMA-ES algorithm performs unconstrained optimization on Rd , by means of iterative sampling of populations sized l from a d-dimensional Gaussian distribution N(m; s 2C), and uses a given parent number m among the sampled points corresponding to the lowest objective function values, to update the parameters of that distribution. Hence, it updates the expected value m, which is used as the current point estimate of the function optimum, the matrix C and the step-size s . The CMA-ES is invariant with respect to strictly increasing transformations of the objective function. Hence, to make use of the evaluations of the objective function in a set of points, it needs to know only the ordering of those evaluations. Details of the algorithm can be found in [ 19, 21 ].

During the more than 20 years of CMA-ES existence, a number of surrogate-assisted variants of this algorithm have been proposed, a survey can be found in [ 6, 43 ]. Here, we pay attention only to the most recent GP-based among them, the Doubly Trained Surrogate CMA-ES (DTS-CMA-ES) [ 6 ], surrogate-assisted variant of CMAES. It employs two GPs f1 and f2, trained consecutively, to find an evaluation of the population x1; : : : ; xl , with f1 used for active learning of training data for f2. Due to the CMA-ES invariance with respect to strictly increasing transformations, it evaluates the difference between predictions only according to the difference in the ordering of those predictions, using the ranking difference error (RDE). The RDE of y 2 Rl with respect to y0 2 Rl considering k best components is defined:

RDE(y; y0) = k åi;(r(y0))i k(jr(y0))i (r(y))ij maxp2P(l ) åik=1 ji p 1(i)j ; (1) where P(l ) denotes the set of permutaions of f1; : : : ; l g and r(y) denotes the ordering of the components of y, i.e., (8y 2 Rl ) r(y) 2 P(l ) & (r(y))i < (r(y)) j ) yi y j. Because the CMA-ES algorithm uses the surrogate model to select the most promising candidates for true evaluation, the metric considers only k best samples. The range of RDE metric is [0; 1], it equals 0 for the exact ordering and 1 for thereverse order.

The algorithm DTS-CMA-ES is described in Algorithm 1, using the following notation: • A for an archive – a set of points that have already been evaluated by the true black-box objective function BB; • ds2C for the Mahalanobis distance given by s 2C: ds2C(x; x0) = q (x x0)>s 2C 1(x x0); (2) • Nk(x; A) for the set of k of ds2C-nearest neighbours of x 2 Rd with respect to the archive A; • fi(x1; : : : ; xl ) = ( fi(x1); : : : ; fi(xl )), for i = 1; 2; • Th for the training set selection:

l Th = [ fx 2 Nh(x j; A)jds2C(x; x j) < rmaxg j=1 (3) with rmax > 0 for h = 1; : : : ; jAj; • k(A) = maxfhjjThj

Nmaxg, with Nmax 2 N; • rPoI for decreasing ordering of f1(x1); : : : ; f1(xl ) according to the probability of improvement with respect to the lowest BB value found so far, i < j ) PoI(rPoI( f1(x1; : : : ; xl )))i;V )

PoI(rPoI( f1(x1; : : : ; xl ))) j;V ); (4) where V = minx2A BB(x). Algorithm 1 Algorithm DTS-CMA-ES Require: x1; : : : ; xl 2 Rd , m, A, s and C – step size and matrix from the CMA-ES distribution, Nmax 2 N such that Nmax l , rmax > 0, b ; emin; emax; einitial; amin; amax; ainitial 2 (0; 1) 1: if this is the 1st call of the algorithm in the current CMA-ES run then a = ainitial; e = einitial take over the returned values of a; e from its previous call in the run 5: end if 6: Train a Gaussian process f1 on Tk(A), estimating mGP, sn; s f ; ` through maximization of the likelihood (7) 7: Evaluate BB(x j) for x j not yet BB-evaluated and such that (rPoI( f1(x1; : : : ; xl ))) j dal e 8: Update A to A [ f(x jj(rPoI( f1(x1); : : : ; f1(xl ))) j dal eg 9: Train a Gaussian process f2 on Tk(A), estimating mGP, sn, s f , ` through maximization of the likelihood (7) 10: For x j such that (rPoI( f1(x1; : : : ; xl ))) j dal e, update f2(x j) = BB(x j) 11: Update e to b RDEm ( f1(x1; : : : ; xl ); ( f2(x1; : : : ; xl )) + (1 b )e and a to amin + max(0; min(1; emeax emeimnin )) 12: For j fulfilling (rPoI( f1(x1; : : : ; xl )) j > dal e, update f2(x j) to f2(x j) minf f2(x j0)j(rPoI( f1(x1; : : : ; xl)) j0 > dal eg + minf f2(x j0 )j(rPoI( f1(x1; : : : ; xl )) j0 dal eg 13: Return f2(x1); : : : ; f2(xl ); e; a 3

Gaussian Processes and Their

Integration with Neural Networks

3.1

Gaussian processes

A Gaussian process on a set X Rd ; d 2 N is a collection of random variables ( f (x))x2X , any finite number of which has a joint Gaussian distribution [ 47 ]. It is completely specified by a mean function m : X ! R, typically assumed constant, and by a covariance function k : X X ! R such that for x; x0 2 X ,

E f (x) = m; cov( f (x); f (x0)) = k(x; x0): (5) Therefore, a GP is often denoted GP(m; k(x; x0)) or GP(m; k).

The value of f (x) is typically accessible only as a noisy observation y = f (x) + e, where e is a zero-mean Gausssian noise with a variance sn > 0. Then cov(y; y0) = k(x; x0) + sn2I(x = x0); (6) where I(proposition) equals for a true proposition 1, for a false proposition 0.

Consider now the prediction of the random variable f (x?) in a point x? 2 X if we already know the observations y1; : : : ; yn in points x1; : : : ; xn. Introduce the vectors x = (x1; : : : ; xn)>, y = (y1; : : : ; yn)> = ( f (x1) + e; : : : f (xn) + e)>, k? = (k(x1; x?); : : : ; k(xn; x?))> and the matrix K 2 Rn n such that (K)i; j = k(xi; x j) + sn2I(i = j). Then the probability density of the vector y of observations is p(y; m; k; sn2) = exp 12 (y m)>K 1(y

m) p(2p)2 det(K + sn2In) ; (7) where det(M) denotes the determinant of a matrix M. Furthermore, as a consequence of the assumption of Gaussian joint distribution, also the conditional distribution of f (x?) conditioned on y is Gaussian, namely

N(m(x?) + k?K 1(y m); k(x?; x?) k?>K 1k?): (8)

According to (6), the relationship between the observations y and y0 is determined by the covariance function k. In the reported research, we have considered 6 kinds of covariance functions, listed below. In their definitions, the notation r = kx0 xk is used, and among the parameters of k, aka hyperparameters of the GP, frequently encountered are s 2f; ` > 0, called signal variance and characteristic length scale, respectively. Other parameters are introduced for each covariance function separately. (i) Linear: kLIN(x; x0) = s02 + x>x0, with a bias s02. (ii) Quadratic is the square of the linear covariance: kQ(x; x0) = (s02 + x>x0)2. (iii) Rational quadratic: kRQ(x; x0) = s 2f 1 + 2ar2`2 with a > 0. (iv) Squared exponential: kSE(x; x0) = s 2f exp (v) Matérn 52 :

kM52 atérn(x; x0) = s 2f 1 + p`5r + 35`r22 exp p`5r . (vi) One composite covariance function, namely the sum of kSE and kQ: kSE+Q(x; x0) = kSE(x; x0) + kQ(x; x0). a 3.2

GP as the Output Layer of a Neural Network

The approach integrating a GP into an ANN as its output layer has been independently proposed in [ 9 ] and [ 52 ]. It relies on the following two assumptions:

1. If nI denotes the number of the ANN input neurons, then the ANN computes a mapping net of nI -dimensional input values into the set X on which is the GP. Consequently, the number of neurons in the last hidden layer equals the dimension d, and the ANN maps an input v into a point x = net(v) 2 X , corresponding to an observation f (x + e) governed by GP (Figure 1). From the point of view of the ANN inputs, the GP is now GP(m(net(v)); k(net(v); net(v0))).

2. The GP mean m is assumed to be a known constant, thus not contributing to the GP hyperparameters and independent of net.

Due to the assumption 2., the only hyperparameters of the GP are the parameters q k of the covariance function. outputGPlayer:y

(+distributionofy) lasthiddenlayer:x . .

.

1sthiddenlayer inputlayer:inputvaluesv

I

I

As to the ANN, it depends on the one hand on the vector qW of its weights and biases, on the other hand on the parameters q A of its architecture (such as the number of hidden layers, the numbers of neurons in each of them, the activation functions). Altogether, the integrated GP-ANN model depends on the vector (q k ; qW ; q A), which we in accordance with [ 9, 52 ] and with the terminology common for GPs also call hyperparameters, although in the context of ANNs, this term is typically used only for q A.

Consider now n inputs to the neural network, v1; : : : ; vn, mapped to the inputs x1 = net(v1); : : : ; xn = net(vn) of the GP, corresponding to the observations y = (y1; : : : ; yn)>. Then the log-likelihood of q is where m is the constant assumed in Assumption 2., and

L(q ) = ln p(y; m; k; sn2); (K)i; j = k(net(vi); net(v j)): (9) (10) To find the combination of values of q with the maximal likelihood, methods for smooth optimization can be applied only with respect to q k and qW , with respect to which is L continuous. With respect to q A, we have two possibilities. First, to fix q A in advance, thus effectively restricting the hyperparameters of the integrated model to (q k ; qW ). Second, to perform the optimization with respect to q A in an outer loop (using, e.g. grid search combined with cross-validation) and for each considered combination of values of q A perform the optimization with respect to (q k ; qW ) in an inner loop.

Finally, let the smooth optimization be performed in the most simple but, in the context of neural networks, also most frequent way – as gradient descent. The partial derivatives forming Ñ(qk ;qW )L can be computed as: ¶ L n ¶ L

i; j=1 ¶¶KLi; j ¶¶Kqi`k;j ; ¶ q`W = ¶ q`k = å

n = å i; j;k=1 ¶ Ki; j ¶ xk ¶ L ¶ Ki; j ¶ net(vk) : ¶ q`W (11) In (11), the partial derivatives ¶¶KLi; j , i; j = 1; : : : ; n, are components of the matrix derivative ¶¶KL , for which the calculations of matrix differential calculus [ 38 ] together with (7) and (9) yield ¶ L ¶ K = As a dataset to compare different configurations of the combined GP-ANN model, we used recorded DTS-CMAES runs from previous experiments. This allowed us to effectively evaluate the surrogate model on its own without having to perform the whole optimization. The opensource Matlab implementation of DTS-CMA-ES, used to obtain the data, is available at [ 4 ]. The underlying surrogate models were Gaussian process with 8 different covariance functions implemented using the GPML Toolbox [ 46 ]. The optimization runs were collected on the platform for comparing continuous optimizers COCO. The employed black-box optimization benchmark set contains 24 noiseless functions scalable to any input dimension. We used dimensions 2, 3, 5, 10, and 20 in 5 different instances, which consist in random rotation and translation in the input space. Therefore, for each combination of benchmark function and dimension, 40 runs are available. More important than the number of runs, however, is the number of their generations because data from each generation apart from the first can be used for testing all those surrogate models that could be trained with data from the previous generations. The number of generations in a particular run of CMA-ES algorithm is unknown before the run and depends on the objective function and its particular instance. The benchmarks and their total numbers of available generations for individual dimensions are listed in Table 1. We skipped the linear function (benchmark function number 5), which is easy to optimize, and the recorded runs did not provide enough data samples to train and evaluate the surrogate models. 4.2

Experimental Setup

We compare the combinations of neural networks with Gaussian processes using six different covariance functions listed in Subsection 3.1. All models were trained on the same set of training data Tk(A) as was used in steps 6 and 9 of Algorithm 1. Due to the condition (3), this set is rather restricted and allows training only a restricted ANN to prevent overfitting. Therefore, we decided to use a multilayer perceptron with a single hidden layer, thus a topology (nI ; nH ; nO), where nI is the dimension of the training data, i.e. nI 2 f2; 3; 5; 10; 20g, and nH = nO = 82 if nI = 2; > <

3 if nI = 3; 5; >:5 if nI = 10; 20: (13) As the activation function for both the hidden and output layer, we chose the logistic sigmoid. We trained the weights and biases of the neural network together with the parameters of the Gaussian process as proposed in [ 52 ] and outlined in Subsection 3.2. As a loss function, we used the Gaussian log-likelihood and optimized the parameters with Adam [ 30 ] for a maximum of 1000 iterations. We also kept a 10 % validation set out of the training data to monitor overfitting, and we selected the model with the lowest L2 validation error during the training. 4.3

Results

To evaluate different covariance functions, we used the recorded data in a similar way as it would be used by DTSCMA-ES algorithm. We took each generation of points except the first one and we used it as testing data. Every point evaluated by the true objective function in all previous generations is then available as a training sample. We filter these samples using the training set selection method Tk(A) and trained the surrogate model on it.

This way, we evaluated six different models on every generation of samples listed in Table 1. The results are presented first in a function-method view in Table 4, then in a dimension-method view in Table 5. The metric used in both tables is RDE, which shows how precise the model is in ordering of the predicted values, therefore the lower is the error value, the better. The first table contains results for every benchmark function separately, averaged over every input dimension, function instance and DTSCMA-ES generation as well as the average over the whole group of functions. The second table provides a view on how the dimension of the input space affects the approximation error. The results for a particular dimension are averaged over all functions in a specific group, instances, and generations. We also visualized with boxplots the summarized RDE values for the function groups in Figure 2.

Moreover, we verified the results by performing multiple comparison tests. A Non-parametric Friedman test was conducted on RDEs across all results for particular functions and function types in Table 4, and for a particular combination of dimensions and function types in Table 5. If the null hypothesis of the equality of all six considered methods was declined, the Wilcoxon signed-rank test was performed for all pairs of covariance functions, and its results were corrected for multiple hypotheses testing using the Holm method. We summarized the results of statistical testing in Tables 2 and 3. The results show that the combined model performs best with the simplest linear kernel. We compared the GPANN with a linear kernel with pure Gaussian processes in paper [ 32 ]. We found out that if we compare the combined GP-ANN with GP both with linear kernels, the neural extensions can bring better results in some cases. However, using more complex covariance functions with GP is still better. Therefore, we tried to apply it also to the GP-ANN combination. Unfortunately, the results with other kernels ended up much worse.

In our future research, we would like to try to overcome this. We want to systematically investigate different ANN topologies, including the direction of deep Gaussian processes [ 8, 12, 22, 23 ], in which only the topology is used from an ANN, but all neurons are replaced by GPs. Moreover, in addition to the selection of the training set used in DTS-CMA-ES and described in Algorithm 1, we want to consider also alternative ways of training set selection, allowing to train larger networks. Finally, we intend to perform research into transfer learning of surrogate models: An ANN-GP model with a deep neural network will be trained on data from many optimization runs, such as those employed in this paper, and then the model used in a new run of the same optimizer will be obtained through additional fine tuning restricted only to the GP and last 1-2 layers of the ANN.

We would also like to try to change the way how the model is trained. In paper [ 52 ] the parameters of Gaussian process are learned together with the parameters of the neural network. We think that this might not work well in the domain of surrogate modeling in black-box optimization. Therefore, we will try to train the network and Gaussian process separately. 5

Conclusion

In this paper, we examined an extension of Gaussian processes used in surrogate modeling. At the beginning, we described the CMA-ES algorithm for black-box optimization and its surrogate-assisted variant DTS-CMA-ES. Then we outlined Gaussian processes, various covariance functions, and their neural network extension. The presented research is our first attempt in the application of the combination of Gaussian processes with neural networks as a surrogate model in black-box optimization. We implemented the combined model and compared the results obtained with six different covariance functions. By looking at the results presented in the previous section, we can conclude that the combined model yields the best results with the simplest linear kernel. The other covariance functions produce much larger errors and the worst performing one seems to be the composite kernel kSE+Q. Unfortunatelly the results showed up to be worse than using Gaussian processes alone, therefore we discuss possible ideas of further improvements at the end of the paper.

Acknowledgement

The research reported in this paper has been supported by SVV project number 260 575 and was also partially supported by Czech Science Foundation (GA CˇR) grant 18-18080S. Computational resources were supplied by the project "e-Infrastruktura CZ" (e-INFRA LM2018140) provided within the program Projects of Large Research, Development and Innovations Infrastructures. n d u

o fu e r

c g e a r c p a n s l ia t cu r u i a p tr v n a o i x io un is t f

c t r n n o u e f f r s e

k f eu ra i

f l d a

m r v h e R b d

e e f g g o a a r r e e e p v v y a a t

e f e r o h a n t s E c v

n o D e p c e t h in n t e s r d n e n o f f a it i d n c Moderately Ill-Conditioned κ LIN κ Q κ RQ κ SE κ Mat κ SE+Q Highly Ill-Conditioned κ LIN κ Q κ RQ κ SE κ Mat κ SE+Q

Multi-modal with Global Structure

Multi-modal Weakly Structured

References

[1] A. Auger , D. Brockhoff , and N. Hansen . Benchmarking the local metamodel cma-es on the noiseless BBOB'2013 test bed . In GECCO'13 , pages 1225 - 1232 , 2013 .

[2] A. Auger , M. Schoenauer , and N. Vanhaecke . LS-CMAES: A second-order algorithm for covariance matrix adaptation . In Parallel Problem Solving from Nature - PPSN VIII , pages 182 - 191 , 2004 .

[3] M. Baerns and M. Holenˇa. Combinatorial Development of Solid Catalytic Materials . Design of High-Throughput Experiments, Data Analysis, Data Mining . Imperial College Press / World Scientific, London, 2009 .

[4] L. Bajer and Z. Pitra. Surrogate CMA-ES. https:// github.com/bajeluk/surrogate-cmaes, 2021 .

[5] L. Bajer , Z. Pitra , and M. Holenˇa. Benchmarking Gaussian processes and random forests surrogate models on the BBOB noiseless testbed . In GECCO'15 Companion , pages 1143 - 1150 , 2015 .

[6] L. Bajer , Z. Pitra , J. Repický , and M. Holenˇa . Gaussian process surrogate models for the CMA evolution strategy . Evolutionary Computation , 27 : 665 - 697 , 2019 .

[7] A.J. Booker , J. Dennis, P.D. Frank , D.B. Serafini , Torczon V., and M. Trosset . A rigorous framework for optimization by surrogates . Structural and Multidisciplinary Optimization , 17 : 1 - 13 , 1999 .

[8] T. Bui , D. Hernandez-Lobato , J. Hernandez-Lobato , Y. Li , and R. Turner . Deep Gaussian processes for regression using approximate expectation propagation . In ICML , pages 1472 - 1481 , 2016 .

[9] R. Calandra , J. Peters , C.E. Rasmussen , and M.P. Deisenroth . Manifold Gaussian processes for regression . In IJCNN , pages 3338 - 3345 , 2016 .

[10] S.M. Clarke , J.H. Griebisch , and T.W. Simpson . Analysis of support vector regression for approximation of complex engineering analyses . Journal of Mechanical Design , 127 : 1077 - 1087 , 2005 .

[11] The COCO platform , 2016 . http://coco.gforge.inria.fr.

[12] K. Cutajar , E.V. Bonilla , P. Michiardi , and M. Filippone . Random feature expansions for deep Gaussian processes . In ICML, pages 884 - 893 , 2017 .

[13] M.A. El-Beltagy , P.B. Nair , and A.J. Keane . Metamodeling techniques for evolutionary optimization of computaitonally expensive problems: Promises and limitations . In Proceedings of the Genetic and Evolutionary Computation Conference , pages 196 - 203 . Morgan Kaufmann Publishers, 1999 .

[14] M. Emmerich , K. Giannakoglou , and B. Naujoks . Singleand multi-objective evolutionary optimization assisted by Gaussian random field metamodels, . IEEE Transactions on Evolutionary Computation , 10 : 421 - 439 , 2006 .

[15] M. Emmerich , A. Giotis , M. Özdemir , T. Bäck , and K. Giannakoglou . Metamodel-assisted evolution strategies . In PPSN VII , pages 361 - 370 . ACM, 2002 .

[16] S. Finck , N. Hansen , R. Ros , and A. Auger . Real-parameter black-box optimization benchmarking 2010: Presentation of the noisy functions . Technical report, INRIA , Paris Saclay, 2010 .

[17] A. Forrester , A. Sobester , and A. Keane . Engineering Design via Surrogate Modelling: A Practical Guide . John Wiley and Sons, New York, 2008 .

[18] H.M. Gutmann . A radial basis function method for global optimization . Journal of Global Optimization , 19 : 201 - 227 , 2001 .

[19] N. Hansen . The CMA evolution strategy: A comparing review . In Towards a New Evolutionary Computation , pages 75 - 102 . Springer, 2006 .

[20] N. Hansen . A global surrogate assisted CMA-ES . In GECCO'19 , pages 664 - 672 , 2019 .

[21] N. Hansen and A. Ostermaier . Completely derandomized self-adaptation in evolution strategies . Evolutionary Computation , 9 : 159 - 195 , 2001 .

[22] A. Hebbal , L. Brevault , M. Balesdent , E.G. Talbi , and N. Melab . Efficient global optimization using deep Gaussian processes . In IEEE CEC , pages 1 - 12 , 2018 .

[23] G. Hernández-Muñoz , C. Villacampa-Calvo , and D. Hernández Lobato . Deep Gaussian processes using expectation propagation and Monte Carlo methods . In ECML PKDD , pages 1 - 17 , paper no. 128 , 2020 .

[24] S. Hosder , L. Watson , and B. Grossman . Polynomial response surface approximations for the multidisciplinary design optimization of a high speed civil transport . Optimization and Engineering , 2 : 431 - 452 , 2001 .

[25] Y. Jin , M. Hüsken , M. Olhofer , and Sendhoff B. Neural networks for fitness approximation in evolutionary optimization . In Y. Jin, editor, Knowledge Incorporation in Evolutionary Computation , pages 281 - 306 . Springer, 2005 .

[26] Y. Jin , M. Olhofer , and B. Sendhoff . Managing approximate models in evolutionary aerodynamic design optimization . In CEC 2001 , pages 592 - 599 , 2001 .

[27] Y. Jin , M. Olhofer , and B. Sendhoff . A framework for evolutionary optimization with approximate fitness functions . IEEE Transactions on Evolutionary Computation , 6 : 481 - 494 , 2002 .

[28] D.R. Jones , M. Schonlau , and W.J. Welch . Efficient global optimization of expensive black-box functions . Journal of Global Optimization , 13 : 455 - 492 , 1998 .

[29] S. Kern , N. Hansen , and P. Koumoutsakos . Local metamodels for optimization using evolution strategies . In PPSN IX , pages 939 - 948 , 2006 .

[30] D.P. Kingma and J. Ba . Adam: A method for stochastic optimization . Preprint arXiv:1412.6980 , 2014 .

[31] J. Koza and J. Tumpach . Surrogate networks . https:// github.com/c0zzy/surrogate-networks, 2021 .

[32] J. Koza , J. Tumpach , Z. Pitra , and M. Holenˇa . Using past experience for configuration of Gaussian processes in black-box optimization . In 15th Learning and Intelligent Optimization Conference, page accepted for publication , 2021 .

[33] J.W. Kruisselbrink , M.T.M. Emmerich , A.H. Deutz , and T. Bäck . A robust optimization approach using kriging metamodels for robustness approximation in the CMA-ES . In IEEE CEC , pages 1 - 8 , 2010 .

[34] S.J. Leary , A. Bhaskar , and A.J. Keane . A derivative based surrogate model for approximating and optimizing the output of an expensive computer simulation . Journal of Global Optimization , 30 : 39 - 58 , 2004 .

[35] J. Lee , Y. Bahri , R. Novak , S.S. Schoenholz , J. Pennington , et al. Deep neural networks as Gaussian processes . In ICLR , pages 1 - 17 , 2018 .

[36] I. Loshchilov , M. Schoenauer , and M. Sebag . Intensive surrogate model exploitation in self-adaptive surrogateassisted CMA-ES (saACM-ES) . In GECCO'13 , pages 439 - 446 , 2013 .

[37] J. Lu , B. Li , and Y. Jin . An evolution strategy assisted by an ensemble of local Gaussian process models . In GECCO'13 , pages 447 - 454 , 2013 .

[38] J.R. Magnus and H. Neudecker. Matrix Differential Calculus with Applications in Statistics and Econometrics . John Wiley and Sons, Chichester, 2007 .

[39] A.G.G. Matthews , J. Hron , M. Rowland , and R.E. Turner . Gaussian process behaviour in wide deep neural networks . In ICLR , pages 1 - 15 , 2019 .

[40] R.H. Myers , D.C. Montgomery , and C.M. Anderson-Cook . Response Surface Methodology: Proces and Product Optimization Using Designed Experiments . John Wiley and Sons, Hoboken, 2009 .

[41] R. Novak , L. Xiao , J. Lee , Y. Bahri , G. Yang , et al. Bayesian deep convolutional networks with many channels are Gaussian processes . In ICLR , pages 1 - 35 , 2019 .

[42] Y.S. Ong , P.B. Nair , A.J. Keane , and K.W. Wong . Surrogate-assisted evolutionary optimization frameworks for high-fidelity engineering design problems . In Y. Jin, editor, Knowledge Incorporation in Evolutionary Computation , pages 307 - 331 . Springer, 2005 .

[43] Z. Pitra , M. Hanuš , J. Koza , J. Tumpach , and M. Holenˇa . Interaction between model and its evolution control in surrogate-assisted CMA evolution strategy . In GECCO'21, page paper no. 358 , 2021 .

[44] Z. Pitra , J. Repický , and M. Holenˇa . Boosted regression forest for the doubly trained surrogate covariance matrix adaptation evolution strategy . In ITAT 2018 , pages 72 - 79 , 2018 .

[45] K. Rasheed , X. Ni , and S. Vattam . Methods for using surrogate modesl to speed up genetic algorithm oprimization: Informed operators and genetic engineering . In Y. Jin, editor, Knowledge Incorporation in Evolutionary Computation , pages 103 - 123 . Springer, 2005 .

[46] C.E Rasmussen and N. Hansen . GPML 4 .0. matlab toolbox. http://www.gaussianprocess.org/gpml/code/ matlab/doc/.

[47] E. Rasmussen and C. Williams . Gaussian Processes for Machine Learning . MIT Press, Cambridge, 2006 .

[48] A. Ratle . Kriging as a surrogate fitness landscape in evolutionary optimization . Artificial Intelligence for Engineering Design, Analysis and Manufacturing , 15 : 37 - 49 , 2001 .

[49] H. Ulmer , F. Streichert , and A. Zell . Evolution strategies assisted by Gaussian processes with improved pre-selection criterion . In IEEE CEC , pages 692 - 699 , 2003 .

[50] University of California in Irvine. Repository of machine learning databases . http://www.ics.uci.edu/ mlearn, 2016 .

[51] V. Volz , G. Rudolph, and B. Naujoks . Investigating uncertainty propagation in surrogate-assisted evolutionary algorithms . In GECCO'17 , pages 881 - 888 , 2017 .

[52] A.G. Wilson , Z. Hu , R. Salakhutdinov , and E.P. Xing . Deep kernel learning . In ICAIS , pages 370 - 378 , 2016 .

[53] Z.Z. Zhou , Y.S. Ong , P.B. Nair , A.J. Keane , and K.Y. Lum . Combining global and local surrogate models to accellerate evolutionary optimization . IEEE Transactions on Systems, Man and Cybernetics . Part C: Applications and Reviews, 37 : 66 - 76 , 2007 .