Landscape Analysis for Surrogate Models in the Evolutionary Black-Box Context Free

When solving a real-world optimization problem, we often have no information about the analytic form of the objective function. Evaluation of such black-box functions is frequently expensive in terms of time and money (Baerns and Holeňa, 2009; Lee et al., 2016; Thomaser et al., 2023; Zaefferer et al., 2016), which has been for two decades the driving force of research into surrogate modeling of black-box objective functions (Büche et al., 2005; Forrester and Keane, 2009; Jin, 2011). Given a set of observations, a surrogate model can be fitted to approximate the landscape of the black-box function.

The covariance matrix adaptation evolution strategy (CMA-ES) by Hansen (2006), which we consider the state-of-the-art evolutionary black-box optimizer, has been frequently combined with surrogate models. However, different surrogate models can show significant differences in the performance of the same optimizer on different data (Pitra et al., 2021). Such performance can also change during the algorithm run due to the varying landscape of the black-box function in different parts of the search space. Moreover, the predictive accuracy of individual models is strongly influenced by the choice of their parameters.

Our central goal in this investigation is to uncover and comprehend the intricate relationships between the configuration of individual models and data about properties and behavior of the objective function gathered during optimization runs and employed for surrogate model training. This data, representing a subset of solutions explored by the algorithm, constitutes the core of our analysis. The primary objective behind this endeavor is to facilitate the development of specific rules and automated approaches for selecting and tuning surrogate models. By gaining a deep understanding of these relationships, we aim to provide a foundation upon which adaptive and intelligent optimization strategies can be built.

Moreover, these relationships hold the promise of being valuable inputs for metalearning. Metalearning leverages insights from the performance of multiple optimization tasks to adapt and improve the optimization process. By characterizing the connections between surrogate model settings, model performance, and optimization outcomes, we lay the groundwork for the integration of these insights into the realm of metalearning, enabling the development of highly adaptable and efficient optimization strategies that learn and adapt from their experiences.

In recent years, research into landscape analysis of objective functions (cf. the overview by Kerschke, 2017) has emerged in the context of algorithm selection and algorithm configuration. Recently, Thomaser et al. (2023) have selected benchmark functions for tuning algorithm configuration based on their landscape features' similarity to real-world optimization problems in vehicle dynamics. However, to our knowledge, such features have been investigated in connection with surrogate models in the context of black-box optimization using static settings only, where the model is selected once at the beginning of the optimization process (Saini et al., 2019). Moreover, the analysis of landscape features also had less attention so far than it deserves and not in the context of surrogate models (Renau et al., 2019, 2020). Additionally, recent research by Saleem and Gallagher (2022) has proposed an alternative approach. While our emphasis lies in representing the objective function through landscape features, their approach suggests representing the same objective function using a GP regression model.

In this paper, we study properties of features representing the fitness landscape in the context of the data from actual runs of a surrogate-assisted version of the CMA-ES on the noiseless part of the comparing continuous optimizers (COCO) benchmark function testbed (Hansen et al., 2016). From the large number of available features, we select a small number of representatives for subsequent research, based on their robustness to sampling and similarity to other features. Basic investigation in connection with landscape features of CMA-ES assisted by a surrogate model based on Gaussian processes (Rasmussen and Williams, 2006) was presented in a conference paper by Pitra et al. (2019). Here, we substantially more thoroughly investigate the relationships between selected representatives of landscape features and the error of several kinds of surrogate models using different settings of their parameters and the criteria for selecting points for their training. Such study is crucial due to completely different properties of data from runs of a surrogate-assisted algorithm compared with generally utilized sampling strategies, which imply different values of landscape features as emphasized in Renau et al. (2020).

The paper is structured as follows. Section 2 briefly summarizes the basics of surrogate modeling in the context of the CMA-ES, surrogate models used in this investigation, and landscape features from the literature. In Section 3, properties of landscape features are thoroughly tested and their connection to surrogate models is analyzed. Section 4 concludes the paper and suggests a few future research objectives.

Black-box optimization problems frequently appear in the real world, where the values of a fitness function $B(x)$ can be obtained only empirically through experiments or via computer simulations. As recalled in the Introduction, such empirical evaluation is sometimes very time-consuming or expensive. Thus, we assume that evaluating the fitness represents a considerably higher cost than training a regression model.

As a surrogate model, we understand any regression model $B^:X→R$ that is trained on the already available input–output value pairs stored in an archive$A={(xi,yi)|yi=B(xi),i=1,...,N}$⁠, and is used instead of the original expensive fitness to evaluate some of the points needed by the optimization algorithm.

The CMA-ES proposed by Hansen and Ostermeier (1996) has become one of the most successful algorithms in the field of continuous black-box optimization. Its pseudocode is outlined in Algorithm 1.

