Raphaël Cosson
email: raphael.cosson@univ-lille.fr

Bilel Derbel
email: bilel.derbel@univ-lille.fr

Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille.fr

Sébastien Vérel

Hernán E Aguirre

Qingfu Zhang
email: qingfu.zhang@cityu.edu.hk

Kiyoshi Tanaka
email: ktanaka@shinshu-u.ac.jp

Sébastien Verel
email: verel@univ-littoral.fr

Hernan Aguirre

Cost-vs-accuracy of sampling in multi-objective combinatorial exploratory landscape analysis

Keywords: Applied computing → Multi-criterion optimization Multi-objective optimization, landscape analysis, automated algorithm selection, NK-landscapes

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

Context and motivations. Since the seminal work of Rice [START_REF] Rice | The Algorithm Selection Problem[END_REF], the design of high-level automated algorithm selection (AAS) approaches has received a growing attention. In the last decade, exploratory landscape analysis (ELA) [START_REF] Kerschke | Automated Algorithm Selection: Survey and Perspectives[END_REF][START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems by Combining Exploratory Landscape Analysis and Machine Learning[END_REF][START_REF] Mersmann | Exploratory Landscape Analysis[END_REF] has emerged as a state-of-the-art methodology supporting the combination of landscape feature extraction with the development of automated recommendation systems integrating the so-considered features on the basis of statistical or machine learning (ML) models. Such a methodology is especially sound in the context of blackbox optimization problems, having unknown characteristics, and for which different algorithms can expose different performance profiles as a function of the problem instance being solved. In this context, fitness landscape analysis [START_REF] Kauffman | The Origins of Order[END_REF] allows one to compute meaningful features characterizing the search difficulty of the (blackbox) landscape underlying a given instance; e.g., in terms of difficulty, ruggedness, or multi-modality. Following a standard supervised ML methodology, a model can then be trained on the basis of extensive experiments on a set of known instances. The training phase leads to a mapping of (pre-computed) landscape feature values to the performance of (pre-executed) algorithms. Given a new unseen instance, features can thus be computed and provided as input to the trained model in order to obtain a prediction about the most suitable algorithm for solving the instance. A key ingredient for the success of such an ELA approach relies on the ability to design cheap and meaningful features. In this paper, we are interested in studying the impact of feature extraction for multi-objective combinatorial optimization.

Independently of the nature of the optimization domain (i.e., continuous or discrete, single-or multi-objective), extracting blackbox features requires a particular sampling of solutions from the search space. This consists in evaluating the objective values of a carefullychosen set of solutions, on the basis of which a numerical statistic, constituting the feature value, is computed. The sampling has two main critical implications. First, the feature value can be different depending on the type and the size of the sampling. This has a direct consequence on the quality of the ML model, and hence on the accuracy of an AAS approach. Second, one has to accommodate the cost of sampling, since the sample size impacts the overall budget, in terms of the number of evaluations affordable to solve a new unseen instance. Consequently, computing features that are both as informative as possible, and as cheap as possible, is a critically important issue. Our work aims at pushing a step towards a better understanding of the impact of the sampling type and cost on the design of a successful multi-objective combinatorial ELA approach.

Related work. In the continuous domain, eliciting the impact of sampling is relatively well studied. On the one hand, different sampling techniques such as uniform sampling, latin hypercube, or Sobol sequence were studied in the past [START_REF] Kerschke | Automated Algorithm Selection on Continuous Black-Box Problems by Combining Exploratory Landscape Analysis and Machine Learning[END_REF][START_REF] Mary | A Survey of Advances in Landscape Analysis for Optimisation[END_REF][START_REF] Morgan | Sampling Techniques and Distance Metrics in High Dimensional Continuous Landscape Analysis: Limitations and Improvements[END_REF]. Very recently, it was shown that, quoting the authors in [START_REF] Renau | Exploratory Landscape Analysis is Strongly Sensitive to the Sampling Strategy[END_REF][START_REF] Renau | Expressiveness and robustness of landscape features[END_REF], "feature values are not absolute . . . [and] cannot be interpreted as stand-alone measures" independently of the method being used for sampling. Such a statement was supported by an analysis of the high-level properties of features, such as their expressiveness and their robustness, with respect to sampling. On the other hand, the computational effort needed for feature computation was studied in the past [START_REF] Kerschke | Low-Budget Exploratory Landscape Analysis on Multiple Peaks Models[END_REF][START_REF] Liefooghe | Landscape features and automated algorithm selection for multi-objective interpolated continuous optimisation problems[END_REF][START_REF] Renau | Exploratory Landscape Analysis is Strongly Sensitive to the Sampling Strategy[END_REF][START_REF] Renau | Expressiveness and robustness of landscape features[END_REF][START_REF] Saleem | Direct Feature Evaluation in Black-Box Optimization Using Problem Transformations[END_REF]. For instance, recommendations for computing cheap features are given in [START_REF] Kerschke | Low-Budget Exploratory Landscape Analysis on Multiple Peaks Models[END_REF]. Reviewing all the literature is out of the scope of this paper, since our focus is on multi-objective combinatorial optimization, but let us however remark that very few features can be interchangeably considered for continuous and discrete domains. In fact, sampling in discrete domains is fundamentally different from sampling in continuous domains, and so are the existing features.

In combinatorial domains, a fitness landscape is generally defined by specifying a neighborhood relation among solutions; e.g., bitflip for bit-strings or swap for permutations. One can find different features defined by assuming a full enumeration of the search space. Although this allows for a more fundamental understanding of combinatorial landscapes, it can only be performed for small landscapes, and it is not a realistic option for AAS. Hence, relying on a sampling method is mandatory. The most common sampling technique is to perform some particular walk over the landscape. That is, starting from an initial solution, some neighbors are evaluated, a new current solution is chosen, and so on until some condition is satisfied. The solutions collected during a walk, and their neighbors, are used to define the landscape features. State-of-the-art features are typically based on solutions sampled by either random or adaptive walks [START_REF] Kauffman | The Origins of Order[END_REF][START_REF]Recent Advances in the Theory and Application of Fitness Landscapes[END_REF][START_REF] Weinberger | Correlated and uncorrelatated fitness landscapes and how to tell the difference[END_REF]. Roughly speaking, this is intended to inform about the challenges that search algorithms have to face when exploring the landscape underlying a blackbox optimization problems, such as the landscape ruggedness or the distribution of local optima.

Although we can find a number of analyses [START_REF] Thomson | Clarifying the Difference in Local Optima Network Sampling Algorithms[END_REF][START_REF] Thomson | Inferring Future Landscapes: Sampling the Local Optima Level[END_REF][START_REF] Vanneschi | NK Landscapes Difficulty and Negative Slope Coefficient: How Sampling Influences the Results[END_REF] with respect to feature sampling in single-objective combinatorial optimization, very few lessons can be learnt when it comes to integrate those features in an effective AAS approach. In multi-objective combinatorial optimization, which is our main focus, few studies can be reported. In a state-of-the-art feature-based approach [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF], the cost of sampling is shown to represent a relatively low proportion of the overall budget dedicated to running the automatically-selected algorithm. In [START_REF] Liefooghe | Dominance, Indicator and Decomposition Based Search for Multiobjective QAP: Landscape Analysis and Automated Algorithm Selection[END_REF], a similar observation is also reported. Despite such observations, the question of how the sampling influences the accuracy of a feature-based ELA approach remains open.

Contribution and methodology. In this paper, we provide the first in-depth investigations on the impact of sampling in multiobjective combinatorial landscape analysis. We rely on the fact that the computation of multi-objective landscape features are mainly influenced by: the type of the walk, the length of the walk, and the explored portion of the neighborhood. We hence propose to conduct a systematic empirical analysis of the impact of each factor, with a particular focus on the sampling cost issues. More precisely, our contribution can be summarized as follows:

• Considering two conventional sampling methods from the literature, we adjust the amount of budget devoted to the feature computation by: (i) constraining the number of neighbors to be evaluated along the walk, and (ii) controlling the length of the walk when the sampling method allows us to do so. Such settings are then analyzed using an extensive number of features taken from two state-of-theart sets: dominance-and indicator-based features [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF], and decomposition-based features [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF]. • These cost-adjustable features are analyzed in a fine-grained manner by eliciting their correlation with global problem characteristics, and in a coarse-grained manner when integrated within an AAS approach. This allows us to contrast the amount of budget used for feature computation with respect to: (i) the ability of feature values to grasp the problem characteristics, and (ii) the gain in approximation quality when solving an unknown problem instance. • We report a comprehensive study using a broad range of 𝜌mnk-landscapes and a portfolio of three variants of the Moea/d algorithm. We show that different global trends of feature values and non-trivial cost-vs-accuracy trade-off(s) can be distinguished. Depending on the feature under consideration, increasing the walk length or increasing the number of visited neighbors do not always lead to more informative feature values. Besides, controlling the cost of features computation can improve the accuracy of the AAS task. Interestingly, this holds independently of whether the sampling cost is taken into account in the available budget, or not.

Outline. In Section 2, we provide the necessary background. In Section 3, we discuss the sampling type and cost. In Section 4, we provide a descriptive analysis of feature values. In Section 5, we report the impact of sampling on high-level ELA tasks. In Section 6, we conclude the paper.

MULTI-OBJECTIVE LANDSCAPE ANALYSIS 2.1 Multi-objective Optimization

A multi-objective combinatorial optimization problem (MCOP) can be defined by a set of m objective functions 𝑓 = (𝑓 1 , . . . , 𝑓 m ), and a discrete set 𝑋 of feasible solutions in the decision space. Let 𝑍 = 𝑓 (𝑋 ) ⊆ R m be the set of feasible outcome vectors in the objective space. To each solution 𝑥 ∈ 𝑋 is assigned an objective vector 𝑧 ∈ 𝑍 , on the basis of the vector function 𝑓 : 𝑋 → 𝑍 . When maximizing, an objective vector 𝑧 ∈ 𝑍 is dominated by a vector 𝑧 ′ ∈ 𝑍 iff ∀𝑚 ∈ {1, . . . , m}, 𝑧 𝑚 ⩽ 𝑧 ′ 𝑚 and ∃𝑚 ∈ {1, . . . , m} s.t. 𝑧 𝑚 < 𝑧 ′ 𝑚 . A solution 𝑥 ∈ 𝑋 is dominated by a solution 𝑥 ′ ∈ 𝑋 iff 𝑓 (𝑥) is dominated by 𝑓 (𝑥 ′ ). A solution 𝑥 ★ ∈ 𝑋 is Pareto optimal if there does not exist any other solution 𝑥 ∈ 𝑋 such that 𝑥 ★ is dominated by 𝑥. The set of all Pareto optimal solutions is the Pareto set. Its mapping in the objective space is the Pareto front. Computing a Pareto set approximation is a difficult task, for which multi-objective evolutionary algorithms (MOEAs) constitute an effective option [START_REF] Deb | Multi-Objective Optimization using Evolutionary Algorithms[END_REF].

Multi-objective Landscape Features

Characterizing the landscape of a MCOP and its impact on the performance of MOEAs is a well-established research area since several decades; see, e.g., [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF][START_REF] Garrett | Multiobjective Landscape Analysis and the Generalized Assignment Problem[END_REF][START_REF] Knowles | Instance Generators and Test Suites for the Multiobjective Quadratic Assignment Problem[END_REF][START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF]. One can distinguish three general categories of landscape features: (i) features based on problemspecific knowledge, (ii) blackbox global feature requiring the full enumeration of the search space, and (iii) blackbox local features relying on sampling the search space. We are specifically interested in the latter, since it enables the effective development of high-level feature-based ELA approaches.

We consider an extensive number of blackbox local features taken from two recent state-of-the-art studies [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF][START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF]. These features are based on the specification of a neighborhood relation N : 𝑋 → 2 𝑋 among solutions. Features are defined on the basis of a walk on the landscape. More formally, a walk is an ordered sequence of neighboring solutions denoted W |ℓ | := (𝑥 0 , 𝑥 1 , • • • , 𝑥 ℓ ), such that 𝑥 0 ∈ 𝑋 , and 𝑥 𝑡 ∈ N (𝑥 𝑡 -1 ) for every 𝑡 ∈ {1, • • • , ℓ } [START_REF] Kauffman | The Origins of Order[END_REF][START_REF]Recent Advances in the Theory and Application of Fitness Landscapes[END_REF][START_REF] Weinberger | Correlated and uncorrelatated fitness landscapes and how to tell the difference[END_REF]. For the sake of clarity, the description of the different methods that we consider to perform a walk is delayed to Section 4. We first provide an overview of the considered features, as summarized in Table 1.

Decomposition-based Features

Decomposition-based landscape features [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF] are defined by aggregating single-objective features obtained using objective space decomposition [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF]. Given 𝜇 weight vectors in the objective space and a scalarizing function, such as the Chebychev function [START_REF] Ehrgott | Multicriteria Optimization[END_REF], 𝜇 scalarizing (single-objective) sub-problems are defined. For each sub-problem, a single-objective feature is first computed on the basis of some walk W |ℓ | . For each single-objective feature, the so-obtained 𝜇 values are then aggregated to compute a statistic constituting the multi-objective feature value. Hence, one needs to specify how the single-objective features are defined and how their values are aggregated over sub-problems [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF]. For aggregating, four simple statistics are considered: the mean (avg), the standard deviation (sd), the first (p1) and second (p2) coefficients of a second order polynomial regression to fit the single-objective features as 1.

Four categories of single-objective features with respect to each sub-problem are considered. The first one, denoted by fv. ★ .★ in Table 1, informs about the distribution of the scalarized singleobjective fitness values along the walk W |ℓ | . More precisely, the mean (avg), the standard deviation (sd), the first auto-correlation coefficient (r1) [START_REF] Kauffman | The Origins of Order[END_REF], the kurtosis (kr) and the skewness (sk) of the fitness values of solutions collected during the walk W |ℓ | is computed. This is indicated using the letter 𝑠 in Table 1. For example, the feature fv.r1.avg refers to the average (over the sub-problems) of the first auto-correlation coefficient (over the sequence of solutions of the walk). The second category, denoted fd.★.★, computes for every solution 𝑥 𝑖 ∈ W |ℓ | , the mean fitness difference with its neighbors N (𝑥 𝑖 ). The mean fitness difference values are then combined using the same 𝑠 ∈ {avg, sd, r1, kr, sk} statistics. Similarly, in the third category, denoted in. ★ .★ for improving neighbors, for every solution in the walk, the proportion of neighbors that improve the solution's scalarized fitness value (w.r.t. the considered single-objective sub-problem) is computed. The last category, denoted spd. ★ .★, computes the maximum distance (normalized by 𝜇) between every two sub-problems improved by a solution in the neighborhood N (𝑥 𝑖 ) of every 𝑥 𝑖 ∈ W |ℓ | . The interested reader is referred to [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF] for more details. In total, we consider 4 (categories) × 5 (𝑠) × 4 (𝑎) = 80 decomposition-based landscape features.

Dominance-and Indicator-based Features

The second considered set contains landscape features based on the dominance relation or on the hypervolume indicator [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF]. They are organized in 8 categories in Table 1. For each category, a measure is computed with respect to the solutions collected by means of a walk W |ℓ | . Akin to decomposition-based features, the similar 𝑠 ∈ {avg, sd, r1, kr, sk} statistic is then applied to this measure in order to obtain the feature value. The categories from [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] are as follows: the proportion of dominated (inf.★), dominating (sup.★), and incomparable (inc.★) neighbors with respect to the current solution, the proportion of locally non-dominated neighboring solutions (lnd.★), the proportion of supported locally non-dominated neighboring solutions (lsupp.★), the hypervolume of a solution (hv.★), the solution's neighborhood hypervolume (nhv.★), and the difference of hypervolume between a solution and its neighbors (hvd.★). The interested reader can refer [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] for a detailed description. All these features are based on computing the corresponding measure with respect to every solution in the input walk W |ℓ | , and then combining the so-obtained measures using the 𝑠 statistic, exactly as described previously. In total, we consider 8 (categories) × 5 (𝑠) = 40 dominance-and indicator-based landscape features.

Finally, let us notice that we include one feature which depends explicitly on the type of the walk, namely the length of an adaptive walk (law in Table 1) as will be discussed later. This feature is known to inform about the multi-modality of the landscape [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF][START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF].

COST-ADJUSTABLE FEATURES 3.1 Random and Adaptive Walks

Two methods are generally used in the literature [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF][START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] in order to compute the input walk W |ℓ | , namely random walk and adaptive walk. In a random walk, there is no particular criterion to pick the next neighboring solution. Starting with a first (randomly chosen) solution 𝑥 𝑡 =0 , at each iteration 𝑡 ∈ {0, • • • , ℓ -1}, the next solution 𝑥 𝑡 +1 is chosen uniformly at random among the neighbors N (𝑥 𝑡 ) of solution 𝑥 𝑡 . The walk length ℓ is a user-defined parameter.

In an adaptive walk, the next solution is selected among the improving neighbors in N (𝑥 𝑡 ). In a multi-objective setting, improving neighbors are considered with respect to the dominance relation [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF][START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. At each step 𝑡 of the walk, 𝑥 𝑡 +1 is picked among neighbors that dominate 𝑥 𝑡 . If no such a solution exists, then the walk ends. Besides, an adaptive walk requires to pick one dominating solution. A common strategy is to pick the first neighbor that dominates the current solution. Unlike a random walk, the length ℓ of an adaptive walk is determined by the number of steps required to fall into a Pareto local optimal solution 𝑥 ℓ . As commented previously, the actual value of ℓ is known to depend on landscape characteristics, and is considered as a feature on its own (law).

All features, with the exception of fv.★.★ and hv.★.★, depend not only on the sequence of solutions visited during the walk W |ℓ | , but their computation also requires to evaluate the whole neighborhood N (𝑥 𝑖 ) for every 𝑥 𝑖 ∈ W |ℓ | . Hence, the feature computation requires Θ(ℓ • |N |) calls to the evaluation function, where |N | denotes the size of the neighborhood, and ℓ depends on the type of the walk as discussed previously.

Truncated Neighborhood-based Features

We are interested in studying the impact of the amount of budget devoted to feature computation. We apply simple modifications to the previous mechanisms in order to control more finely the feature cost. We consider two options: (i) to control the walk length ℓ, and (ii) to control the number of neighbors that need to be evaluated.

The first option is only possible when using a random walk, since the termination step of an adaptive walk cannot be controlled explicitly. This is however a natural option in order to adjust the overall budget which was not studied in a systematic manner. The second option applies for both random and adaptive walks. More specifically, for any integer value 𝑟 ≤ |N |, let us call a truncated neighborhood, denoted N 𝑟 , a neighborhood obtained from the original neighborhood N by considering solely 𝑟 neighbors, i.e., ∀𝑥 ∈ 𝑋 ,

N 𝑟 (𝑥) = {𝑦 1 , 𝑦 2 , • • • , 𝑦 𝑟 | ∀𝑗 ∈ [1.
.𝑟 ], 𝑦 𝑗 ∈ N (𝑥)}. As such, we use a truncated neighborhood having 𝑟 solutions sampled uniformly at random from N . Then, we compute the different features on the basis of the so-defined truncated neighborhood.

On the one hand, these modifications do not affect the sequence of solutions of a random walk. By contrast, an adaptive walk using N 𝑟 can stop earlier at any step 𝑡, if the 𝑟 randomly-selected neighbors do not include any locally dominating solution, and even if there exist some solutions in N (𝑥 𝑡 ) \ N 𝑟 (𝑥 𝑡 ) that dominate 𝑥 𝑡 . On the other hand, independently of the type of the walk (i.e., random or adaptive), for every solution 𝑥 𝑡 in the so-computed walk, solely the 𝑟 solutions of N 𝑟 (𝑥 𝑡 ) are evaluated when computing the statistic required for the final feature value. In other words, although we follow exactly the same feature specification as discussed previously, by replacing the (full) neighborhood by its 𝑟 -truncated variant, the so-computed feature values can be different.

Additionally, the previous simple modification reduces the feature budget to Θ(ℓ • 𝑟 ). However, it is not clear how this impacts the feature values, nor it is clear how much the so-extracted values are still meaningful. Besides, one may wonder whether a specific choice of the budget is to be preferred to accurately characterize the underlying landscapes. In fact, different combinations of 𝑟 and ℓ values could be considered, while still constraining the overall budget to the same pre-fixed value.

A DESCRIPTIVE PER-FEATURE ANALYSIS 4.1 Experimental Setup

4.1.1 Benchmarks. Following previous works [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF][START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF], we use 𝜌mnklandscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF] as a problem-independent model for constructing diverse multi-objective multi-modal benchmarks with objective correlation. The function vector 𝑓 = (𝑓 1 , . . . , 𝑓 𝑖 , . . . , 𝑓 m ) is defined as 𝑓 : {0, 1} n ↦ → [0, 1] m such that each objective 𝑓 𝑖 is to be maximized. The objective value 𝑓 𝑖 (𝑥) of a solution 𝑥 = (𝑥 1 , . . . , 𝑥 𝑗 , . . . , 𝑥 n ) is an average value of the contributions associated with each variable 𝑥 𝑗 . Given objective 𝑓 𝑖 , and variable 𝑥 𝑗 , a contribution function 𝑓 𝑖 𝑗 : {0, 1} k+1 ↦ → [0, 1] assigns a real value for every combination of 𝑥 𝑗 and its k epistatic interactions 𝑥 𝑗 1 , . . . , 𝑥 𝑗 k . The functions are defined by:

𝑓 𝑖 (𝑥) = 1 n n ∑ 𝑗=1 𝑓 𝑖 𝑗 (𝑥 𝑗 , 𝑥 𝑗 1 , . . . , 𝑥 𝑗 k ), ∀𝑖 ∈ {1, . . . , m}
The k interactions are set uniformly at random among the (n -1) variables other than 𝑥 𝑗 . By increasing the value of k, problem instances can be tuned from smooth to rugged. The 𝑓 𝑖 𝑗 -values follow a multivariate uniform distribution of dimension m, defined by an m × m positive-definite symmetric covariance matrix (𝑐 𝑝𝑞 ) s.t. 𝑐 𝑝𝑝 = 1 and 𝑐 𝑝𝑞 = 𝜌 for all 𝑝 ≠ 𝑞 where 𝜌 > -1 m-1 defines the correlation among the objectives. 4.1.2 Parameter setting. We use standard latin hypercube sampling to generate a set of 1 000 bi-objective instances spanning the following ranges: n ∈ {50, 51, . . . , 200}, k ∈ {0, 1, 2, . . . , 8} and 𝜌 ∈] -1, 1]. A standard bit-flip neighborhood N is considered, i.e., two solutions are neighbors if their Hamming distance is 1. We consider a truncated neighborhood N 𝑟 of size 𝑟 = 𝛼 • n, where 𝛼 ∈ {0.05, 0.1, 0.25, 0.5, 1}. Notice that 𝛼 = 1 corresponds to the full original neighborhood, whereas for 𝛼 = 0.05 a sample of 5% of the original neighborhood is explored. The two types of walks are experimented. For a random walk, the length is set in the range ℓ ∈ {5, 10, 25, 50, 100}. For each instance, we thus obtain 1 (rnd. walk) × 5 (𝑟 ) × 5 (ℓ) + 1 (adapt. walk) × 5 (𝑟 ) = 30 possible values for each of the 121 features.

