From Real World Data to Test Functions

CIplus Band 6/2016 From Real World Data to Test Functions Andreas Fischbach, Martin Zaefferer, Jörg Stork, Martina Friese, Thomas Bartz-Beielstein From Real World Data to Test Functions Andreas Fischbach, Martin Zaefferer, Jörg Stork, Martina Friese, Thomas Bartz-Beielstein SPOTSeven Lab, Dept. of Comp. Sci. and Eng. Sci. TH Köln E-Mail: {andreas.fischbach, martin.zaefferer, joerg.stork, martina.friese, thomas.bartz-beielstein}@th-koeln.de www.spotseven.de 1 Introduction When researchers and practitioners in the field of computational intelligence are confronted with real-world problems, the question arises which method is the best to apply. Nowadays, there are several, well established test suites and well known artificial benchmark functions available. However, relevance and applicability of these methods to real-world problems remains an open question in many situations. Furthermore, the generalizability of these methods cannot be taken for granted. Some preliminary ideas about generalizability are discussed in [1, 2]. This paper describes a data-driven approach for the generation of test instances, based on real-world data, as depicted in Figure 1. The test instance generation uses data-preprocessing, feature extraction, modeling, and parameterization. It was applied to several real-world scenarios, e.g., in the context of genetic programming [ 3 ]. In this work we apply this concept on a classical design of experiment real-world project and generate test instances for benchmarking, i.e., design of experiment methods and model fitness. But it can also be used to compare and analyze several surrogate techniques and optimization algorithms as well. In most cases, complex and expensive real-world problems do not provide sufficient data for comparison of methods. Thus, our goal is our goal is to create a toolbox containing multiple data sets of real-world projects. With that toolbox, researchers are granted access to both the data sets and the derived test functions. This work mainly focuses on the following questions: Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 1 Createinitialdesign Evaluateobjectivefunction Buildsurrogatemodel Evaluatemodelfitness Real‐worldDataset Createparameterizedmodel Generatetestinstances 2 3 4 5 1 Figure 1: Simplified test process of surrogate models fitted upon evaluation on initial designs with real-world data based test instances. An optional validation step can be added before the test instances are used in step 3. (Q-1) What is the characteristic of a variation of a certain model parameter in terms of the models fitness landscape? (Q-2) How can similarity between models be computed and what are useful thresholds to separate similar models from almost equal models on the one hand, and completely different models on the other hand? Considering an example application for the test function generator, the design of experiment used to gather the reference data set used in this work will be further analyzed: (Q-3) Which design of experiment works best for the underlying real-world problem? The remainder of this work is organized as follows. Section 2 gives a literature review. Section 3 describes the process of the test function generation and applied modeling techniques as well as model similarity measures used. Section 4 illustrates the reference data set and the example application, the evaluation of design of experiment methods. Section 5 concludes the work and gives a short outlook for future work. 2Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 2 Related Work In the research field of benchmarking real-parameter problems many contributions address the random generation of problem instances based on user defined parameters, like the Max-Set of Gaussian Landscape Generator [ 4 ], or the Krigifier [5], both using Gaussian Processes. The Max-Set of Gaussian Landscape Generator computes the upper envelope of 𝑚 weighted Gaussian process realizations and can be used to generate continuous, bound-constrained optimization problems. The landscape generator is parameterized to control, e.g., the number of Gaussian components and implicitly, the number of local optima, the occurrences and variations of hills and the peaks, and the global optimum. The Krigifier realizes a procedure for generating nonlinear objective functions. Its idea is based on the convenient supposition that objective functions are realizations of stochastic processes. The user specifies an underlying trend, a stochastic process and a finite number of points at which the process will be observed. The Krigifier creates a noise term and uses the trend and the noise term to produce an objective function. A completely different approach is addressed by the Real-Parameter BlackBox Optimization Benchmarking [ 6 ]. The organizers of this benchmarking challenge choose and implement a benchmarking function testbed, typically covering artificial test functions like Sphere, Rosenbrock, Rastrigin, etc. Participants then have to apply their black-box problem solvers and their results will be gathered and compared. All these approaches do not rely directly on real-world problems. Our work realizes one step towards closing a gap in means of validating methods regarding their applicability and generalizability in practical deployment. 3 Test Function Generation The simplified process of real-world data based test function generation is depicted in Figure 2. It consists of the following steps: 1. Process a data set of a real-world problem, X is the design matrix and Y is a vector of corresponding outcome values of the underlying process. Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 3 Industrial Data xi,y i xi,y i Level 0 model e.g. Kriging or linear model generate ↵↵ Level 1 model1 model1 ↵1 ↵1 Level 1 model2 model2 ↵2 ↵2 Level 1 modeln modeln ↵n ↵n define bounds generation Figure 2: Resulting hierarchy of the generation process of different models to be used as test function instances. 2. Build a model, further denoted as level 0 model, that is suitable for regression and interpolation purposes of the data. In this work the Kriging technique is used to build the model. 3. Create a parameter 𝛼 to vary the previously fitted level 0 model. The parameter 𝛼 is a scalar or vector, that perturbs the generated model. It may, e.g., define a change in parameters or other variables of the derived model. First, bounds are generated for 𝛼 , ensuring numerical robustness. Then, 𝑛 instances 𝛼𝑖 with 𝑖 = 1 , ..., 𝑛 are randomly created within the chosen bounds. Finally, a randomly selected subset of the desired number of 𝑚 < 𝑛 instances is chosen. This ensures that the instances which need to be solved by the evaluated methods, are never known in advance. 4. Apply the selected 𝛼 instances on the level 0 model to retrieve the desired number of 𝑚 level 1 model instances, each coupled with an 𝛼𝑗 , 𝑗 = 1 , ..., 𝑚 . The process ensures the fulfillment of the requirements on the similarity demands of the models by computing similarity measures and discarding infeasible 𝛼instances. The remainder of this section describes the Kriging modeling technique, 4Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 the model variation and the model similarity evaluation. 3.1 Kriging Often used for the purpose of regression and interpolation, Kriging is a modeling method based on Gaussian processes. In the following we will stick closely to the descriptions by Forrester et al. [ 7 ]. Further details can be found in their book. Given a set of 𝑛 solutions X = {x(𝑖)}𝑖=1...𝑛 in a 𝑘 -dimensional continuous search space with observations y = {𝑦(𝑖)}𝑖=1...𝑛 , Kriging tries to determine an expression for a predicted value at an unknown location by interpreting the observations y as realizations of a stochastic process. The stochastic process is defined by the set of random vectors Y = {𝑌 ( x(𝑖) ) }𝑖=1...𝑛 . The correlation of the random variables 𝑌 ( · )is modeled as follows [7]: cor [︁𝑌(x(𝑖)), 𝑌 (x(𝑙))]︁= exp ⎛ ⎝− 𝑘 ∑︁ 𝑗=1 𝜃𝑗|𝑥(𝑖) 𝑗−𝑥(𝑙) 𝑗|𝑝𝑗⎞ ⎠.(1) The matrix that collects correlations of all pairs { ( 𝑖, 𝑙 ) } is called the correlation matrix Ψ. It is used in the Kriging predictor ^𝑦(x) = ^𝜇+𝜓𝑇Ψ−1(y−1^𝜇),(2) where ^𝑦 ( x )is the predicted function value of a new sample x , ^𝜇 is the maximum likelihood estimate (MLE) of the mean and 𝜓 is the vector of correlations between training samples X and the new sample x . The width parameter 𝜃 = (𝜃1, . . . , 𝜃𝑗, . . . , 𝜃𝑘)𝑇 determines how far the influence of each sample point x spreads. In detail, the larger the width parameter is, the faster are the potential changes in the predicted value. The smaller the width parameter is, the slower are the potential changes in the prediction. Since there is one 𝜃𝑖 for each dimension, this parameter can control the activity in each dimension. The parameter 𝑝𝑗 is usually fixed at 𝑝𝑗 = 2, and defines the shape of the correlation function: At 𝑝𝑗 = 2, the correlation function is more smooth, whereas 𝑝𝑗 = 1 is less smooth. In case of noise, the parameter 𝜆 is added to the diagonal of the correlation matrix Ψ . This allows the model a more smooth fit through observations (regression), in contrast to the default which reproduces all training data exactly (interpolation). Classically, 𝜆 is used to deal with noisy data. But it can as well be used to smoothen more rugged fitness landscapes. Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 5 All model parameters are determined by Maximum Likelihood Estimation (MLE). For 𝜃, 𝜆 and 𝑝, MLE requires numerical optimization. 3.2 Model variations The main parameters controlling the behavior of the model are 𝜃 and 𝜆 . The major goal of the test function generator is the deployment upon real-world data, which is usually noisy. So the variation of the parameter 𝜆 is a natural choice at a first glance. In addition the variation of the width parameter 𝜃 seems important to change the model (slightly) by maintaining the general characteristic of the fitness landscape under certain circumstances. The bounds of the variations of the parameters have to be defined carefully, otherwise the model can show signs of degeneration. The test function generator will compute lower and upper bounds for 𝛼 according to the fitted level 0 model. The first component of 𝛼 represents the 𝜆 value and the remaining components represent the corresponding 𝜃 values. They will be added to their corresponding values of the level 0 model to retrieve the altered level 1 model. Afterwards the correlation Matrix Ψ is recalculated, so that the changes take effect before predictions are made. The test function generation returns the desired number of similar functions, randomly drawn from the search space defined by the level 0 model and the bounds for 𝛼 . If the given bounds do not allow the creation of sufficient feasible instances, an error message is produced. An example is shown in Fig. 3. Here, the one-dimensional function 𝑓(𝑥)=(−18𝑥−2)2sin(20𝑥−4) is first sampled by 11 equidistant points. The derived level 0 Kriging model has 𝜆 = 0 and 𝜃 = 100. To derive the level 1 models, the bounds for 𝑎𝑙𝑝ℎ𝑎 are set to 𝑎𝑙𝑝ℎ𝑎𝑙𝑜𝑤 = [0 ,− 90] and 𝑎𝑙𝑝ℎ𝑎ℎ𝑖𝑔ℎ = [1 , 900] so that 𝜆 will be set to values between zero and one, and 𝜃 between ten and one thousand. For demonstration purposes, the extreme values for 𝑎𝑙𝑝ℎ𝑎 are chosen, as shown above each plot in Fig. 3. The plot shows that different 𝛼 values affect the ruggedness of the function and vary the number of local optima. 6Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 0.0 0.2 0.4 0.6 0.8 1.0 −200 0 100 300 Level 0 Model x f(x) −200 0 100 300 y (x) 0.0 0.2 0.4 0.6 0.8 1.0 −200 0 100 300 Level 1 Model, α = [0.2,0] x f(x) −100 0 100 200 y (x) 0.0 0.2 0.4 0.6 0.8 1.0 −200 0 100 300 Level 1 Model, α = [1,0] x f(x) −50 0 50 100 y (x) 0.0 0.2 0.4 0.6 0.8 1.0 −200 0 100 300 Level 1 Model, α = [0.2,−90] x f(x) −10 0 10 20 30 y (x) 0.0 0.2 0.4 0.6 0.8 1.0 −200 0 100 300 Level 1 Model, α = [1,−90] x f(x) 0 5 10 15 20 y (x) 0.0 0.2 0.4 0.6 0.8 1.0 −200 0 100 300 Level 1 Model, α = [1,900] x f(x) −50 0 50 100 150 y (x) Figure 3: A level 0 Kriging model and several derived level 1 models. The solid line is the true function 𝑓(𝑥), circles indicate observations used to fit the level 0 model and the dashed line indicates the level 1 model ^𝑦(𝑥). Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 7 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Lambda variation Measure value 0.00 0.25 0.50 0.75 1.00 −2 0 2 4 Theta 1 variation Measure value Measure linetype mae pearson.r rmse spearman t.test.p Measure symbol mae pearson.r rmse spearman t.test.p Figure 4: Effects of the variation of model parameters on the chosen similarity measures MAE, RMSE, Pearsons r, Spearman and the t-test p value. Left: Effects of 𝜆changes. Right: Effects of 𝜃1changes. A detailed look at the result table from the screening experiment was taken, to decide which thresholds for each of the measures should be taken to discard test function instances. In the summary of the measures, depicted in Table 4, the distribution of the measures can be seen. The correlation coefficients does not become very low, with their mean at around 0 . 95. The p value of the t-tests distributes on the whole scale from 0to +1. Table 4: Summary of statistical measures of similarity gathered on screening experiments for 𝛼. rmse pearson.r spearman t.test.p Min. :0.0000 Min. :0.8732 Min. :0.8560 Min. :0.0000006 1st Qu.:0.1819 1st Qu.:0.9268 1st Qu.:0.9283 1st Qu.:0.0189555 Median :0.2430 Median :0.9524 Median :0.9564 Median :0.1431840 Mean :0.2376 Mean :0.9486 Mean :0.9506 Mean :0.3314180 3rd Qu.:0.2935 3rd Qu.:0.9729 3rd Qu.:0.9749 3rd Qu.:0.6384715 Max. :0.4120 Max. :1.0000 Max. :1.0000 Max. :1.0000000 A full linear model with interaction terms between all parameters and main 14 Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 0.00 0.25 0.50 0.75 1.00 −2 0 2 4 Theta 4 variation Measure value Measure linetype mae pearson.r rmse spearman t.test.p Measure symbol mae pearson.r rmse spearman t.test.p Figure 5: Effects of the variation of 𝜃4. Table 5: Thresholds for measures preventing too similar and degenerated test functions. rmse pearson.r spearman t.test.p lower bound 0.2 0.85 0.85 0.5 upper bound 0.35 0.92 0.92 0.9 effects up to second degree predicting the RMSE and Spearman correlation coefficient reveals an adjusted coefficient of determination of about 0 . 85 and 0.74 respectively. For the design of experiment application, the bounds for the measures were set es shown in Table 5. With these thresholds 20 random test functions were drawn. The impact of the chosen designs on the RMSE of the derived linear models compared to the test function are shown in Figure 6. It can be seen that the base design, applied to create the training data in the real-world project, could not be outperformed by any other design type. At a first glance it might be surprising, that designs with a larger design size, e.g. the Uniform 81 or even the LHD 81, perform worse. But these designs were not set up to fit second order linear regression models and would therefore not be the first choice in such a setup. The results would surely Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 15 base FFD 2 FFD 3 LHD 16 LHD 67 LHD 81 Uniform 16 Uniform 67 Uniform 81 0.00 0.05 0.10 0.15 0.20 RMSE Design type Figure 6: Resulting RMSE values of models based on different design of experiment methods each applied on the same set of 20 random test functions. look different if the level 2 model would as well be a Kriging model. 5 Conclusion In this work the generation process of test function instances based on a real-world industrial data set is described. Based on this data a Kriging model is fitted and altered according to a variation parameter 𝛼 , by adding its components to the corresponding model parameter. The generated instances will be discarded or kept according to thresholds for several model similarity measures. This leads finally to a test function instance pool that can be used to benchmark, i.e, design of experiment and modeling methods, for their practical use on the underlying problem. The major research Questions addressed in this work were: (Q-1) What is the characteristic of a variation of a certain model parameter in terms of the models fitness landscape? The variation of the Kriging model parameters 𝜆 and 𝜃 lead to altered models that were still correlated at rather high coefficient rates around 0 . 9. This was ensured by the definition of lower and upper bounds relative to the parameter values of the base level 0 16 Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 model. This kind of limit the amount of changes applicable to the model and can be problematic for example when a large number of test instances is needed. (Q-2) How can similarity between models be computed and what are useful thresholds to separate similar models from almost equal models on the one hand, and completely different models on the other hand? This can kind of simple be visually analyzed for one or two dimensional problems, but is getting harder in general for larger number of dimensions. Statistical measures can help to distinguish models that are too similar or too different. Further analysis of the resulting fitness landscapes in sense of, e.g., number of local optima or gradients, can help to judge the similarity of two different models. (Q-3) Which design of experiment works best for the underlying real-world problem The application of different design of experiment methods has shown, that the base design, used to gather the training data set, was not outperformed by any other design method. Even larger designs like a Full factorial designs at three levels per factor and Latin hypercube designs with 81 points could not dominate the results. This of course has to be further analyzed by altering the model techniques for the level 2 models. Even slight changes to the base design can now be analyzed and might lead to interesting new design points for the customer in future real-world experiments. Interesting future work include the shift of the similarity computation to the beginning of the instance generation. The Matrix computations to alter Kriging models can be regarded as time consuming, especially with high dimensional data. So it would be beneficial to compute an expectational value for the similarity of two models. In addition, fitness landscape analysis is an interesting topic to include in the work. The fitness landscapes comparison based on the RMSE or correlations of function values can be extended by interesting features to add or to avoid, e.g., the number of local optima, plane areas or large gradients. Finally, different variation methods on the models, e.g. rotating, scaling, distortion of the input space should be applied. Considering Kriging models, conditional simulation, could be applied to deliver possible realizations of a Gaussian process with a certain probability and therefore be very suitable Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 17 to deliver model variations based on real-world data. References [1] M. Chiarandini and Y. Goegebeur, “Mixed Models for the Analysis of Optimization Algorithms,” in Experimental Methods for the Analysis of Optimization Algorithms (T. Bartz-Beielstein, M. Chiarandini, L. Paquete, and M. Preuss, eds.), pp. 225–264, Germany: Springer, 2010. [2] T. Bartz-Beielstein, “How to Create Generalizable Results,” in Springer Handbook of Computational Intelligence, pp. 1127–1142, Berlin, Heidelberg: Springer Berlin Heidelberg, 2015. [3] O. Flasch, “A Modular Genetic Programming System,” May 2015. Dissertation zur Erlangung des Grades eines Doktors der Ingenieurwissenschaften der Technischen Universität Dortmund, Fakultät für Informatik. [4] M. Gallagher and B. Yuan, “A general-purpose tunable landscape generator,” IEEE Transactions on Evolutionary Computation, vol. 10, pp. 590–603, Oct. 2006. [5] M. W. Trosset, I. for Computer Applications in Science, and Engineering., The Krigifier [microform] : a procedure for generating pseudorandom nonlinear objective functions for computational experimentation / Michael W. Trosset. Institute for Computer Applications in Science and Engineering, NASA Langley Research Center ; National Technical Information Service, distributor Hampton, VA : Springfield, VA, 1999. [6] N. Hansen, A. Auger, S. Finck, and R. Ros, “Real-parameter black-box optimization benchmarking 2009: Experimental setup,” Tech. Rep. RR-6828, INRIA, 2009. [7] A. Forrester, A. Sobester, and A. Keane, Engineering Design via Surrogate Modelling. Wiley, 2008. [8] T. Chai and R. R. R. Draxler, “Root mean square error (RMSE) or mean absolute error (MAE)? – Arguments against avoiding RMSE in the literature,” Geoscientific Model Development, vol. 7, pp. 1247–1250, June 2014. 18 Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 [9] M. A. Rebolledo Coy, S. Krey, T. Bartz-Beielstein, O. Flasch, A. Fischbach, and J. Stork, “Modeling and Optimization of a Robust Gas Sensor,” in Bioinspired Optimization Methods and their Applications (G. Papa and M. Mernik, eds.), pp. 267–278, May 2016. [10] R. Tibshirani, “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society Series B (Methodolgical), vol. 58, no. 1, pp. 267–288, 1996. [11] A. Fischbach, J. Stork, M. Zaefferer, S. Krey, and T. Bartz-Beielstein, “Analyzing Capabilities of Latin Hypercube Designs Compared to Classical Experimental Design Methods ,” in 25. Workshop Computational Intelligence (F. Hoffmann and E. Hüllermeier, eds.), pp. 255–270, 2015. Proc. 26. Workshop Computational Intelligence, Dortmund, 24.-25.11.2016 19 Kontakt/Impressum Diese Veröffentlichungen erscheinen im Rahmen der Schriftenreihe "CIplus". Alle Veröffentlichungen dieser Reihe können unter https://cos.bibl.th-koeln.de/home abgerufen werden. Köln, Januar 2012 Die Verantwortung für den Inhalt dieser Veröffentlichung liegt beim Autor. Datum der Veröffentlichung: 14.11.2016 Herausgeber / Editorship Prof. Dr. Thomas Bartz-Beielstein, Prof. Dr. Wolfgang Konen, Prof. Dr. Boris Naujoks, Prof. Dr. Horst Stenzel Institute of Computer Science, Faculty of Computer Science and Engineering Science, TH Köln, Steinmüllerallee 1, 51643 Gummersbach url: www.ciplus-research.de Schriftleitung und Ansprechpartner/ Contact editor’s office Prof. Dr. Thomas Bartz-Beielstein, Institute of Computer Science, Faculty of Computer Science and Engineering Science, TH Köln, Steinmüllerallee 1, 51643 Gummersbach phone: +49 2261 8196 6391 url: http://www.spotseven.de eMail: [email protected] ISSN (online) 2194-2870