After generating $λ$ new candidate solutions using the mean $m(g)$ of the mutation distribution added to a random Gaussian mutation with covariance matrix $C(g)$ in generation $g$ (Step 3), the fitness function $B$ is evaluated for the new offspring (Step 4). The new mean $m(g+1)$ of the mutation distribution is computed as the weighted sum of the $μ$ best points among the $λ$ ordered offspring $x1,...,xλ$ (Step 5).

The sum of consecutive successful mutation steps of the algorithm in the search space $(m(g+1)-m(g))/σ(g)$ is utilized to compute two evolution path vectors $pσ$ and $pc$ (Step 6). Successful steps are tracked in the sampling space and stored in $pσ$ using the transformation $C(g)-1/2$⁠. The evolution path $pc$ is calculated similarly to $pσ$⁠; however, the coordinate system is not changed (Step 8). The two remaining evolution path elements are a decay factor $cσ$ decreasing the impact of successful steps with increasing generations, and $μw$ used to normalize the variance of $pσ$⁠.

The covariance matrix adaptation (Step 10) is performed using rank-one and rank-$μ$ updates. The rank-one update utilizes $pc$ to calculate the covariance matrix $pcpc⊤$⁠. The rank-$μ$ update cumulates successful steps of $μ$ best individuals in matrix $Cμ$ (Step 9).

The step-size $σ$ (Step 11) is updated according to the ratio between the evolution path length $∥pσ∥$ and the expected length of a random evolution path. If the ratio is greater than 1, the step-size is increasing, and decreasing otherwise.

The CMA-ES uses restart strategies to deal with multimodal fitness landscapes and to avoid being trapped in local optima. A multi-start strategy where the population size is doubled in each restart $nr$ is referred to as IPOP-CMA-ES (Auger and Hansen, 2005).

The level of the model precision is highly determined by the selection of the points from the archive to the training set $T$⁠. Considering the surrogate models we use in this paper, we list the following training set selection (TSS) methods:

• TSS full, taking all the already evaluated points; that is, $T=A$⁠,

• TSS knn, selecting the union of the sets of $k$-nearest neighbors of all points for which the fitness should be predicted, where $k$ is user defined (e.g., in Kern et al., 2006),

• TSS nearest, selecting the union of the sets of $k$-nearest neighbors of all points for which the fitness should be predicted, where $k$ is maximal such that the total number of selected points does not exceed a given maximum number of points $N max$ and no point is further from current mean $m(g)$ than a given maximal distance $r max$ (e.g., in Bajer et al., 2019).

The models in this paper use the Mahalanobis distance given by the CMA-ES matrix $σ2C$ (see Equation 4) for selecting points into training sets, similarly to, for example, Kruisselbrink et al. (2010). Considering the fact that TSS knn is the original TSS method of the model from Kern et al. (2006) (see Subsection 2.1.3), we list it here and utilize it only in connection with this particular model. The particular parameters of the TSS methods employed in this study can be found in Subsection 3.3.1.

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 (Hosder et al., 2001; Myers et al., 2009). Therefore, it is not surprising that typical response surface models, that is, low order polynomials, also belong to most traditional and most successful surrogate models (Rasheed et al., 2005; Kern et al., 2006; Auger et al., 2013; Hansen, 2019). Here, we present two surrogate models using at most a quadratic polynomial utilized in two successful surrogate-assisted variants of the CMA-ES: the local-metamodel-CMA (lmm-CMA) by Kern et al. (2006) and the linear-quadratic CMA-ES (lq-CMA-ES) by Hansen (2019).

The linear-quadratic model in the algorithm lq-CMA-ES differs from the lmm-CMA in an important aspect:

Whereas surrogate models in the lmm-CMA are always full quadratic, the lq-CMA-ES admits also pure quadratic or linear models. The employed kind of model depends on the number of points in $T$ (according to the employed TSS method).

Regression forest (Breiman, 2001) (RF) is an ensemble of regression decision trees (Breiman, 1984). Recently, the gradient tree boosting (Friedman, 2001) has been shown useful in connection with the CMA-ES (Pitra et al., 2018). Therefore, we will focus only on this method.

Let us consider a sample set $S$ of $N$ pairs of observations in the context of continuous black-box optimization $S=(xi,yi)∈RD×R∪{∘}|i=1,...,N$⁠, where $∘$ denotes missing $yi$ value (e.g., $xi$ was not evaluated yet). Then the sample set can be utilized to describe landscape properties using a landscape feature$ϕ:⋃N∈NRN,D×R∪{∘}N,1↦R∪{±∞,•}$⁠, where $•$ denotes impossibility of feature computation.