Random Walk Analysis

We start by studying the impact of the different sampling configurations when a random walk is considered. We use two measures to elicit the relationship between the setting of ℓ and 𝑟 , and the behavior of the so-computed feature values. Firstly, we compute the deviation of feature values with respect to the highest budgetdemanding setting of the sampling. More precisely, for every instance, we compute the value 𝑉 max for every feature f, obtained when 𝛼 = 1 and ℓ = 100. Then, we compute for every other value 𝑉 𝑟,ℓ (f) of feature f, obtained with other settings of 𝑟 and ℓ, the relative deviation to 𝑉 max (f), that is, |𝑉 max (f) -𝑉 𝑟,ℓ (f)|/𝑉 max (f). Secondly, for each feature (in each sampling setting), we compute the Spearman correlation between the feature value and the value of 𝜌 and k over the considered 𝜌mnk-landscapes .

Preliminaries.

Before going into further details, we found that 30% of features are neither correlated to 𝜌 nor 𝑘, hence being of limited interest in the following analysis, given the extensive number of considered configurations. They mostly consist of features relating to the first auto-correlation coefficient ★.r1 or the kurtosis statistic ★.kr. Among the remaining features, we are able to elicit 4 classes where we observe specific trends of feature values and correlation as a function of the setting of ℓ and 𝑟 . In the first and second classes, containing respectively 20% and 15% of the overall features, feature values depend exclusively on either ℓ or 𝑟 . In the third and fourth classes, containing respectively 15% and 20% of the features, they expose a more complex dependency on the values of ℓ and 𝑟 . These four classes are illustrated in Fig. 1, reporting the average values of the relative deviation over the different instances (top), and the evolution of the Spearman correlation coefficients (bottom), as a function of the sampling cost. Due to space restrictions, only one selected feature from each class is shown in Fig. 1. Two features are taken from the decomposition-based set [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF], namely the average maximal distance between improving sub-problems spd.avg.avg and the average proportion of improving neighbors for each sub-problem in.avg.p2. The two other features are taken from the dominance-and indicator-based set [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF], namely the average number of incomparable solutions inc.avg and the standard deviation of the neighborhood hypervolume nhv.sd. This is discussed in more details below.