A large number of various landscape features have been proposed in recent years in literature. Mersmann et al. (2011) proposed six easy to compute feature sets (each containing a number of individual features for continuous domain) representing different properties: y-Distribution set with measures related to the distribution of the objective function values, Levelset features capturing the relative position of each value with respect to objective quantiles, Meta-Model features extracting the information from linear or quadratic regression models fitted to the sampled data, Convexity set describing the level of function landscape convexity, Curvature set with gradient and Hessian approximation statistics, and Local Search features related to local searches conducted from sampled points. The last three feature sets require additional objective function evaluations.

The cell-mapping (CM) approach proposed by Kerschke et al. (2014) discretizes the input space to a user-defined number of blocks (i.e., cells) per dimension. Afterwards, the corresponding features are based on the relations between the cells and points within. Six cell-mapping feature sets were defined: CM Angle, CM Convexity, CM Gradient Homogeneity, Generalized CM, Barrier, and SOO Trees by Flamm et al. (2002) and Derbel et al. (2019). It should be noted that CM approach is less useful in higher dimensions where the majority of cells is empty and feature computation can require a lot of time and memory.

Nearest better clustering (NBC) features proposed by Kerschke et al. (2015) are based on the detection of funnel structures. The calculation of such features relies on the comparison of distances from observations to their nearest neighbors and their nearest better neighbors, which are the nearest neighbors among the set of all observations with a better objective value. Lunacek and Whitley (2006) proposed the set of dispersion features comparing the dispersion among the data points and among subsets of these points from the sample set. The information content features of a continuous landscape by Muñoz, Kirley et al. (2015) are derived from methods for calculating the information content of discrete landscapes. In Kerschke (2017), three other feature sets were proposed: Basic set with features such as the number of points, search space boundaries, or dimension; the proportion of Principal Components for a given percentage of variance; and coefficients of Linear Model fitted in each cell. A comprehensive survey of landscape analysis can be found, for example, in Muñoz, Sun et al. (2015). A more detailed description of feature sets used in our research can be found in the online available supplementary material.¹

In recent years, several surrogate-model approaches have been developed to increase the performance of the CMA-ES (cf. the overview in Pitra et al., 2017). Each such approach has two complementary aspects: the employed regression model and the so called evolution control (Jin et al., 2001) managing when to evaluate the model and the true objective function while replacing the fitness evaluation step (Step 4) in Algorithm 1. In Pitra et al. (2021), we have shown a significant influence of the surrogate model on the CMA-ES performance. Thus, we are interested in relationships between the prediction error of surrogate models using various features of the fitness function landscape.

First, we state the problem and research question connected with the relations between surrogate models and landscape features. Then, we present our set of new landscape features based on the CMA-ES state variables. Afterwards, we investigate the properties of landscape features in the context of the training set selection methods and select the most convenient features for further research. Finally, we analyze relationships between the selected features and measured errors of the surrogate models with various settings.

The problem can be formalized as follows: In a generation $g$ of a surrogate-assisted version of the CMA-ES, a set of surrogate models $M$ with hyperparameters $θ$ are trained utilizing particular choices of the training set $T$⁠. The training set $T$ is selected out of an archive$A$ (⁠$T⊂A$⁠) containing all points in which the fitness has been evaluated so far, using some TSS method (see Section 2.1.2). Afterwards, $M$ is tested on the set of points sampled using the CMA-ES distribution $X te =xk|xk∼Nm(g),σ(g)2C(g),k=1,...,α$⁠, where $α∈N$ depends on the evolution control. The research question connected to this problem is: What relationships between the suitability of different models for predicting the fitness and the considered landscape features do the testing results indicate?

To utilize additional information comprised in CMA-ES state variables, we have proposed a set of features based on the CMA-ES (Pitra et al., 2019).

In each CMA-ES generation $g$ during the fitness evaluation step (Step 4), the following additional features $ϕ$ are obtained for the set of points $X=xii=1N$⁠:

• Generation number$ϕgeneration CMA =g$ indicates the phase of the optimization process.

• Step-size$ϕstep_size CMA =σ(g)$ provides information about the extent of the approximated region.

• Number of restarts$ϕrestart CMA =nr(g)$ performed till generation $g$ may indicate landscape difficulty.

• Mahalanobis mean distance of the CMA-ES mean $m(g)$ to the sample mean $μX$ of $X$

  11$ϕmean_dist CMA (X)=(m(g)-μX)⊤ΣX-1(m(g)-μX),$

  where $ΣX$ is the sample covariance of $X$⁠. This feature indicates suitability of $X$ for model training from the point of view of the current state of the CMA-ES algorithm.

• Square of the $pc$evolution path length $ϕevopath_c_norm CMA =∥pc(g)∥2$ is the only possible non-zero eigenvalue of rank-one update covariance matrix $pc(g+1)pc(g+1)⊤$ (see Subsection 2.1.1). That feature providing information about the correlations between consecutive CMA-ES steps indicates a similarity of function landscapes among subsequent generations.

• $pσ$ evolution path ratio, that is, the ratio between the evolution path length $∥pσ(g)∥$ and the expected length of a random evolution path used to update step-size.

  12$ϕevopath_s_norm CMA =∥pσ(g)∥EN(0,I)=∥pσ(g)∥ΓD22ΓD+12.$

  It provides useful information about distribution changes.
• CMA similarity likelihood. This feature quantifies the degree of agreement between the distribution of the set of points $X$ and the current CMA-ES distribution in the present generation, expressed on the logarithmic scale. In particular, it expresses on the logarithmic scale the likelihood of the set of points $X$ concerning the CMA-ES distribution.

  13$*-12ptϕcma_lik CMA (X)=-N2Dln2πσ(g)2+lndetC(g)-12∑x∈Xx-m(g)σ(g)⊤C(g)-1x-m(g)σ(g).$

Based on the studies in Bajer et al. (2019) and Pitra et al. (2021), we have selected the Doubly Trained Surrogate CMA-ES (DTS-CMA-ES) (Pitra et al., 2016) as a successful representative of surrogate-assisted CMA-ES algorithms.

During the model training procedure in the DTS-CMA-ES, shown in general in Algorithm 2, a set of training points $T$ is selected and transformations of $X tr$ using matrix $1σC-1/2$ and $y tr$ to zero mean and unit variance (⁠$T=(X tr ,y tr )$⁠) are calculated in Steps 2 and 3 before the model hyperparameters $θ$ are fitted. In case of a successful fitting procedure, the resulting model is tested for constancy on an extra generated population due to the CMA-ES restraints against stagnation. However, if the fitting procedure is unsuccessful, such as encountering an error due to method failure, the model is considered not trained. The GP model utilizes the Matlab fmincon optimizer for hyperparameter fitting, and if fmincon fails due to an infeasible starting point or numerical problems, the CMA-ES is employed as an alternative (Bajer et al., 2019). Using the same transformation, the points in a testing set $X te$ are transformed before obtaining the model prediction of fitness values $y^ te$⁠, which is then inversely transformed to the original output space.

We accept hyperparameter fitting failures in our research because we prefer adopting established and published configurations to refining individual models. Specifically, for GPs and polynomial models, configurations were drawn from publications explicitly exploring model setups (Auger et al., 2013; Bajer et al., 2019; Hansen, 2019). For RFs, we followed Pitra et al. (2018), supplementing certain settings while actively investigating others (see Subsection 3.4.1). Although rerunning GP hyperparameter fitting could mitigate some failures, our intentional treatment of training as an automatic step aligns with our overarching goal: to comprehensively explore intricate relationships between the objective function features and the model.

To investigate landscape features in the context of surrogate-assisted optimization, we have generated² a dataset $D$ of sample sets using independent runs of the 8 model settings from Pitra et al. (2019) for the DTS-CMA-ES algorithm (Bajer et al., 2019; Pitra et al., 2016) on the 24 noiseless single-objective benchmark functions from the COCO framework (Hansen et al., 2016). All runs were performed in dimensions 2, 3, 5, 10, and 20 on instances 11–15. These instances all stem from the same base problem, and are obtained by translation and/or rotation in the input space and also translation of the objective values. The algorithm runs were terminated if the target fitness value $10-8$ was reached or the budget of 250 function evaluations per dimension was depleted. Taking into account the double model training in DTS-CMA-ES, we have extracted archives $A$ and testing sets $X te$ only from the first model training. The DTS-CMA-ES was employed in the overall best non-adaptive settings from Bajer et al. (2019). The TSS methods were set with parameter values: $k=⌈min{2m dim ,|A|m dim }⌉$⁠, where $m dim =D(D+3)/2+1$⁠, for TSS knn as suggested by Auger et al. (2013), and $r max =4Qχ2(0.99,D)$ and $N max =20·D$ for TSS nearest according to the best settings for higher dimensions in (Bajer et al., 2019).

To obtain 100 comparable archives and testing sets for landscape features investigation, we have generated a new collection of sample sets $D100$⁠, where points for new archives and new populations are created using the weighted sum of original archive distributions from $D$⁠. The $g$-th generated dataset uses the weight vector $w(g)=19(0,...,0g-3,1g-2,2g-1,3g,2g+1,1g+2,0g+3,...,0)⊤$⁠, which provides distribution smoothing across the available generations.³ The data from well-known benchmarks were also used by Saini et al. (2019) and Renau et al. (2020).