4.2.2

First and second feature classes. In the first class, illustrated by feature nhv.sd in Fig. 1 (first column), the feature values are mostly impacted by the length of the walk ℓ, while being insensitive to the value of 𝑟 . More precisely, the higher the value of ℓ, the smaller the relative deviation, independently of 𝑟 . This indicates that the features converge to some fixed values as a function of ℓ. However, the proportion of evaluated neighbors 𝑟 has almost no impact on the feature values, while implying a substantially higher costs. For instance, computing nhv.sd using a 5%-truncated neighborhood and ℓ = 100 achieves a relative deviation lower than 0.05%, while being 20 times less expensive than a full neighborhood exploration. Besides, our observations with respect to the relative deviation values stay consistent when looking at the Spearman correlations. For instance, the nhv.sd feature is mostly correlated with the parameter 𝜌 of the considered landscapes. The correlation changes in a consistent manner with respect to the relative feature deviation. In particular, the larger ℓ, the higher the correlation with 𝜌. However, the neighborhood proportion has almost no impact. Interestingly, exactly the opposite behavior can be reported for other features, which constitute our second class of features represented by inc.avg in Fig. 1 (second column). In this second class, feature values appear to converge to some value as a function of the neighborhood proportion 𝑟 . However, longer walks have a very small impact on the feature values while leading to a substantially higher cost. For instance, computing the inc.avg feature using the full neighborhood and a very small random walk of size ℓ = 5 achieves less then 2.5% deviation compared to the most expensive setting, while requiring 20 times less evaluations. We also observe that the evolution of the correlation coefficient with respect to 𝜌 follows the same trend as a function of ℓ and 𝑟 .

Third and fourth feature classes.

In the third class of features, represented by in.avg.p2 in Fig. 1 (third column), both ℓ and 𝑟 influences substantially the feature value. Roughly speaking, a larger budget, i.e., larger values of ℓ and 𝑟 , leads to lower relative deviations. However, this trend is not linear with ℓ and 𝑟 . For instance, for a small neighborhood proportion 𝛼 ∈ {0.1, 0.05}, a clear gap is observed with the other values of 𝛼 ∈ {0.25, 0.5, 1}. A similar observation can be made for relatively small walks of length ℓ ∈ {5, 10}. However, combining a random walk of length ℓ = 25 and a neighborhood proportion of 𝑟 = 25%, is enough to achieve a relative deviation lower than 3% while using 16 times less evaluations compared against the most expensive setting. Notice that the correlation coefficient increases consistently as commented previously for the other classes. Finally, in the fourth class of features, represented by spd.avg.avg in Fig. 1 (forth column), the relative deviation decreases consistently with the cost of sampling. However, the correlation coefficient does not follow the same trend in the sense that intermediate settings can provide substantially higher correlations. For instance, computing the spd.avg.avg feature using a small walk of length ℓ ∈ {5, 10} and the full neighborhood shows significantly larger correlation values compared against the most expensive setting, while requiring a significantly lower budget.