Because each TSS method in Subsection 2.1.2 results in a different training set $T$ using identical $A$⁠, we have performed all the feature investigations for each TSS method separately. By combining the two basic sample sets for feature calculation $A$ and $T$ with a population $P$ consisting of the points without a known value of the original fitness to be evaluated by the surrogate model, we have obtained two new sets $AP=A∪P$ and $TP=T∪P$⁠. Step 2 of Algorithm 2 performing the transformation of the input could also influence the landscape features. Thus, we have utilized either transformed and non-transformed sets for feature calculations, resulting in 8 different sample sets (4 in case of TSS full due to $T=A$⁠): $A$⁠, $A⊤$⁠, $AP$⁠, $AP⊤$⁠, $T$⁠, $T⊤$⁠, $TP$⁠, $TP⊤$⁠, where $⊤$ denotes the transformation.

These sets do not require additional evaluations of the objective function and can be computed even in higher dimensions. Some features were not computed due to the following reasons (see Section 2 in the supplementary material for feature definitions): From $Φ Basic$⁠, we have used only $ϕdim$ and $ϕobs$ because the remaining ones were constant on the whole $D100$ dataset. $ϕdim$ is identical regardless the sample set. $ϕobs$ and all features from $Φy-D$ were not computed using the transformation matrix because it does not influence their resulting values. The feature $ϕlin_simple_intercept MM ∈Φ MM$ was excluded because it is useless if fitness normalization is performed (see Step 3 in Algorithm 2). The remaining features from $Φ MM$⁠, all the features from $Φ NBC$⁠, and all three features from $Φy-D$ were not computed on sample sets with $P$ because it also does not influence their resulting values.

For the rest of the paper, we will consider features which are independent of the sample set (i.e., $ϕdim$ and 5 $Φ CMA$ features) as a part of $A$-based features only. This results in the total numbers of landscape features equal to 197 for TSS full (all were $A$-based) and 388 for TSS nearest and TSS knn (from which 191 features were $T$-based).

The results of experiments concerning landscape features are presented in Tables S1–S6 in the supplementary material. In the paper, we report only main highlights of the results in Tables 1 and 2.

Table 2:

Results of a medoid clustering into 14 clusters of the features for individual TSS methods, according to their Schweizer-Wolf measure $σ$⁠. The clusters are separated by horizontal lines and their medoid representatives are marked as gray lines. $N•$ denotes the minimal number of points for which the feature calculation resulted in $•$ in at most 1% of cases. The third column shows the robustness—the proportion of measured cases such that the feature values of samples from the same CMA-ES distribution did not differ more than 0.05.

View Large

First, we have investigated the impossibility of feature calculation (i.e., the feature value $•$⁠) for each feature. Such information can be valuable and we will consider it as a valid output of any feature. On the other hand, a large amount of $•$ values on the tested dataset suggests low usability of the respective feature. Therefore, we have excluded features yielding $•$ in more than 25% of all measured values, which were for all TSS methods the feature $ϕeps_s Inf$ on $AP$⁠, $AP⊤$⁠, $TP$⁠, and $TP⊤$ and for the TSS knn also the features $ϕratio_mean_02 Dis$⁠, $ϕratio_median_02 Dis$⁠, $ϕdiff_mean_02 Dis$⁠, $ϕdiff_median_02 Dis$ on $T$⁠, $T⊤$⁠, $TP$⁠, and $TP⊤$⁠, as well as $ϕquad_w_interact_adj_r2 MM$ on $T$ and $T⊤$⁠. This decreased the numbers of features to 195 for TSS full, 384 for TSS nearest, and 366 for TSS knn. Many features are difficult to calculate using low numbers of points. Therefore, for each feature we have measured the minimal number of points $N•$ in a particular combination of feature and sample set, for which the calculation resulted in $•$ in at most 1% of cases. All measured values can be found in Tables S1–S6 in the supplementary material and values of selected robust features (see below) in Table 2. For the calculation of most of the features, the CMA-ES default population size value in $2D$⁠: $N•=λ def =4+⌊3ln2⌋=6$⁠, or initial point plus doubled default population size in $2D$⁠: $N•=1+2λ def =13$ was sufficient in all sample sets indexed with $P$⁠. Considering sample-set-independent features, no points are needed, because the values concern the CMA-ES iteration or CMA-ES run as a whole, not particular points. As can be seen from Tables S1–S6 in the supplementary material, the quantile of function values used for splitting the sample set decreases as the dispersion features $Φ Dis$ require more points to be computable (⁠$ϕ≠•$⁠). For quantile values, see Section 2.3 in the supplementary material. Finally, $ϕlin_w_interact_adj_r2 MM$ and $ϕquad_w_interact_adj_r2 MM$ have also high values of $N•$⁠, which is plausible taking into account their descriptions (see Section 2.6 in the supplementary material).

To enhance the comparability of the investigated features, we opted for a tailored normalization approach. Standardization, which scales features to have a mean of 0 and a standard deviation of 1, was deemed unsuitable for our features due to several reasons. Firstly, the use of standardization would result in the loss of information about the overall measure of centrality and overall variability of each landscape feature. This information is preserved when we choose normalization over standardization. Secondly, the features encompass a wide spectrum of ranges, from relatively narrow intervals such as $[-1,1]$ for $ϕnb_cor NBC$ to broader spans like $[1.5·10-10,1.4·1015]$ for $ϕstep_size CMA$⁠. Additionally, some features even contain infinite values.