Discussion.

From the previous observations, we can say that there is a complex interaction between the feature values and the sampling configuration. In particular, a higher budget does not systematically lead to more consistent values, independently of the considered feature. Moreover, the settings of ℓ and 𝑟 can lead to seemingly different features exposing different trade-offs both in terms of feature values, feature cost, and feature ability to grasp the global characteristics of the underlying landscape. Hopefully, our analysis suggests that there could be some setting of ℓ and 𝑟 leading to reasonably cheap and accurate features. This will be studied in more details in Section 5. In the following, we first complement our analysis with respect to the adaptive walk.

Adaptive Walk Analysis

4.3.1 Feature values. Firstly, we report in Fig. 2 (left) the relative deviation of features with respect to the values obtained using an adaptive walk with the largest budget (ℓ = 100 and 𝑟 = 100%). We observe that the relative deviation spans a wide range (from 2% to more than 60%), as a function of 𝑟 . Actually, the largest deviations are with respect to very small 𝑟 values. This is to contrast with random walks where, although 𝑟 was found to have a significant impact, the range of the relative deviation is relatively small, which indicates that the size of the 𝑟 -truncated neighborhood is critical for an adaptive walk. In general, the feature values converge as 𝑟 increases. However, some features are found to be less sensitive to 𝑟 , as illustrated by the in.avg.p2 and the spd.avg features, for which a small value of 𝑟 allows to obtain a significantly lower deviation range. Overall, an adaptive walk visiting 25% of the neighborhood provides feature values which are relatively close to a full neighborhood (up to 5%), while requiring a significantly lower budget. One should however keep in mind that given the range of feature values, even such a small deviation is not necessarily negligible.

Correlation values.

Secondly, in Fig. 2 (middle and right), we report the evolution of the correlation between feature values and the benchmark parameters 𝜌 and k. In the situations where a feature shows a significant (positive or negative) correlation with either 𝜌 (e.g., rhv.sd or inc.avg) or k (e.g., spd.avg) when using a full neighborhood, the strength of the correlation decreases when extracting the feature with decreasing values of 𝑟 . For the other situations, no general tendency can be reported. From our collected data, we can state that computing features with an adaptive walk using the original neighborhood seems to expose the strongest correlation with the benchmark characteristics for the most part.

HIGH-LEVEL ELA TASKS

In the second part of our analysis, we consider two high-level ELA tasks: (i) predicting the global benchmark parameters 𝜌 and k, and (ii) selecting the (per-instance) best performing algorithm in a portfolio. For both tasks, we use the exact same set of 𝜌mnklandscapes described previously.

Task #1: Predicting Benchmark Parameters