Unlike standardization, the proposed way of nonlinear normalization, employing the sigmoid function, does not depend on the particular distribution of the feature, with only one exception—the 0.01 and 0.99 quantiles, which are mapped to $ϕ norm (Q0.01)=0.01$ and $ϕ norm (Q0.99)=0.99$⁠. By using quantiles, we effectively cut off the tails of the feature distribution, mitigating the impact of outliers and infinite values. This approach transforms infinity values to 0 and 1, fostering a more comparable representation of features with substantial differences in value sets. This nonlinear transformation achieves a balance between trimming extreme values and retaining the core data, aligning with our objective to enhance feature comparability while accommodating the unique characteristics of the features.

We have tested the dependency of individual features on the dimension using feature medians from 100 samples for each distribution from $D100$⁠. The Friedman's test rejected the hypothesis that the feature medians are independent of the dimension for all features at the family-wise significance level 0.05 using the Bonferroni-Holm correction. Moreover, for most of the features, the subsequently performed pairwise tests rejected the hypothesis of equality of feature medians for all pairs of dimensions. There were only several features for which the hypothesis was not rejected for some pairs of dimensions (see Tables S1–S6 in the supplementary material). Therefore, the influence of the dimension on the vast majority of features is essential.

Our investigation into the impact of multiple landscape features on the predictive error of surrogate models necessitates a high level of robustness of these features against random sampling of points. Consistency is paramount in our analysis, as we aim to ensure that when the algorithm is in the same location under similar circumstances, the features provide comparable values or, at the very least, exhibit substantial similarity. This requirement stems from the problem statement, emphasizing the need for consistent feature behavior in identical or similar situations.

Additionally, it is essential to note that the features were initially designed for datasets generated using variations of uniform sampling and Latin hypercube design. Our current dataset, however, originates from a mixture of Gaussian distributions, each representing a distinct generation. This diversification in the sampling approach underscores the importance of robust features capable of maintaining consistency across different situations.

To have a robust set of $k$ independent features, which return identical values for input in 95% cases, we would like all features to be identical in $1000.95k%$ cases. Thus, for our dataset $D100$ with 100 samples for each distribution even a small $k$ requires all values to be identical, which is almost impossible to achieve for most of the investigated features. Therefore, we define feature robustness as a proportion of cases for which the difference between the 1st and 100th percentile calculated after standardization on samples from the same CMA-ES distribution is $≤0.05$⁠. Table 1 lists numbers of features achieving different levels of robustness. We have selected the robustness $≥0.9$ to be used for subsequent analyses. The robustness calculated for individual features is listed in Tables S1–S6 in the supplementary material and its values for features with robustness $≥0.9$ in Table 2. The chosen level of robustness excluded from further computations all features from $Φ NBC$ and $Φy-D$ for every TSS and all features from $Φ Inf$ for TSS full and TSS nearest. Both $Φ NBC$ and $Φy-D$ are rather sensitive to the input data, where the influence of non-uniform sampling of data is probably not negligible. Features from $Φ Inf$ have very varied robustness. Whereas $ϕeps_max Inf$ provides high robustness around 0.85 on transformed sample sets (up to 0.97 for TSS knn), computations of $ϕm0 Inf$ and $ϕeps_s Inf$ resulted in the robustness 0.004 and 0.016, respectively, which were the two lowest among all features. The majority of $Φ CMA$ features provided high robustness caused mainly due to the independence of most of the features on the sample set. Features from $Φ Lvl$ based on quadratic discriminant analysis (qda) showed nearly double robustness compared to the rest of $Φ Lvl$ features. The transformation of the input increased robustness of specific types of features from $Φ Dis$ and $Φ MM$⁠. In particular, it increased the robustness of difference-based features from $Φ Dis$ to more than 0.99 and also of features based on model coefficients from $Φ MM$⁠. TSS knn is more sample dependent, therefore, the number of points in a sample set can vary. This also decreases the robustness of some ratio-based features. On the other hand, coefficients of simple models from $Φ MM$ show robustness over 0.99 and $ϕcma_mean_dist CMA$ mostly over 0.9. The limited number of robust features supports the conclusion of Renau et al. (2020) that the distribution of the initial design notably influences landscape analysis. A noticeable dependence of robustness on which of the sets $A$⁠, $T$⁠, or $P$ is used was not observed.