5.1.1 Task description. This task informs about the power of a predictive model as a function of the sampling. We train one model with each possible configuration of the sampling in order to fit respectively the k and 𝜌 parameters of the input benchmarks. As in [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF], we use a random forest [START_REF] Breiman | Random Forests[END_REF] classification model for k, and a random forest regression model for 𝜌. Each model is trained on the basis of the input feature values using a default value of 100 trees [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. We then focus on the models' 𝑅 2 values, as summarized in Fig. 3. The higher the 𝑅 2 value of a model, the more explainable the variance, in the benchmark parameter 𝜌 or k, using the underlying sampling.

5.1.2 Results and Discussion. From Fig. 3, we observe that in comparison to a random walk, an adaptive walk provides a poor tradeoff in terms of the 𝑅 2 value and the implied budget for 𝑘. We also notice that the 𝑅 2 has a global tendency to increase as the budget underlying the sampling increases. However, reasonably high 𝑅 2 values are obtained at lower costs by adjusting the sampling parameters. For instance, an adaptive walk with half of the neighborhood (𝑟 = 50%) obtains almost the same 𝑅 2 values (0.94 for 𝜌 and 0.84 for k) while reducing the feature cost by half in average. Looking more carefully at the random walk configurations, we find a number of settings of ℓ and 𝑟 that obtain the same high level of accuracy, but with a drastic saving in the budget when compared against the most expensive setting. Interestingly, we found that the configuration providing the best trade-offs are different depending on the target prediction task. For k, a very small random walk (ℓ = 5) and the full neighborhood exploration (𝑟 = 100%) provides among the highest 𝑅 2 values (0.97) while requiring 20 times less budget than the most expensive setting. In fact, it appears that the most important sampling parameter for predicting k is the proportion 𝑟 of the neighborhood. For 𝜌, a configuration using a walk length ℓ = 50 and very small proportion of the neighborhood 𝑟 = 10% provides among the highest 𝑅 2 values (0.94) while requiring 20 times less budget than the most expensive setting. Actually, the 𝑅 2 value with respect to 𝜌 is mostly impacted by the length of the walk, provided that the proportion of evaluated neighbors is at least 𝑟 = 10%.

5.2 Task #2: Automated Algorithm Selection 5.2.1 Task description. The second and more sophisticated task aims at studying the impact of the sampling cost when tackling the AAS problem [START_REF] Rice | The Algorithm Selection Problem[END_REF] using a feature-informed ELA approach. For this purpose, we consider the exact same portfolio of algorithms as in [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF]. The reader is referred there for more details. For completeness, we comment that this portfolio is composed of three variants of the state-of-the-art Moea/d algorithm that are shown to have different performance profiles depending on the considered 𝜌mnklandscapes [START_REF] Marquet | Shake Them All! -Rethinking Selection and Replacement in MOEA/D[END_REF]. We measure algorithm performance in terms of hypervolume relative deviation [START_REF] Zitzler | Performance Assessment of Multiobjective Optimizers: An Analysis and Review[END_REF] to the best-known Pareto front collected over all runs for every instance. Following a standard supervised ML approach, we split the 1000 𝜌mnk instances into a training set T and a testing set I. We then train a multi-output random forest regression model [START_REF] Breiman | Random Forests[END_REF] in order to fit the relative hypervolume deviation of the different algorithms on the basis of the feature extracted with respect to each instance in the training set T . Afterwards, given an unseen instance from the testing set I, the features are first computed, and the performance of each algorithm is then predicted on the basis of the so-trained model. The algorithm with the best predicted performance is selected as the recommended one. Since extracting the features has a cost, we distinguish two scenarios for the purpose of our analysis. In the first scenario, the selected algorithm is run with the maximum budget B, hence not counting the sampling cost. In the second more realistic scenario, the algorithm is run with a budget of B -FB, where FB is the number of evaluations used to compute the features.

For performance assessment, we consider an adaptation of the merit measure [START_REF] Lindauer | The Algorithm Selection Competition Series[END_REF]. Let rhv(𝐴, 𝑖) be the average relative hypervolume deviation of a given algorithm 𝐴 using the maximum allowed budget B on an instance 𝑖. Similarly, let rhvc(𝐴, 𝑖) be the average relative hypervolume deviation of algorithm 𝐴 when subtracting the cost of feature computation, that is when running the algorithm with a budget of B -FB. Let rhvc(𝐴, 𝐽 ) and rhv(𝐴, 𝐽 ) be the average over a given set of instances 𝐽 of the relative hypervolume deviation, respectively when the feature cost is or not taken into account. A merit of 𝑚 = 0 corresponds to the perfect oracle. A merit m < 1 (respectively m > 1) indicates that using the model performs better (respectively worse) than SBS. The maximum budget is set to 10 6 evaluations. The parameters of the three Moea/d algorithms are exactly the same as in [START_REF] Cosson | Decomposition-Based Multi-objective Landscape Features and Automated Algorithm Selection[END_REF]. Each algorithm is executed for 20 independent runs on each instance. We adopt a standard repeated random hold-out strategy with a 90/10% split for training/testing, and we report the average merit value over 100 independent folds.

Results and Discussion

. In Table 2, we report the merit values obtained for the different configurations. First, we observe that the merit value m, taking into account the feature cost, is always worse than m ′ . With no surprise, this indicates that reducing the budget devoted to feature computation can help improving the overall performance. More importantly, the merit m is always lower than 1. This means that a feature-informed AAS approach is always better than SBS, independently of the experimented sampling configuration. This is especially interesting for the cheapest configuration (ℓ = 5, 𝑟 = 5%), implying an extremely constrained feature budget of at most 50 evaluations. Besides, a random walk provides better merit values than an adaptive walk. The merit increases when the sampling size decreases. More precisely, the loss of accuracy is slightly more important when the length of random walk is reduced. We also find that using the most expensive sampling configuration (ℓ = 100, 𝑟 = 100%) is not the most efficient in terms of merit. In fact, a random walk with ℓ = 100 and 𝑟 = 25% requires four time less budget while providing the best merit values. Surprisingly, this holds independently of whether the cost of sampling is taken into account or not. This gain can be attributed to the fact that a relatively small budget (i.e., 0.5% of the global budget) still allows to compute informative features leading to accurate predictions, while leaving more chance for the algorithm to converge. Which is consistent with the results from our first task, where we were able to elicit some cost-vs-accuracy trade-offs when varying the sampling parameters ℓ and 𝑟 for predicting the values of k and 𝜌.

Finally, we use the mean Gini impurity measure to extract the importance of each feature from the trained random forest models. We then compute the importance rank of each feature normalized in [0, 1], with 1 being the most important feature. In Fig. 4, we report the 10 most important features when using the most expensive setting (ℓ = 100, 𝑟 = 100%, left) and the sampling providing the best merit (ℓ = 100, 𝑟 = 25%, right). We clearly see that the ten most important features are the same in both settings. This indicates that using a restricted budget does not necessary lead to a major change in feature importance. However, we can observe small variations in the respective ranks of features. Besides, we found that five features fall into the first category elicited in the analysis of Section 4. In particular, the values of the four most important features hvd.sd, lsupp.sd, inc.sd, hv.sd are highly impacted by the length of the walk ℓ. One feature fd.avg.p2 falls into the second category where the proportion of neighbors 𝑟 is found to have the greatest impact on this feature value. Another feature inf.sd belongs to the third category, increasing its accuracy according to the sampling size. Finally, three features, namely hv.r1, in.sd.avg, and fv.r1.avg fall into the fourth category, where the highest correlation with the landscape characteristics is observed with a reduced sampling size.

CONCLUSION

In this paper, we conducted the first systematic analysis on the impact of sampling on the extraction of multi-objective combinatorial landscape features, and their integration into feature-informed performance prediction and algorithm selection approaches. From extensive experiments using state-of-the-art features, we provide evidence that feature values, and their correlation with unknown benchmark characteristics, expose complex dependencies with the sampling type and cost. In particular, we showed that a random walk using a reasonably small proportion of neighbors leads to cheap and informative feature values. Different questions are however left open. For instance, it would be interesting to extend our analysis to other algorithms, combinatorial domains and objective space dimensions, eventually including real-world problems, where the development of successful landscape-aware approaches is still not fully addressed.

Figure 1 :

 1 Figure 1: Average feature relative deviation (top), and feature correlation (bottom) with 𝜌 (first and second subfigure) and k (third and fourth subfigures) for random walks. From left to right: nhv.sd, inc.avg, in.avg.p2 and spd.avg.avg. The x-axis (in log-scale) refers the number of evaluations of the sampling.

Figure 2 :

 2 Figure 2: Average feature relative deviation (left), and feature correlation with 𝜌 (middle) and k (right) for adaptive walks. The x-axis (in log-scale) refers the number of evaluations of the sampling.

Figure 3 :

 3 Figure 3: Average 𝑅 2 (over 50 repetitions) of random forest models for predicting 𝜌 (left) and 𝑘 (right). The x-axis (in log-scale) refers the number of evaluations.

  The so-called single best solver (SBS) is the algorithm with the best rhv value on the training set T ; i.e., 𝑆𝐵𝑆 = arg min 𝐴 rhv(𝐴, T ). Notice that the SBS does not need any feature computation. The virtual best solver (VBS) is the oracle providing the best rhv(•, 𝑖, B) value for each testing instance 𝑖 ∈ I. Let us call the Automated Solver (AS) the algorithm selected by the trained model. The merit measures the gain in quality of AS relatively to SBS and VBS. For the purpose of our analysis, we consider two variants of the merit where the sampling cost is either taken into account (m) or not (m ′ ): m = rhvc(AS, I)rhv(VBS, I) rhv(SBS, I)rhv(VBS, I) ; m ′ = rhv(AS, I)rhv(VBS, I) rhv(SBS, I)rhv(VBS, I)

Figure 4 :

 4 Figure 4: Top 10 most important features. Line styles and shapes refer to the feature categories.

Table 1 :

 1 Summary of the 121 considered multi-objective landscape features (𝑠 ∈ {avg, sd, r1, kr, sk}, 𝑎 ∈ {avg, sd, p1, p2}).

	decomp. [3]	fv.𝑠.𝑎 fd.𝑠.𝑎 in.𝑠.𝑎 spd.𝑠.𝑎	scalar fitness values fitness difference between a solution and its neighbors proportion of improving neighbors maximum distance between improving sub-problems
	dominance [12]	inf.𝑠 sup.𝑠 inc.𝑠 lnd.𝑠 lsupp.𝑠	proportion of dominated neighbors proportion of dominating neighbors proportion of incomparable neighbors proportion of locally non-dominated neigh-bors proportion of supported locally non-dominated neighbors
	indic. [12]	hv.𝑠 hvd.𝑠 nhv.𝑠	solution's hypervolume neighborhood's hypervolume difference of hypervolume between a solution and its neighbors
	[3, 12] law	length of adaptive walk

a function of the ordered indices of sub-problem weight vectors. This is indicated using the letter 𝑎 in the notation of Table

Table 2 :

 2 Merit average values for the different sampling configurations.

			Merit with feature cost (m)			Merit without feature cost (m ′ )
			Random Walk		Adaptive		Random Walk		Adaptive
	𝑟										Walk
	100%	0.537	0.529	0.562	0.577 0.576	0.624	0.513	0.517	0.555	0.573 0.573	0.617
	50%	0.521	0.552	0.580	0.620 0.577	0.610	0.507	0.544	0.575	0.616 0.574	0.606
	25%	0.515	0.566	0.569	0.640 0.634	0.629	0.507	0.561	0.566	0.637 0.631	0.626
	10%	0.544	0.583	0.566	0.609 0.585	0.672	0.538	0.579	0.563	0.606 0.583	0.670
	5%	0.551	0.575	0.552	0.639 0.646	0.692	0.547	0.572	0.550	0.637 0.644	0.690

ℓ = 100 ℓ = 50 ℓ = 25 ℓ = 10 ℓ = 5 Walk ℓ = 100 ℓ = 50 ℓ = 25 ℓ = 10 ℓ = 5


Roberto Santana
email: roberto.santana@ehu.eus

Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille.fr

Bilel Derbel
email: bilel.derbel@univ-lille.fr

Boomerang-shaped Neural Embeddings for NK Landscapes

Keywords: • Computing methodologies → Discrete space search; Neural networks NK landscapes, fitness landscape, local optima, adaptive walk, neural embedding

Understanding the landscape underlying NK models is of fundamental interest. Different representations have been proposed to better understand how the ruggedness of the landscape is influenced by the model parameters, such as the problem dimension, the degree of non-linearity and the structure of variable interactions. In this paper, we propose to use neural embedding, that is a continuous vectorial representation obtained as a result of applying a neural network to a prediction task, in order to investigate the characteristics of NK landscapes. The main assumption is that neural embeddings are able to capture important features that reflect the difficulty of the landscape. We propose a method for constructing NK embeddings, together with metrics for evaluating to what extent this embedding space encodes valuable information from the original NK landscape. Furthermore, we study how the embedding dimensionality and the parameters of the NK model influence the characteristics of the NK embedding space. Finally, we evaluate the performance of optimizers that solve the continuous representations of NK models by searching for solutions in the embedding space.

CCS CONCEPTS

• Theory of computation → Optimization with randomized search heuristics;

INTRODUCTION

Neural embeddings [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF][START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF] have been mostly used as a vectorial representation of text data, amenable for the application of machine learning techniques. In this context, a word embedding is an -dimensional real-valued vector that represents a word from a vocabulary. An advantage of word embedding is that it makes possible to encode semantic relationships between words [START_REF] Antoniak | Evaluating the stability of embedding-based word similarities[END_REF][START_REF] Bolukbasi | Man is to computer programmer as woman is to homemaker? debiasing word embeddings[END_REF][START_REF] Faruqui | Problems with evaluation of word embeddings using word similarity tasks[END_REF]. That is, the proximity between vectors in the embedding space can be interpreted as a semantic similarity of the corresponding words these vectors codify. Word embeddings can also be used to solve supervised and unsupervised machine learning problems, for instance, as inputs to a model that predicts the class of a given word, or to estimate a continuous value associated with the word. Recently, neural embeddings have also been applied to other types of objects such as graphs [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF] and biological sequences [START_REF] Bepler | Learning protein sequence embeddings using information from structure[END_REF][START_REF] Melidis | dom2vec: Unsupervised protein domain embeddings capture domains structure and function providing data-driven insights into collocations in domain architectures[END_REF].

The NK model for fitness landscape has been extensively applied in evolutionary biology and evolutionary optimization [START_REF] Aguirre | Genetic algorithms on NK-landscapes: Effects of selection, drift, mutation, and recombination[END_REF][START_REF] Manderick | The genetic algorithm and the structure of the fitness landscape[END_REF][START_REF] Pelikan | Analysis of estimation of distribution algorithms and genetic algorithms on NK landscapes[END_REF][START_REF] Verel | Where are bottlenecks in NK fitness landscapes?[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF]. It provides a relevant benchmark to investigate how the ruggedness of the fitness landscape, and consequently the difficulty of the optimization process, are governed by the model parameters.

The main goal of this paper is to introduce alternative landscape representations of NK models, namely, by learning neural embedding representations of NK landscapes, and to study whether and how the features of these embedding spaces provide valuable information about the characteristics of the original landscape. We thereby consider neural networks to associate an -dimensional continuous representation with each binary solution from an NK landscape. We evaluate how the properties of this space are related to the characteristics of the original NK landscape. More particularly, we evaluate how neighbors in the discrete space are located in the embedding space, and the way in which a path that goes from a solution to the local optimum of its basin of attraction is projected into the embedding space. We also explore the way in which the embedding dimensionality ( ) and the degree of nonlinearity of NK landscapes ( ) influence the properties of the soobtained NK embedding spaces. Finally, we investigate the application of numerical optimization methods in the embedding space, as a way to asses the difficulty of the original problem. Extensive experiments are conducted using the CMA-ES algorithm [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF].

The remainder of the paper is organized as follows. Section 2 presents some background related to NK landscapes and neural embeddings. Section 3 gives an overview of related work. Section 4 explains the steps we follow to learn the NK model embedding representations, and the metrics we consider to evaluate them. Section 5 lists the experimental set-up, describes the numerical results and discusses the results of the experiments. Finally, Section 6 concludes the paper and discusses future work.

BACKGROUND

In this section we briefly present the main concepts related to NK landscapes and neural embeddings.

NK Landscapes

Let X = ( 1 , . . . ,

) denote a vector of binary variables. We will use x = ( 1 , . . . , ) to denote an assignment to the variables.

will denote a set of indices in {1, . . . , }, and (respectively ) a subset of the variables of X (respectively x) determined by the indices in .

The NK model was originally introduced to study the effect of different epistatic patterns in the landscape of the function. An NK landscape [START_REF] Kauffman | Origins of Order[END_REF] is defined by the following components:

• The number of variables, , that defines the problem dimension.

• The number of interactions per variable, , that defines the problem non-linearity. • A set of interactions, Π( ) ⊆ X, for every , ∈ {1, . . . , }.

• A sub-function defining a real value for each combination of values of and Π( ), ∈ {1, . . . , }.

The fitness function , to be maximized, is then defined as follows:

(x) = 1 =1 ( , Π( )). (1)
By increasing the number of interactions per variable, NK landscapes can be gradually tuned from smooth to rugged.

Neural Embedding

Neural embedding appeared as a straightforward way to apply machine learning algorithms to problems with a sparse representation [START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF]. In natural language processing, the problem arises from the issue of how to represent words and sentences so that they could be processed by machine learning models such as neural networks. A natural way to represent a sentence is using a binary vector where each word from a vocabulary is either present ( = 1) or not ( = 0). Since only a very small fraction of the words in a vocabulary are typically present in a sentence, this bag-of-words representation is extremely sparse. By contrast, word embeddings are continuous vector representations of words which in comparison with the bag-of-words representation has a lower dimensionality, and are denser. In a word embedding vector space, each word in a vocabulary has an associated real-valued vector.

To learn word embeddings, neural networks are used. First, a vocabulary of words is defined, and a corpus of sentences that contains words from this vocabulary is selected. Then, a machine learning task receives the bag-of-words representation of the sentence as input. Once the network is trained with the corpus, the set or word embeddings is extracted from the weights of the neural network.

For instance, in the word2vec approach [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF], a feed-forward neural network language model [START_REF] Bengio | A neural probabilistic language model[END_REF] with a number of added changes is used to learn two types of embedding models, the Skip-gram, and the continuous bags-of-words (CBOW) models. Both models use a bag-of-words representation for each word; i.e., a binary vector of dimension equal to the vocabulary size and where only one element is set to 1. The difference between these two models lies in the prediction task they address: Skip-gram learns to predict the surrounding words of a given word ( ) in a sentence, while CBOW learns to predict the word that is most likely to be in the center, given the surrounding words. The surrounding words of ( ) are those fromto + , where is an algorithm parameter that determines the window size [START_REF] Akiba | Optuna: A next-generation hyperparameter optimization framework[END_REF].

Figure 1 shows a Skip-gram model for a vocabulary of size 6, an embedding dimension = 3, and a window size of 2. Each unit of the input layer, corresponding to the -th word in the vocabulary, is connected to all units of the (single) hidden layer. The word embedding corresponding to word , extracted after the neural network training is completed, has in position the -th weight associated to the connection from to the -th hidden unit. In the example shown in Figure 1, the input word is the third word of the vocabulary, and its context in the sentence is formed by the sixth and second words of the vocabulary, highlighted in the bag-of-words representations of words ( -1) and ( + 1).

Beyond constituting a more compact representation, word embeddings have been shown to capture a number of regularities from the corpus under which they were trained [START_REF] Goldberg | word2vec explained: Deriving Mikolov et al. 's negativesampling word-embedding method[END_REF][START_REF] Lastra-Díaz | A reproducible survey on word embeddings and ontology-based methods for word similarity: linear combinations outperform the state of the art[END_REF][START_REF] Zhai | Intrinsic and extrinsic evaluations of word embeddings[END_REF]. For instance, words whose vectors are close in the embedding space are usually related semantically. Furthermore, some algebraic operations between vectors can be used to infer similarity and semantic relationships between words. For instance, since the distance between embeddings for words man and woman is similar to the distance between embeddings for words king and queen, it is possible to define algebraic equations such as

# » = # » -# » + # » ,
where #» indicates the embedding representation for word w.

RELATED WORK

In this section, we review related work in three areas: representations proposed to study NK landscapes, applications of evolutionary algorithms on embedding representations, and proposals of methods and algorithms on embedding representations for other domains.

Representation for NK Landscapes

Buzas and Dinitz [START_REF] Buzas | An analysis of NK landscapes: Interaction structure, statistical properties and expected number of local optima[END_REF] proposed a characterization of NK landscapes using parametric linear models comprised of main effects and interaction effects. Based on this "analytical" representation, they computed the statistical distribution of the model effect coefficients, and showed that the model coefficients have meaningful interpretations.

In their seminal work, Ochoa et al. [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF] introduced a network characterization of NK landscapes. In this network representation, each vertex corresponds to a local optimum, and edges represent basin adjacency between two local optima. The representation is compact since non-optimal solutions are discarded. The analysis of local optima networks (LONs) shows the existence of small-world characteristics. The network representation additionally allows for the extraction of general network metrics that serve to compare how different characteristics of the fitness landscapes change according to the parameters of the NK model. A number of extensions and variations of LONs have been proposed in the literature [START_REF] Ochoa | First-improvement vs. best-improvement local optima networks of NK landscapes[END_REF][START_REF] Tomassini | Complex-network analysis of combinatorial spaces: The NK landscape case[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF]. Most relevantly, these models have been instrumental in the identification of features inferring on problem difficulty and on algorithm performance prediction [START_REF] Chicano | Local optima networks, landscape autocorrelation and heuristic search performance[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF].

Such a network representation approach has been later extended to study multi-objective NK landscapes. The Pareto local optimal solutions network (PLOS-net) is introduced in [START_REF] Liefooghe | On Pareto local optimal solutions networks[END_REF]. Similarly to LONs, sub-optimal solutions are not represented in the network, vertices correspond to PLOs and edges are constructed between mutually non-dominated neighbors. Network metrics are extracted from PLOS-nets and used to explain the performance of different types of multi-objective algorithms.

The work presented in this paper have a number of similarities and differences from network-based approaches. Similarly to LONs, and as it will be discussed in Section 4, local optima play an important role in the construction of the embedding representations. Exhaustive search or sampling are required to construct a set of paths that, in each case, lead to a local optimum. In contrast with LONs, the representation of solutions lies in a real-valued space.

Evolutionary Algorithms and Embedding Representations

A number of works have investigated the combination of evolutionary computation and embedding representations. In [START_REF] Manzoni | Towards an evolutionary-based approach for natural language processing[END_REF][START_REF] Santana | Reproducing and learning new algebraic operations on word embeddings using genetic programming[END_REF][START_REF] Santana | Semantic composition of word-embeddings with genetic programming[END_REF], word embedding representations learned by means of word2vec are used as inputs to evolve genetic programs which are able to solve different natural language processing tasks. The dimension of the embeddings and the choice of the genetic programming fitness function were shown to play an important role in algorithm performance. Roman et al. [START_REF] Roman | Sentiment analysis with genetically evolved Gaussian kernels[END_REF] applied genetic programming to sentence embeddings to evolve Gaussian kernels for sentiment classification. Similarly to the case when word embeddings are used, genetic programming shows its ability to efficiently evolve programs that process the sentence embedding representation to solve the prediction tasks.

Also related to our research are proposals that apply embeddings for a better visualization and understanding of evolutionary search spaces. Michalak [START_REF] Michalak | Low-dimensional Euclidean embedding for visualization of search spaces in combinatorial optimization[END_REF] proposes to use a different class of embeddings, known as the t-Distributed Stochastic Neighbor Embedding (t-SNE). At last, Fyvie et al. [START_REF] Fyvie | Towards explainable metaheuristics: PCA for trajectory mining in evolutionary algorithms[END_REF] apply principal components analysis for trajectory mining of estimation of distribution algorithms. They used this approach to identify new methods of population diversity. To our knowledge, neural embeddings have not been previously used as an analytical or visualization tool of search spaces in evolutionary computation.

Embedding Representations for Other Domains

While neural embeddings have been mainly proposed for text representation, an increasing number of works report representations learned for other domains, including biological sequences [START_REF] Asgari | Continuous distributed representation of biological sequences for deep proteomics and genomics[END_REF][START_REF] Bepler | Learning protein sequence embeddings using information from structure[END_REF][START_REF] Yang | Learned protein embeddings for machine learning[END_REF], program source codes [START_REF] Alon | code2vec: Learning distributed representations of code[END_REF] and graphs [START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF]. Particularly relevant for our work is research on learning graph embeddings. Some of these approaches construct the graph node descriptions by first conducting a random walk on the graph [START_REF] Grover | node2vec: Scalable feature learning for networks[END_REF][START_REF] Perozzi | Deepwalk: Online learning of social representations[END_REF]. In this paper, we take inspiration from these previous approaches to create sequences of search trajectories from which we construct the corpora required to learn the embeddings.

INTRODUCING NEURAL EMBEDDING REPRESENTATIONS OF NK LANDSCAPES

As discussed in Section 2.2, word embeddings are learned from a corpus of sentences comprising words from a vocabulary. The embedding captures the semantic relationships between the words as expressed in the sequences. In our approach, we want to represent the relationship between solutions from an NK landscape (either local optima or not), in terms of their fitness values. Therefore, our vocabulary will be integrated by NK solutions; i.e., each word from our vocabulary is a binary solution from the search space. We define a sentence as a sequence of NK model solutions that are related according to some predefined semantics. In this paper, the chosen semantics establishes that each solution in the sequence, except the first one, is the best improving neighbor for the previous solution in the sequence. The neighborhood structure N can be defined in different ways, we here consider the usual 1-bit flip neighborhood structure; i.e., two solutions are neighbors if their Hamming distance is 1. Sequences of solutions can be constructed by means of adaptive walks [START_REF] Kauffman | Origins of Order[END_REF][START_REF] Klemm | Geometry and coarse-grained representations of landscapes[END_REF][START_REF] Nowak | Analysis of adaptive walks on NK fitness landscapes with different interaction schemes[END_REF]. For small search spaces, an alternative is to conduct an exhaustive enumeration and to run an adaptive walk from each possible solution of the search space, as done for constructing LONs in [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF]. For the purpose of our analysis, we follow this exhaustive approach in the paper.

For each solution from the binary search space 0 ∈ , we generate the sequence ( 0 , 1 , . . . , ℓ ) such that ∈ N ( -1 ) and ( ) > ( -1 ) for all ∈ {1, . . . , ℓ }. We follow a best improvement adaptive walk such that = arg max ∈N ( -1 ) ( ), ∈ {1, . . . , ℓ }, until no further improvement is possible. The adaptive walk determines the local optimum ℓ , and therefore the basins of attraction of 0 , required for constructing the sequences for inferring the embedding space. Notice that, when there is no equivalent solution in the neighborhood, there is one single local optima per starting solution for such an adaptive walk. Figure 2 shows an example of the way two sentences are constructed for a toy NK landscape with 6 variables. In Figure 2, is the order of the words in the sentence, the binary representation of the solution, ( ) the index of the binary solution in the vocabulary, and ( ) the NK model fitness value.

The embedding representation is learned from the set of sequences by means of gradient optimization of the neural network loss function. Since the neural network is shallow, e.g., the network in Figure 1 has only one hidden layer, the learning process is extremely fast and is scalable to datasets comprising of millions of sentences. In this paper, we rely on the word2vec algorithm [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF], that represents sentences as a sequence of indices of the words in the vocabulary. This means that the model does not actually use the particular encoding associated with each index, being a word for word embeddings or a binary vector for NK embeddings. What is relevant for the algorithm is the way in which these components are organized in the sequence.

Once the neural network is trained, the learned embeddings are extracted as explained in Section 2.2 and will serve as a basis to scrutinize the characteristics of the considered NK landscape, as captured by the embeddings. More details on the methods for learning word embeddings using the word2vec method can be found in [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF][START_REF] Mikolov | Distributed representations of words and phrases and their compositionality[END_REF].

Parameters for Learning NK Embeddings

An important parameter to be specified before using the word2vec algorithm is the size of the embedding representation. In this paper we consider to study the impact of using different embedding dimensions. Furthermore, the word2vec algorithm assumes that for a word to be represented as an embedding, it has to appear a minimum number of times in the corpus (5 times is the default parameter). This represents an issue in our case since some NK model solutions might only appear in a single sequence. To sort out this obstacle, we simply include five copies of the complete set of sequences in the input dataset. This guarantees that each solution is represented by one embedding, without biasing the embedding representation. The incurred computational cost can be neglected since the embedding space is only learned once.

In addition to the number of times each sentence is included in the corpus, the way we generate sequences of NK solutions determine that they may have a different length. For instance, a sequence that starts from a local optimum has a single element. Since this sort of short sequences does not allow the algorithm to learn any pattern from the sequence, we simply complete short sequences to a minimum length of 5 elements. Completion is done by repeating the last element from the original sequence (i.e., the local optimum) until reaching the required minimum length.

Another relevant parameter for the word2vec algorithm is the window size. We set this parameter to = 1, and we use the skipgram variant of word2vec; i.e., the task to be solved by the neural network is to predict the two words that surround every single word in a sequence.

Evaluating the Information Content of NK Embeddings

The usability of an embedding space lies in its capacity to encode the particular characteristics of the NK landscape. We propose to compute a number of metrics in order to evaluate the information from the NK model that the embedding encodes. In this section, we introduce the metrics and the rationale behind using them.

All the metrics are based on the information about the distance between solutions in the embedding space, given that this distance provides an idea of the "landscape" of the NK model in this space. Commonly, in the analysis of embedding spaces, the cosine distance between embeddings is used. We computed both the Euclidean and the cosine distances in our analysis. However, due to space restriction, we focus on the cosine distance below, as emphasized by the use of in the metrics described below. The considered metrics are as follows:

(1) nsteps_vs_cos_last: Correlation between the number of steps (from 0 to ℓ ) performed by adaptive walks to reach the local optimum (i.e.; the value of ℓ) and the (cosine) distance between 0 and ℓ in the embedding space. (2) fitness_diff_vs_cos_last: Correlation between the fitness difference among each starting and ending points of the sequence (i.e.; ( ℓ ) -( 0)) and the cosine distance between them in the embedding space. (3) max_fitness_diff_vs_cos_last: Correlation between the fitness difference ( ( ★ )-( 0)) of each starting point 0 with the global optimum ★ and the cosine distance between them in the embedding space. (4) basin_size_vs_cos_last: Correlation between the size of the basin of the attraction of 0 and the cosine distance between them in the embedding space.

For all considered metrics, we use the Spearman correlation coefficient.

The nsteps_vs_cos_last metric serves to determine whether local optima that require more steps to be reached from the starting point are similarly more distant from the starting solution in the embedding space. The same rationale is behind the fitness_diff_vs _cos_last metric, but in this case we consider the actual difference between the fitness values. The max_fitness_diff_vs_cos_last metric is more related to the original fitness-distance correlation measure of problem difficulty [START_REF] Jones | Fitness distance correlation as a measure of problem difficulty for genetic algorithms[END_REF]. A high correlation would indicate that, as we move closer to the global optimum in the embedding space, the fitness value also increases. At last, the basin_size_vs_cos _last metric serves to determine whether there is any relationship between the size of a local optimum's basin of the attraction and the distance of each point to its local optimum. It addresses the following question: Is a point closer to its corresponding local optimum if the latter has a larger basin of attraction?

EXPERIMENTAL ANALYSIS

In this section, we report experimental results addressing the following research questions:

• How is the landscape of solutions in the embedding space?

Are solutions with similar fitness clustered together, or is the distance among solutions in the embedding space determined by their corresponding Hamming distance? • What is the distribution of local optima in the embedding space? • Is the landscape in the embedding space smoother than in the original binary representation? • Which is the influence of in the characteristics of the embedding space? • Can optimization be conducted in the space of the embedding representation by mapping the continuous representation to the fitness values of encoded solutions?

We start by introducing the considered benchmark. Thus, a visual exploration of the learned embedding spaces is presented. Subsequently, we provide an analysis of the metrics introduced in Section 4.2. Finally, we report results on the application of a numerical optimizer in the (continuous) embedding space.

Benchmark

In terms of benchmark problems, we generate NK landscapes with = 14, ∈ {1, . . . , 12}, and random variables interactions. Fitness tables are generated following a uniform distribution in [0, 1]. For each parameter setting, 50 instances are independently generated. For each of the 2 14 solutions in a given NK landscape, we construct the path from each solution to its corresponding local optimum using a best-improvement adaptive walk. To evaluate the effect of the embedding dimension, embedding spaces are generated for each set of sequences for embeddings of dimension ∈ {2, 3, 5, 10, 25, 50, 100}. Therefore, the total number of embedding spaces considered in our experiments is 12 × 50 × 7 = 420.

Visual Inspections of the Embedding Spaces

Figure 3 gives a representation of the obtained embedding space for different values and an embedding dimension of = 2 for selected NK landscapes. Local optima are represented as circles, and other solutions as squares. The color of each solution renders its fitness values. We observe that solutions are distributed in a boomerang-shaped area, where local optima are concentrated to the extremes of the boomerang arms. As the number of local optima increases with , their location gets closer to the base of the boomerang. The shapes also display different orientations for different NK models. However, for a specific instance, the shapes are similar under multiple repetitions of the embedding learning algorithm (data not shown due to space constraints). Another important characteristic of the embedding space is that solutions are not uniformly distributed, and there are large areas where no solution is represented.

Using the embedded representations of the complete space of 2 solutions, we can map the NK fitness function to the continuous space. However, there are large areas of the continuous space that are empty. In order to "fill the gaps", we propose to simply assign any -dimensional vector for which a pre-image in the (binary) domain does not exist with the fitness of the closest point in the embedding space for which a pre-image exists. This amounts to a reconstruction or estimation of the NK model in the continuous space. We use an efficient implementation of the nearest neighbor searching procedure based on KDtrees [START_REF] Maneewongvatana | It's okay to be skinny, if your friends are fat[END_REF] to determine the fitness value of any point within the boundaries of the embedding space. Figure 4 shows such reconstructions for the same embedding spaces displayed in Figure 3. For the reconstructions, a grid of 120 × 120 points in the continuous space is used. Displayed in the figure as a black dot is the location of the global optimum. Remarkably, for the three different values, the global optimum is located in the yellow area that covers most local optima.

Metric-based Analysis

In the next set of experiments, we evaluate the metrics defined in Section 4.2 for increasing values and embeddings of different dimensions. Figure 5 reports the correlation between the cosine distance, in the embedding space, from an initial point to its local optimum and: (left) the number of steps performed by the adaptive walk to reach the local optimum, (center) the fitness difference between the starting point and its corresponding local optimum, (right) the fitness difference between the starting point and the global optimum.

We remark that the strongest correlation arises between the number of steps and the distance. The more the intermediate solutions visited by the adaptive walk to reach the local optimum in the original space, the higher the distance from the initial solution to the local optimum in the embedding space. This correlation increases with but it does differently depending on the dimensionality of the embeddings. The fitness difference between a solution and its local optimum is, for any and values, always positively correlated with the distance between the two solutions in the embedding space. However, the strength of the correlation depends on , and there is an important gap in the strength of these correlations for embeddings of dimension = 2 and those of higher dimensions. These results indicate that the distribution of solutions according to the fitness, as reported in Figure 3, significantly changes in higher dimensions.

The fitness-distance correlation was one of the first metrics proposed for fitness landscape analysis [START_REF] Jones | Fitness distance correlation as a measure of problem difficulty for genetic algorithms[END_REF]. We report the fitnessdistance correlation using the distance in the embedding space in Figure 5 (right). It can be seen in the figure that the correlation is always negative but the strength of the correlation depends both on and . For a high dimensionality the embedding space is highly deceptive, being closer to the global optimum does not implies an improvement in the fitness. It is important to notice that the correlations are computed considering all solutions. Another approach could consider solutions that are within a certain distance to the global optimum

The correlation between the distance and the size of local optima's basins of attraction is reported in Figure 6. It is weaker than the previously-discussed correlations. Similarly, it increases with , although it does not seem to depend on the embedding dimensionality. Regardless of the dimension of the embeddings, the distances of the points to their corresponding local optima is ranked in a similar way with respect to the size of the basin of attraction. We analyze below the relationship between the average distance of each solution to its [START_REF] Goldberg | word2vec explained: Deriving Mikolov et al. 's negativesampling word-embedding method[END_REF] neighbors _ ℎ , and the average distance of each solution to its local optimum _ . Figure 7 shows this information for different values and ∈ {2, 10, 100}. Analyzing in detail the distance to the neighbors is important since neighboring solutions are the set of candidates for the local search algorithm. While in the original domain the distance to neighbors is fixed to 1, it is not clear how this distance is in the embedding space. It can be seen in Figure 7 that, for small values, a solution gets closer to its neighbors in the embedding space as it is further away from the local optima. As increases, and the number of local optima also increases, the average distance of a solution to its neighbors increases relatively to the distance to its local optimum. This effect is visible for all the embedding dimensions. However, the transition is smoother for = 2 and much more abrupt for = 100.

Optimizing in the Embedding Space

An open question is how the difficulty of optimizing the NK model in the original binary space translates to a search in the embedding space, and whether this difficulty is also associated to the ruggedness of the original landscape, as defined by . To look into this issue, we consider the surrogate function ˜ ( ) that assigns to each possible continuous solution , the fitness value of the closest point in the embedding space with a pre-image in the corresponding NK landscape. To optimize this continuous function, we use a variant of the covariance matrix adaptation evolution strategy (CMA-ES) [START_REF] Hansen | Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation[END_REF], with an increasing population size restart [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF] (with an increasing factor of = 2), and the covariance matrix constrained to be diagonal [START_REF] Ros | A simple modification in CMA-ES achieving linear time and space complexity[END_REF]. This configuration was selected for efficiency and after preliminary experiments. The stopping condition for CMA-ES was to reach a maximum number of 500 function evaluations.

We evaluate the relative fitness deviation of the solution returned by CMA-ES to the global optimum; i.e.,

˜ ( ★ )-( ˆ ) ( ★ )
where ˆ and ★ are the optimum, and the best solution found by CMA-ES, respectively. As a second measure of performance, we compute the average rank of the solution found by CMA-ES (in terms of its fitness value), where the optimum has a rank of 0 and the solution with the lowest fitness has rank of 2 14 -1. CMA-ES was executed 30 times for each of the 50 instances corresponding to models with ∈ {1, . . . , 12}, and for embedding dimension ∈ {2, 3, 5, 10}. In total, 30 × 50 × 12 × 4 = 72 000 runs of CMA-ES were performed.

For the experiments, we used the Optuna library [START_REF] Akiba | Optuna: A next-generation hyperparameter optimization framework[END_REF].

Figure 8 reports that, as increases, the relative deviation to the optimum of the solution found by CMA-ES also increases. This effect is more noticeable for = 2 and less clear for higher dimensions ≥ 5. This result suggests that embedding spaces indeed capture the increase of difficulty induced by higher values, although this effect also depends on the embedding dimensionality. Another relevant finding from Figure 8 is that, as the embedding dimensionality increases, so does the approximation quality of the solution found by the optimizer. We speculate this could be due to the fact that, in higher dimensions, it is more likely that the basin of attractions of two different solutions are communicating. In two dimensions, the optimal solution could be isolated in a region completely surrounded by lower-quality solutions.

Figure 9 reports the average rank of the solution found by CMA-ES as increases. This figure provides a different perspective of the algorithm behavior. Instead of increasing, the average rank of the solution decreases with until it reaches a plateau, at different values for different . This indicates that, even if CMA-ES is not able reach the optimal solution, it often identifies high-quality local optima.

Discussion

We summarize the findings from our experimental analysis by revising the questions addressed at the beginning of this section.

Firstly, in embedding landscapes, the neighborhood relationship determined by 1-bit flip does not hold. The solutions that are closest to a solution in the embedded representation do not necessarily include the neighbors in the original binary domain. The location of solutions in the embedding space is more related to whether solutions are local optima, or not. This organization is consistent with the intuition behind the construction of word embeddings, where words with a similar "semantics" are grouped closer in the embedding space. In our model, the primary semantics relates to local optimality.

Secondly, for two dimensions, we corroborated that solutions are distributed in a boomerang-shaped form where local optima are more likely to be located in the arms. This distribution changes as the proportion of local optima increases with . Although it is difficult to determine whether the embedding landscape is smoother than in the original domain, we have seen that the fitness-distance correlation is strongly negative, and that it achieves the highest values for embeddings of high dimensionality. We have also shown that it is indeed possible to conduct the optimization process in the embedding space and observed that, for medium and large values, the average rank of the best-found solutions is not affected by the ruggedness of the original landscape.

CONCLUSIONS

In this paper, we proposed to approach NK landscapes from the perspective of neural embeddings. We addressed the question of how to create embeddings for NK models, and showed that these representations have as a characteristic feature that non-local solutions and local optima are located in different areas of the embedding space. We also corroborated that the distance between solutions in the embedding space is well correlated with the fitness difference between any two solutions and with the number of steps to reach local optima.

The transformation of the original search space to the continuous domain opens new perspectives that are worth being explored. We showed that it is possible to transform the search for the optimal solution from a binary to a continuous search space. This continuous space could be used to learn a model on the continuous representation that serves as a surrogate, or to infer properties of some NK solutions based on the proximity they have in the embedding space. One important direction for future research is to evaluate embedding spaces constructed using only a subset of solutions. In principle, it would be possible to apply similar methods for creating a set of paths from which learning the embedding spaces. However, a requirement would be to have a good covering of solutions from different areas of the search space. Another direction that deserves consideration is how to learn an inverse mapping from the continuous representation to the binary domain (e.g., using a neural network or a simple regression approach). Such inverse mapping would allow us to evaluate if perturbations to high-fitness solutions in the embedding representation can be mapped back to other high-fitness solutions in the binary representation. Finally, the strategies for the construction of the embedding spaces could be applied to other discrete optimization problems.

Figure 1 :

 1 Figure 1: Skipgram model used to learn the embeddings.

Figure 3 :

 3 Figure 3: Embedding spaces for exemplary NK landscapes with different values. The size of the local optima is proportional to their basin size.

Figure 4 :

 4 Figure 4: Estimation of the fitness landscape in the continuous domain from the NK embeddings. The location of the optimal solution is indicated with a black circle.

Figure 5 :

 5 Figure 5: Correlation between the distance, in the embedding space, from an initial point to its local optimum and: (left) the number of steps performed by the adaptive walk to reach the local optimum, (center) the fitness difference between the starting point and the local optimum, (right) the fitness difference between the starting point annd the global optimum.

Figure 6 :

 6 Figure 6: Correlation between the distance, in the embedding space, from an initial point to the size of the local optimum's basin of attraction.

Figure 7 :

 7 Figure 7: Relationship between the average distance of a point to its neighbors and the distance of a point to its local optimum for different values of and ∈ {2, 10, 100}.

Figure 8 :

 8 Figure 8: Relative distance to the optimal fitness value for the different embedding sizes as increases.

Figure 9 :

 9 Figure 9: Rank, in terms of the fitness value, of the best found solution for the different embedding sizes as increases. Rank 0 corresponds to the optimum. Ranks are normalized by the search space size (2 14 ).

ACKNOWLEDGMENTR.

  Santana acknowledges support by the Basque Government, Spain (project KK-2020/00049 through the ELKARTEK program and IT 124419), and by the Spanish Ministry of Economy and Competitiveness MINECO, Spain (projects TIN2016-78365-R and PID2019-104966GB-I00).