The large number of features suggests that for the purpose of investigation of their relationships to surrogate models, they should be clustered into a smaller number of groups of similar features. To this end, we have performed agglomerative hierarchical clustering according to $1-σ SW (ϕi,ϕj)$⁠, where $σ SW$ is the measure $σ$ by Schweizer and Wolff (1981) and $ϕi$⁠, $ϕj$ are the vectors of all medians from $D100$ for the features $i$ and $j$⁠. To compensate for the ordering-dependency of agglomerative hierarchical clustering, we have performed 5 runs of clustering for each TSS method to find the optimal number of clusters. The number of clusters exceeding a threshold 0.9 for $σ SW$⁠, averaged over all 15 runs, was 14; that number was subsequently used as the value of $k$ for subsequent $k$-medoid clustering using again Schweizer-Wolff measure $σ$ as a similarity. The features that are medoids of those 14 clusters are listed in Table 2. Even such a small number of feature representatives can be sufficient for achieving excellent performance in a subsequent investigation (Hoos et al., 2018; Renau et al., 2021). A majority of clusters contain features from the same feature set. Sometimes, the whole cluster is composed of the same features calculated only on different sample sets which suggests that the influence of feature calculation on various sample sets might be negligible in those clusters. On the other hand, feature clusters for TSS knn often all share the base sample set (⁠$A$ or $T$⁠), or the same transformation even if the features are not from the same set of features. Considering the large numbers of available features, it is worth noticing that most of the medoids are identical or at least very similar for all TSS methods. TSS full and TSS nearest medoids share identical features, where only 4 (⁠$ϕobs$⁠, $ϕdiff_median_02 Dis$⁠, $ϕlda_qda_25 Lvl$⁠, $ϕquad_simple_cond MM$⁠) differ in the sample set (⁠$A$ vs. $T$⁠, $A⊤$ vs. $T⊤$⁠, $AP$ vs. $TP$⁠, and $A⊤$ vs. $T⊤$⁠, respectively). Notice that sample sets differ only in using $A$ or $T$⁠. Moreover, 10 out of 14 representatives (⁠$ϕdim$⁠, $ϕobs$⁠, $ϕevopath_s_norm CMA$⁠, $ϕrestart CMA$⁠, $ϕcma_lik CMA$⁠, $ϕdiff_median_02 Dis$⁠, $ϕdiff_mean_05 Dis$⁠, $ϕlda_qda_10 Lvl$⁠, $ϕlda_qda_25 Lvl$⁠, $ϕquad_simple_cond MM$⁠) are the same for all considered TSS methods, whereas sample sets utilized for feature calculation sometimes differ. Such similarity can indicate great importance of those features for characterizing the fitness landscape in the CMA-ES surrogate modeling context.

To investigate the relationships between surrogate model errors and 14 landscape features selected for each TSS in the previous subsection, we have utilized 4 surrogate models: lq, lmm, GPs, and RFs. The latter two models in 8 and 9 different settings respectively. In Saini et al. (2019), $Φ MM$ features were utilized to investigate connection with several surrogate models in default settings in the scenario, where the model is selected once for a specific problem regardless of the possible following optimization process.

We have split the dataset $D$ into validation and testing parts (⁠$D val$ and $D test$⁠) on the covariance function level uniformly at random (considering following levels of $D$⁠: dimension, function, instance, covariance function). More specifically, for each dimension, function, and instance in $D$⁠, we have uniformly selected runs using 7 covariance functions to $D test$ and 1 covariance to $D val$⁠. In each of those runs, data from 25 uniformly selected generations were used.

The regression model from lmm-CMA was used in its improved version published by Auger et al. (2013). The lmm model operates with the transformation matrix in its own way, thus, the transformation step during training (Step 2 in Algorithm 2) is not performed for this model.

The linear-quadratic model was used in the version published by Hansen (2019). The original version utilizes all data without a transformation; therefore, the input data for lq model are not transformed for TSS full.

The GP regression model was employed in the version published by Bajer et al. (2019) using 8 different covariance functions described in Table 3.

The experimental settings of the RF model involved 5 splitting methods obtained from algorithms CART (Breiman, 1984), SCRT (Dobra and Gehrke, 2002), OC1 (Murthy et al., 1994), PAIR (Chaudhuri et al., 1994), and SUPP (Hinton and Revow, 1996). These settings were determined using Latin-hypercube design on 100 samples from all combinations of the values in Table 4 for the number of trees in RF (⁠$n tree$⁠), the number of points bootstrapped out of $N$ training points (⁠$Nt$⁠), and the number of randomly subsampled dimensions used for training each individual tree (⁠$nD$⁠). The maximum tree depth was set to 8, following (Chen and Guestrin, 2016). Additionally, $α$ and $γ$ were set to 0, indicating no regularization is applied, aligning with traditional gradient tree boosting (Friedman, 2001). The remaining decision tree parameters were kept identical to settings from Pitra et al. (2018). The settings winning on the validation dataset $D val$ are shown in Table 5. Preliminary testing showed that RFs performance in connection with DTS-CMA-ES is higher when the input data are not transformed. Thus, the transformation step is omitted during model training for both TSS methods.

We have summed up the cases when the model did not provide the prediction; that is, the error value is not available, in Table 6. Such cases can occur when hyperparameters fitting fails, fitness prediction fails, or the model is considered constant (Steps 4, 7, and 9 in Algorithm 2, respectively).

The numbers in Table 6 basically confirm that more complex methods are more likely to fail. Whereas the lq model, the most simple model among all tested, provided predictions in almost 98% of all cases, where the missing results can be attributed to constant predictions, the GP model with Gibbs covariance function provided on average only 55% of the required predictions. Generally, the models were able to provide predictions more often when using TSS full than when using TSS nearest, probably due to the locality of the training set, considering that it is easier to train a constant model on a smaller area, where the differences between objective values are more likely negligible. The TSS knn was the most successful selection method for lmm, probably because it was designed directly for this model. In the following investigation, the cases where the error values were missing for a particular model were excluded from all comparisons involving that model.

We have tested the statistical significance of the $MSE$ and $RDE$ differences for 19 surrogate model settings using TSS full and TSS nearest methods and also the lmm surrogate model utilizing TSS knn, that is, 39 different combinations of model settings $ψ$ and TSS methods $(ψ, TSS )$⁠, on all available sample sets using the Friedman test. The resulting p-values of the two Friedman tests, one for each error measure, are below the smallest double precision number. A pairwise comparison of the model settings with respect to $MSE$ and $RDE$ revealed significant differences among the vast majority of pairs of model settings according to the non-parametric two-sided Wilcoxon signed rank test with the Holm correction for the family-wise error. To better illustrate the differences between individual settings, we also count the percentage of cases at which one setting had the error lower than the other. The pairwise score and the statistical significance of the pairwise differences are summarized in Table 7 for $RDE$ and Tables S7 and S8 in the supplementary material for both error measures. Results of statistical tests confirmed that the obtained $MSE$ and $RDE$ values are sufficiently diverse for further investigation of model settings suitability. The best overall results were provided by GPs with all covariances except LIN. Especially, GP using SE+Q as a covariance function in TSS nearest significantly outperformed all other $(ψ, TSS )$ combinations. The polynomial models using TSS nearest take the second place in such comparison, being outperformed only by the GP models mentioned above. Generally, models using TSS nearest provided better results than when using TSS full (in case of lmm also better than TSS knn). The percentages of $RDE$ show smaller differences in precision than $MSE$ due to the larger probability of error equality on the limited number of possible $RDE$ values compared to the infinite number of possible $MSE$ values. The significant variations in error values across all model settings, especially when employing different TSS methods for both error measures, emphasize the sensitivity of surrogate model performance to the selected configurations.

To compare the suitability of individual features as descriptors, we separated areas where the surrogate model $M$ with a particular $(ψ, TSS )$ combination has the best performance. Subseqently, we compared the distributions of each feature value on the whole $D test$ with a distribution of each feature value of the selected part using Kolmogorov-Smirnov test (KS test) for MSE and RDE separately. The hypothesis of distribution equality is tested at the family-wise significance level $α=0.05$ after the Holm correction. The resulting p-values summarized in Tables S9–S13 in the supplementary material and visualized in color in Figure 1 show significant differences in distribution among combinations of model settings and TSS method $(ψ, TSS )$ for the vast majority of considered features. Features from $Φ Dis$ and $ϕstep_size CMA$ provided the most significant differences for almost all models. For the GP models and the sample set leading to the lowest $MSE$⁠, the p-values were often even below the smallest double precision number. These features also provided very low p-values for lmm and lq model. The only exception is $ϕdiff_median_02 Dis (T⊤)$ providing notably higher p-values for all model settings. Moreover, features calculated on $T$-based sample sets for TSS nearest provided in most cases higher p-values than when calculated on $A$-based. From the model point of view, differences of RF model settings are much lower than the rest of the settings for all the tested features. This might suggest the lower ability to distinguish between the individual RF settings maybe due to the low performance of these settings on the dataset. The p-values for the $MSE$ and the $RDE$ also differ notably mainly due to the different ranges of values of these two error measures. As its consequence, p-values for $MSE$ more clearly suggest the distinctive power of individual features on model precision, regardless the fact that $RDE$ is more convenient for CMA-ES related problems. Overall, the results of the KS test have shown that there is no tested landscape feature representative for which the selection of the covariance function would be negligible for the vast majority of model settings. In essence, the individualized optimization approach, crucial for adapting model settings to observed distribution patterns, is strengthened by the persistent and robust influence of these patterns across diverse error measures. This consistency increases confidence in their universal applicability, extending across different methodologies and training set selection methods.