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

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

Hernán Aguirre

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

What Makes an Instance Difficult for Black-Box 0-1 Evolutionary Multiobjective Optimizers?

This paper investigates the correlation between the characteristics extracted from the problem instance and the performance of a simple evolutionary multiobjective optimization algorithm. First, a number of features are identified and measured on a large set of enumerable multiobjective NK-landscapes with objective correlation. A correlation analysis is conducted between those attributes, including low-level features extracted from the problem input data as well as high-level features extracted from the Pareto set, the Pareto graph and the fitness landscape. Second, we experimentally analyze the (estimated) running time of the global SEMO algorithm to identify a (1 + ε)-approximation of the Pareto set. By putting this performance measure in relation with problem instance features, we are able to explain the difficulties encountered by the algorithm with respect to the main instance characteristics.

Introduction

In single-objective black-box combinatorial optimization, fitness landscape analysis aims at apprehending the relation between the geometry of a problem instance and the dynamics of randomized search algorithms. Understanding the main problem-related features allows to explain the behavior and the performance of such algorithms, the ultimate goal being to predict this performance and adapt the algorithm setting to the instance being solved. Recently, the performance of single-objective randomized search algorithms has been correlated to fitness landscape features [START_REF] Daolio | Local optima networks and the performance of iterated local search[END_REF]. In this paper, we propose a general methodology to analyze the correlation between problem features and algorithm performance in black-box 0-1 evolutionary multiobjective optimization. To the best of our knowledge, this is the first time that such an analysis is conducted in multiobjective optimization.

We first identify a number of existing and original multiobjective problem features. They include low-level features extracted from the problem input data like variable correlation, objective correlation, and objective space dimension, as well as high-level features from the Pareto set, the Pareto graph and the ruggedness and multimodality of the fitness landscape. Some of them are here proposed for the first time. They consist of a simple autocorrelation function, based on a local hypervolume measure, and allowing to estimate the ruggedness of the fitness landscape. We report all these measures on a large number of enumerable multiobjective NK-landscapes with objective correlation (ρMNKlandscapes), together with a correlation analysis between them.

Next, we conduct an experimental analysis on the correlation between instance features and algorithm performance. To do so, we investigate the estimated running time of a simple evolutionary multiobjective optimization algorithm, namely global SEMO [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF], to identify a (1 + ε)-approximation of the Pareto set. In particular, the original hypervolume-based autocorrelation functions appear to be the features with the highest correlation with the algorithm performance. Overall, the running time of the algorithm is impacted by each of the identified multiobjective problem feature. Our analysis shows their relative importance on the algorithm efficiency. Moreover, taking the features all together allows to better explain the dynamics of randomized search algorithms.

The paper is organized as follows. Section 2 details the background information related to fitness landscape analysis, multiobjective optimization and ρMNK-landscapes. In Section 3, low-level and high-level instance features are identified, and quantitative results, together with a correlation analysis, are reported for ρMNK-landscapes. Section 4 presents the experimental setup of global SEMO and discusses the correlation between the problem features and the estimated running time of global SEMO. Section 5 concludes the paper and discusses further research.

Preliminaries

Fitness Landscape Analysis

In single-objective optimization, fitness landscape analysis allows to study the topology of a combinatorial optimization problem [START_REF] Weinberger | Correlated and uncorrelatated fitness landscapes and how to tell the difference[END_REF], by gathering important information such as ruggedness or multimodality. A fitness landscape is defined by a triplet (X, N , φ), where X is a set of admissible solutions (the search space), N : X → 2 X is a neighborhood relation, and φ : X → IR is a (scalar) fitness function, here assumed to be maximized. A walk over the fitness landscape is an ordered sequence x 0 , x 1 , . . . , x of solutions from the search space such that x 0 ∈ X, and x t ∈ N (x t-1 ) for all t ∈ {1, . . . , }.

An adaptive walk is a walk such that for all t ∈ {1, . . . , }, φ(x t ) > φ(x t-1 ), as performed by a conventional hill-climbing algorithm. The number of iterations, or steps, of the hill-climbing algorithm is the length of the adaptive walk. This length is a good estimator of the average diameter of the local optima basins of attraction, characterizing a problem instance multimodality. The larger the length, the larger the basin diameter. This allows to estimate the number of local optima when the whole search space cannot be enumerated exhaustively.

Let x 0 , x 1 , . . . be an infinite random walk over the search space. The autocorrelation function and the correlation length of such a random walk allow to measure the ruggedness of a fitness landscape [START_REF] Weinberger | Correlated and uncorrelatated fitness landscapes and how to tell the difference[END_REF]. The random walk autocorrelation function r : N → IR of a (scalar) fitness function φ is defined as follows.

r(k) = E[φ(x t ) • φ(x t+k )] -E[φ(x t )] • E[φ(x t+k )] Var(φ(x t )) (1)
where E[φ(x t )] and Var(φ(x t )) are the expected value and the variance of φ(x t ), respectively. The autocorrelation coefficients r(k) can be estimated within a finite random walk x 0 , x 1 , . . . , x of length .

r(k) = -k t=1 (φ(x t ) -φ) • (φ(x t+k ) -φ) t=1 (φ(x t ) -φ) 2 (2)
where φ = 1 t=1 φ(x t ), and 0. The estimation error diminishes with the walk length . The correlation length τ measures how the autocorrelation function decreases. This characterizes the ruggedness of the landscape: the larger the correlation length, the smoother the landscape. Following [START_REF] Weinberger | Correlated and uncorrelatated fitness landscapes and how to tell the difference[END_REF], we define the correlation length by τ = -1 ln(r(1)) , making the assumption that the autocorrelation function decreases exponentially.

Multiobjective Optimization

A multiobjective optimization problem can be defined by an objective vector function f = (f 1 , . . . , f M ) with M 2 objective functions, and a set X of feasible solutions in the decision space. In the combinatorial case, X is a discrete set. Let Z = f (X) ⊆ IR M be the set of feasible outcome vectors in the objective space. To each solution x ∈ X is assigned an objective vector z ∈ Z on the basis of the vector function f : X → Z with z = f (x). The conventional Pareto dominance relation is defined as follows. In a maximization context, an objective vector z ∈ Z is dominated by an objective vector z ∈ Z, denoted by z ≺ z , if and only if ∀m ∈ {1, . . . , M }, z m z m and ∃m ∈ {1, . . . , M } such that z m < z m . By extension, a solution x ∈ X is dominated by a solution x ∈ X, denoted by x ≺ x , if and only if f (x) ≺ f (x ). A solution x ∈ X is said to be Pareto optimal (or efficient, non-dominated), if and only if there does not exist any other solution x ∈ X such that x ≺ x. The set of all Pareto optimal solutions is called the Pareto set X ⊆ X. Its mapping in the objective space is called the Pareto front Z ⊆ Z. One of the most challenging task in multiobjective optimization is to identify a minimal complete Pareto set [START_REF] Ehrgott | Multicriteria optimization[END_REF], i.e. a Pareto set of minimal size, that is one Pareto optimal solution for each point from the Pareto front.

However, in the combinatorial case, generating a complete Pareto set is often infeasible for two main reasons [START_REF] Ehrgott | Multicriteria optimization[END_REF]: (i) the number of Pareto optimal solutions is typically exponential in the size of the problem instance, and (ii) deciding if a feasible solution belongs to the Pareto set may be NP-complete. Therefore, the overall goal is often to identify a good Pareto set approximation. To this end, heuristics in general, and evolutionary algorithms in particular, have received a growing interest since the late eighties.

ρMNK-Landscapes

The family of ρMNK-landscapes constitutes a problem-independent model used for constructing multiobjective multimodal landscapes with objective correlation [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. It extends single-objective NK-landscapes [START_REF] Kauffman | The Origins of Order[END_REF] and multiobjective NKlandscapes with independent objective functions [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF]. Feasible solutions are binary strings of size N , i.e. the decision space is X = {0, 1} N . The parameter N refers to the problem size (the bit-string length), and the parameter K to the number of variables that influence a particular position from the bit-string (the epistatic interactions). The objective vector function f = (f 1 , . . . , f m , . . . , f M ) is defined as f : {0, 1} N → [0, 1) M . Each objective function f m is to be maximized and can be formalized as follows.

f m (x) = 1 N N i=1 c m i (x i , x i1 , . . . , x i K ) , m ∈ {1, . . . , M } (3)
where c m i : {0, 1} K+1 → [0, 1) defines the multidimensional component function associated with each variable x i , i ∈ {1, . . . , N }, and where K < N . By increasing the number of variable interactions K from 0 to (N -1), ρMNK-landscapes can be gradually tuned from smooth to rugged. In this work, we set the position of these epistatic interactions uniformly at random. The same epistatic degree K m = K and the same epistatic interactions are used for all objectives m ∈ {1, . . . , M }. Component values are uniformly distributed in the range [0, 1), and follow a multivariate uniform distribution of dimension M , defined by a correlation coefficient ρ > -1

M -1 , i.e. the same correlation ρ is defined between all pairs of objective functions. As a consequence, it is very unlikely that the same objective vector is assigned to two different solutions. The positive (respectively negative) data correlation allows to decrease (respectively increases) the degree of conflict between the objective function values very precisely [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. An instance generator and the problem instances under study in this paper can be found at the following URL: http://mocobench.sf.net/.

In the following, we investigate ρMNK-landscapes with an epistatic degree K ∈ {2, 4, 6, 8, 10}, an objective space dimension M ∈ {2, 3, 5}, and an objective correlation ρ ∈ {-0.9, -0.7, -0.4, -0.2, 0.0, 0.2, 0.4, 0.7, 0.9} such that ρ > -1

M -1 . The problem size is set to N = 18 in order to enumerate the search space exhaustively. The search space size is then |X| = 2 18 . 30 different landscapes, independently generated at random, are considered for each parameter combination: ρ, M , and K. This leads to a total of 3300 problem instances.

Problem Features and Correlation Analysis

In this section, we identify a number of general-purpose features, either directly extracted from the problem instance itself (low-level features), or computed from the enumerated Pareto set and from the fitness landscape (high-level features). Then, a correlation analysis is conducted on those features in order to highlight the main similarities in characterizing the difficulties of a problem instance.

Low-level Features from Problem Input Data

First, we consider the following low-level features related to the definition of ρMNK-landscapes.

Number of epistatic interactions (K):

This gives the number of variable correlations in the construction of ρMNK-landscapes. As will be detailed later, despite the K-value can generally not be retrieved directly from an unknown instance, it can be precisely estimated within some high-level fitness landscape metrics described below. Number of objective functions (M): This parameter represents the dimension of the objective space in the construction of ρMNK-landscapes. Objective correlation (ρ): This parameter allows to tune the correlation between the objective function values in ρMNK-landscapes. In our analysis, the objective correlation is the same between all pairs of objectives.

High-level Features from the Pareto Set

The high-level fitness landscape metrics considered in our analysis are described below. We start with some general features related to the Pareto set.

Number of Pareto optimal solutions (npo):

The number of Pareto optimal solutions enumerated in the instance under consideration simply corresponds to the cardinality of the (exact) Pareto set, i.e. npo = |X |. The approximation set manipulated by any EMO algorithm is directly related to the cardinality of the Pareto optimal set. For ρMNK-landscapes, the number of Pareto optimal solutions typically grows exponentially with the problem size, the number of objectives and with the degree of conflict between the objectives [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. Hypervolume (hv): The hypervolume value of a the Pareto set X gives the portion of the objective space that is dominated by X [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. We take the origin as a reference point z = (0.0, . . . , 0.0). Average distance between Pareto optimal solutions (avgd): This metric corresponds to the average distance, in terms of Hamming distance, between any pair of Pareto optimal solutions. Maximum distance between Pareto optimal solutions (maxd): This metric is the maximum distance between two Pareto optimal solutions in terms of Hamming distance.

High-level Features from the Pareto Graph

In the following, we describe some high-level features related to the connectedness of the Pareto set [START_REF] Ehrgott | Connectedness of efficient solutions in multiple criteria combinatorial optimization[END_REF][START_REF] Gorski | Connectedness of efficient solutions in multiple objective combinatorial optimization[END_REF]. If all Pareto optimal solutions are connected with respect to a given neighborhood structure, the Pareto set is said to be connected, and local search algorithms would be able to identify many non-dominated solutions by starting with at least one Pareto optimal solution; see e.g. [START_REF] Paquete | A two-phase heuristic for the biobjective 0/1 knapsack problem[END_REF][START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF]. We follow the definition of k-Pareto graph from [START_REF] Paquete | A two-phase heuristic for the biobjective 0/1 knapsack problem[END_REF]. The k-Pareto graph is defined as a graph P G k = (V, E), where the set of vertices V contains all Pareto optimal solutions, and there is an edge e ij ∈ E between two nodes i and j if and only if the shortest distance between solutions x i and x j ∈ X is below a bound k, i.e. d(x i , x j ) k. The distance d(x i , x j ) is taken as the Hamming distance for ρMNK-landscapes. This corresponds to the bit-flip neighborhood operator. Some connectedness-related high-level features under investigation are given below. to the smallest distance k such that the k-Pareto graph is connected, i.e. for all pairs of vertices x i , x j ∈ V in P G k , there exists and edge e ij ∈ E.

High-level Features from the Fitness Landscape

At last, we give some high-level metrics related to the number of local optima, the length of adaptive walks, and the autocorrelation functions.

Number of Pareto local optima (nplo):

A solution x ∈ X is a Pareto local optimum with respect to a neighborhood structure N if there does not exist any neighboring solution x ∈ N (x) such that x ≺ x ; see e.g. [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF].

For ρMNK-landscapes, the neighborhood structure is taken as the 1-bit-flip, which is directly related to a Hamming distance 1. This metric reports the number of Pareto local optima enumerated on the ρMNK-landscape under consideration. Length of a Pareto-based adaptive walk (ladapt): We here compute the length of adaptive walks by means of a very basic single solution-based Pareto-based Hill-Climbing (PHC) algorithm. The PHC algorithm is initialized with a random solution. At each iteration, the current solution is replaced by a random dominating neighboring solution. As a consequence, PHC stops on a Pareto local optimum. The number of iterations, or steps, of the PHC algorithm is the length of the Pareto-based adaptive walk. As in the single-objective case, the number of Pareto local optima is expected to increase exponentially when the adaptive length decreases for ρMNKlandscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. Correlation length of solution hypervolume (corhv): The ruggedness is here measured in terms of the autocorrelation of the hypervolume along a random walk. As explained in Section 2.1, the correlation length τ measures how the autocorrelation function, estimated with a random walk, decreases. The autocorrelation coefficients are here computed with the following scalar fitness function φ : X → IR: φ(x) = hv({x}), where hv({x}) is the hypervolume of solution x ∈ X, the reference point being set to the origin. The random walk length is set to = 10 4 , and the neighborhood operator is the 1-bit-flip.

Correlation length of local hypervolume (corlhv):

This metric is similar to the previous one, except that the fitness function is here based on a local hypervolume measure. The local hypervolume is the portion of the objective space covered by non-dominated neighboring solutions, i.e. for all x ∈ X, φ(x) = hv(N (x) ∪ {x}). Similarly to corhv, the random walk length is set to = 10 4 , and the neighborhood operator N is the 1-bit-flip.

Correlation Analysis

The correlation matrix between each pair of features is reported in Fig. First of all, when taken independently, the number of objective functions M and the objective correlation ρ are both moderately correlated to the cardinality of the Pareto set npo (the absolute correlation coefficient is around 0.5 in both cases). Surprisingly, the objective space dimension does not explain by itself the large amount of nondominated solutions found in many-objective optimization [START_REF] Wagner | Pareto-, aggregation-, and indicator-based methods in many-objective optimization[END_REF]. As pointed out in [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF], this should be put in relation with the degree of conflicts between the objective function values. Indeed, as shown in Fig. 2, it is easy to build a simple multilinear regression model based on M and ρ to predict the value of npo with a very high precision (resulting in a correlation coefficient of 0.87, and explaining 76% of the variance). This highlights that the impact of many-objective fitness landscapes on the search process cannot be analyzed properly without taking the objective correlation into account.

Interestingly, other important remarks can be extracted from the figure. With respect to the Pareto set, the hypervolume value increases with the objective space dimension. Moreover, and unsurprisingly, the Pareto set size and the size of the largest connected component from the Pareto graph are highly correlated. So are the maximum distance between Pareto optimal solutions and the minimum distance for the Pareto set to be connected. As also reported in [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF], there is a high correlation between the number of Pareto optimal solutions npo and of Pareto local optima nplo. More importantly, the number of Pareto local optima nplo can be precisely estimated with the length of a Pareto-based adaptive walk ladapt (the absolute correlation coefficient between log(nplo) and ladapt is 1). As a consequence, this allows to estimate the size of the Pareto set as well. At last, the number of epistatic interactions (decision variable correlations) K can be estimated with hypervolume-based autocorrelation functions along a random walk corhv and corlhv. Since there is not much difference between the correlations coefficients of both functions, the first one corhv should preferably be considered due to its simplicity.

4 Problem Features vs. Algorithm Performance

Experimental Setup

Global SEMO. Global SEMO (G-SEMO for short) [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF] is a simple elitist steadystate EMO algorithm for black-box 0-1 optimization problems dealing with an arbitrary objective vector function defined as f : {0, 1}

N → Z such that Z ⊆ IR M , like ρMNK-landscapes. A pseudo-code is given in Algorithm 1. It maintains an unbounded archive A of non-dominated solutions found so far.

Algorithm 1 Pseudo-code of G-SEMO Input: x 0 ∈ X Output: Archive A 1: A ← x 0 2: loop 3: select x from A at random 4:
create x by flipping each bit of x with a probability 1/N 5:

A ← non-dominated solutions from A ∪ {x } 6: end loop

The archive is initialized with one random solution from the search space. At each iteration, one solution is chosen at random from the archive. Each bit of this solution is independently flipped with a rate r = 1/N , and the obtained solution is checked for insertion in the archive. Within such an independent bitflip mutation, any solution from the search space can potentially be reached by applying the mutation operator to any arbitrary solution. In its general form, the G-SEMO algorithm does not have any explicit stopping rule [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF]. In this paper, we are interested in its running time, in terms of a number of function evaluations, until an (1 + ε)-approximation of the Pareto set has been identified and is contained in the internal memory A of the algorithm, subject to a maximum number of function evaluations.

Performance Evaluation. For any constant value ε 0, the (multiplicative) ε-dominance relation ε can be defined as follows. For all z, z ∈ Z, z ε z if and only if z m • (1 + ε) z m , ∀m ∈ {1, . . . , M }. Similarly, for all x, x ∈ X, x ε x if and only if f (x) ε f (x ). Let ε 0. A set X ε ⊆ X is an (1 + ε)approximation of the Pareto set if and only if, for any solution x ∈ X, there is one solution x ∈ X ε such that x ε x . This is equivalent of finding a Pareto set approximation whose multiplicative epsilon quality indicator value with respect to the exact Pareto set is (1 + ε), see e.g. [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. Interestingly, under some general assumptions, there always exists an (1 + ε)-approximation, for any given ε 0, whose cardinality is both polynomial in the problem size and in 1/ε [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF].

Following a conventional methodology from single-objective continuous blackbox optimization benchmarking [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009[END_REF], the expected number of function evaluations to identify an (1 + ε)-approximation is chosen as a performance measure. However, as any EMO algorithm, G-SEMO can either succeed or fail to reach an accuracy of ε in a single simulation run. In case of a success, the runtime is the number of function evaluations until an (1 + ε)-approximation was found. In case of a failure, we simply restart the algorithm at random. We then obtain a "simulated runtime" [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009[END_REF] from a set of given trials of G-SEMO on a given instance. Such a performance measure allows to take into account both the success rate p s ∈ (0, 1] and the convergence speed of the G-SEMO algorithm. Indeed, after (n -1) failures, each one requiring T f evaluations, and the final successful run with T s evaluations, the total runtime is T = n-1 i=1 T f + T s . By taking the expectation value and by considering that the probability of success after (n -1) failures follows a Bernoulli distribution of parameter p s , we have:

E[T ] = 1 -p s p s E[T f ] + E[T s ] (4)
In our case, the success rate p s is estimated with the ratio of successful runs over the total number of executions ( ps ), the expected runtime for unsuccessful runs E[T f ] is set to a constant function evaluation limit T max , and the expected runtime for successful runs E[T s ] is estimated with the average number of function evaluations performed by successful runs.

ert = 1 -ps ps T max + 1 N s Ns i=1 T i (5)
where N s is the number of successful runs, and T i is the number of evaluations required for successful run i. For more details, we refer to [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB-2009[END_REF].

Parameter Setting. In our analysis, we set ε = 0.1. The time limit is set to T max = 2 N /10 < 26215 function evaluations without identifying an (1 + ε)approximation. The G-SEMO algorithm is executed 100 times per instance. For a given instance, the success rate and the expected number of evaluations for successful runs are estimated from those 100 executions. However, let us note that G-SEMO was not able to identify a (1 + ε)-approximation set for any of the runs on one instance with M = 3, ρ = 0.2 and K = 10, one instance with M = 3, ρ = 0.4 and K = 10, ten instances with M = 5, ρ = 0.2 and K = 10, six instances with M = 5, ρ = 0.4 and K = 10, as well as two instances with M = 5, ρ = 0.7 and K = 10. Moreover, G-SEMO was not able to solve the following instances due to an overload CPU resources available: M = 5 and ρ ∈ {-0.2, 0.0}. Overall, this represents a total amount of 2980 instances times 100 executions, that is 298000 simulation runs.

Computational Results

The correlation between each feature and the running time of G-SEMO is reported in Fig. 3. First, with respect to low-level features, there exists a high correlation between log(ert) and K, which is the highest absolute correlation observed on our data. However, surprisingly, the correlation of the performance measure with M and ρ is not significant. Second, with respect to high-level features from the Pareto set, the size of the Pareto set and its hypervolume does not explain the variance of log(ert). Nevertheless, the larger the distance between Pareto optimal solutions in the decision space, the larger the running time of G-SEMO. Similarly, when the Pareto graph is close to a fully connected graph, G-SEMO is likely to take less time to identify a (1 + ε)-approximation (the absolute correlation value is around 0.3). As a consequence, the number of Pareto Fig. 3. Correlation between log(ert) and each feature. The feature names are reported on the first line, correlation coefficients are reported on the second line, and scatter plots as well as smoothing splines are displayed on the third line.

optimal solutions has a smaller impact on the performance of G-SEMO than the structure existing between those solutions in the decision space.

With respect to high-level fitness landscape features, the number of Pareto local optima nplo and its estimator ladapt both present a significant correlation with the estimated running time of G-SEMO. Indeed, the more Pareto local optima, the longer the running time (the absolute correlation value is close to 0.5). At last, the hypervolume-based autocorrelation functions highly explain the variance of the G-SEMO performance. For both corhv and corlhv, the absolute correlation value is around 0.8. Overall, this correlation analysis gives a "big picture" of a well-suited multiobjective fitness landscape for G-SEMO. This corroborates the impact of the problem instance properties identified in the previous section on the performance of multiobjective evolutionary algorithms.

Discussion

In this paper, we attempted to give a first step towards a better understanding of the evolutionary multiobjective optimization algorithm performance according to the main characteristics of the problem instance. We first presented a number of general problem features, together with a correlation analysis between those features on a large set of enumerable multiobjective NK-landscapes. Then, we put in relation the running time of a simple evolutionary multiobjective optimization algorithm with those features. Our analysis clearly shows the high impact of theses problem-related properties on the performance of the algorithm. In particular, two relevant hypervolume-based autocorrelation functions have been proposed for the first time, allowing to precisely estimate the ruggedness of the instance under consideration, as well as the algorithm running time.

Using the general methodology introduced in the paper applied to larger problem instances would allow to appreciate the impact of the multiobjective features on the performance of evolutionary multiobjective optimizations when tackling large-size instances. This should be possible with features that do not require the complete enumeration of the decision space, including the problem size, the number of objectives, the objective correlation, the length of a Pareto-based adaptive walk, and the hypervolume-based autocorrelation functions proposed in this paper. As well, the impact of the stopping condition, and in particular the approximation quality (the ε-value) should be carefully investigated. At last, a similar study would allow to better understand the structure of the landscape for real-world multiobjective combinatorial optimization problems. This work pushes towards the design of a meta-algorithm able to select the most efficient evolutionary multiobjective algorithm or parameter setting according to a prediction model based on the main problem instance features.

  Number of connected components (nconnec): This metric gives the number of connected components in the 1-Pareto graph, i.e. in P G k with k = 1. Size of the largest connected component (lconnec): This corresponds to the size of the largest connected component in the 1-Pareto graph P G 1 . Minimum distance to be connected (kconnec): This measure corresponds

  1. Correlation matrix between all pairs of features. The feature names are reported on the diagonal. For each pair of features, scatter plots and smoothing splines are displayed below the diagonal, and the corresponding correlation coefficients are reported above the diagonal.

												10
												8
											log(npo)	6
												4
												2
												0
												1	2	3	4	5	6	7
												(beta_1 x M) + (beta_2 x rho) + e
											Fig. 2.	Scatter	plot	of
											the linear regression model
											log(npo) = β1M +β2ρ+e, with
			2.0 4.5		0 6		2 8		0 800		2 8 β1 = 1.30567, β2 = -2.87688, 0 4 8 0.2 0.7
											and e = 0.27735. Residual
	2.0 4.5	K	-9.9e-22 M	0.00 0.26	-0.082 0.57	-0.22 0.90	0.60 0.13	0.56 -0.23	0.32 0.39	-0.35 0.53	0.45 standard error: 1.037 on 3297 0.31 -0.34 -1.00 -0.97 2 6 -0.34 0.35 -0.39 -0.025 0.085 degrees of freedom, multiple R-squared: 0.7629, adjusted R-squared: 0.7627, F-statistic:
				ρ	-0.49	0.06	-0.25	0.17	-0.22	-0.35	0.24 5303 on 2 and 3297 DF, -0.67 0.64 0.0085 0.12 -0.5 p-value: < 2.2e -16.
	0 6				log_npo	0.67	0.28	-0.40	0.59	0.90 -0.57 0.81 -0.81 0.065	0.075
						hv	0.04	-0.41	0.32	0.64 -0.50	0.41	-0.44	0.19	0.25
												-0.7
	2 8						avgd	0.58	0.24	0.021	0.49	0.56	-0.58 -0.60 -0.59
								maxd	-0.17	-0.49 0.89 -0.086	0.068	-0.55	-0.55	2 8
	0 800								nconnec	0.44	-0.28	0.48	-0.48 -0.33	-0.29
										log_lconnec -0.64 0.61 -0.59	0.34	0.33	0 6
	2 8										kconnec	0.25	-0.44 -0.44
												log_nplo -1.00 -0.33	-0.34	4 8
	0 4 8											ladapt	0.36 corhv	0.36 0.97	0.4 0.8
	0.2 0.7
		2 6		-0.5	-0.7		2 8		0 6		4 8	0.4 0.8

corlhv Fig. 1.


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

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

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

Hernán Aguirre

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

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

Decomposition-based Multi-objective Landscape Features and Automated Algorithm Selection

Landscape analysis is of fundamental interest for improving our understanding on the behavior of evolutionary search, and for developing general-purpose automated solvers based on techniques from statistics and machine learning. In this paper, we push a step towards the development of a landscape-aware approach by proposing a set of landscape features for multi-objective combinatorial optimization, by decomposing the original multi-objective problem into a set of single-objective sub-problems. Based on a comprehensive set of bi-objective ρmnk-landscapes and three variants of the state-of-the-art Moea/d algorithm, we study the association between the proposed features, the global properties of the considered landscapes, and algorithm performance. We also show that decomposition-based features can be integrated into an automated approach for predicting algorithm performance and selecting the most accurate one on blind instances. In particular, our study reveals that such a landscape-aware approach is substantially better than the single best solver computed over the three considered Moea/d variants.

Introduction

Context. Evolutionary algorithms have been proven extremely effective for solving a broad range of optimization problems. In the last decades, the community has gained a deep understanding on the key components underlying the design of a successful evolutionary approach for a given problem. However, one of the main challenge remains to automate the process of choosing the most suitable algorithm or configuration for the instance under consideration. In fact, it is well known that the structural properties of an optimization problem highly impact the dynamics and performance of search algorithms, leading to the requirement of adopting a landscape-aware algorithm selection and configuration methodology for the success of evolutionary problem solving. On the one hand, landscape analysis [START_REF]Recent Advances in the Theory and Application of Fitness Landscapes[END_REF] provides a principled approach for studying and analyzing the relation between the underlying search space structure and algorithm behavior. On the other hand, machine learning techniques can be leveraged to perform sophisticated tasks, such as predicting algorithm performance, identifying the best algorithm configuration, or selecting the best algorithm [START_REF] Kerschke | Automated algorithm selection: Survey and perspectives[END_REF]. As a byproduct, and since the pioneer work of Rice [START_REF] Rice | The algorithm selection problem[END_REF], landscape-aware algorithm configuration has emerged as an appealing approach for increasing the effectiveness and efficiency of evolutionary algorithms. In this paper, we contribute to the development of such an approach when specifically dealing with multi-objective and combinatorial optimization problems.

Related Work. Independently of whether the target problems and algorithms aim at optimizing a single or multiple objectives, and of whether they have a combinatorial or continuous nature, every landscape-aware methodology needs to address the following two research challenges: (i) the design of a set of informative and high-level landscape features, and (ii) the development of automated recommendation systems integrating the so-designed features on the basis of statistical or machine learning prediction models. Looking at the specialized literature, a large amount of work has been conducted for single-objective optimization since pioneering works in the field [START_REF] Grefenstette | Predictive models using fitness distributions of genetic operators[END_REF]. For instance, in single-objective continuous optimization, the exploratory landscape analysis (ELA) constitutes one major achievement made by the community to collect and combine existing features under a common tool and methodology [START_REF] Kerschke | Automated algorithm selection: Survey and perspectives[END_REF]. Similarly, a number of features for single-objective combinatorial optimization have been developed over the years, and recent studies integrate them into sophisticated automated approaches for algorithm selection and configuration [START_REF] Beham | Algorithm selection on generalized quadratic assignment problem landscapes[END_REF][START_REF] Liefooghe | Towards landscapeaware automatic algorithm configuration: Preliminary experiments on neutral and rugged landscapes[END_REF]. Those single-objective features are either based on problem-specific characteristics such as the maximum cost between two cities in the traveling salesperson problem [START_REF] Mersmann | A novel feature-based approach to characterize algorithm performance for the traveling salesperson problem[END_REF], or on general descriptors from the underlying landscape. In the latter case, this is achieved by relying on a neighborhood relation over the search space in order to define a (combinatorial) landscape, and by studying its properties and characteristics in terms of multimodality, ruggedness, or neutrality [START_REF]Recent Advances in the Theory and Application of Fitness Landscapes[END_REF].

Despite the significant progress made in the last decades, the existing literature on the development of a unified landscape-aware approach targeting multiobjective optimization problems is more scarce [START_REF] Kerschke | The r-package flacco for exploratory landscape analysis with applications to multi-objective optimization problems[END_REF][START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multi-objective optimization[END_REF]. Although landscape features can in principle be applied to the multi-objective case, the statistical and machine learning models considered in the single-objective case are still to be studied and validated when turning into a multi-objective setting. For multi-objective combinatorial optimization, one can refer to the recent study in [START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multi-objective optimization[END_REF], providing a comprehensive analysis of problem-independent multi-objective landscape features, and showing their effectiveness in predicting algorithm performance and in selecting from an algorithm portfolio. The features described there are mostly based on dominance and (hypervolume) indicator. In fact, they were designed to grasp the landscape characteristics, but also to capture the search behavior of dominance-based multi-objective algorithms. In a subsequent study [START_REF] Liefooghe | Dominance, indicator and decomposition based search for multi-objective QAP: landscape analysis and automated algorithm selection[END_REF], those previous features were shown to provide useful information about the performance of different classes of multi-objective evolutionary algorithms (MOEAs).

However, their correlation with the performance of decomposition-based MOEAs reveal to be less significant [START_REF] Liefooghe | Dominance, indicator and decomposition based search for multi-objective QAP: landscape analysis and automated algorithm selection[END_REF]. Getting inspiration from the fact that successful MOEAs may rely on different search paradigms, our work considers applying other mechanisms, such as decomposition, to design new multi-objective landscape features, hence pushing one step further the development of a landscapeaware automated evolutionary approach.

Methodology and Contribution. We rely on the concept of decomposition [START_REF] Trivedi | A survey of multiobjective evolutionary algorithms based on decomposition[END_REF] to develop a new set of general-purpose multi-objective landscape features and to study their effectiveness when integrated into an automated algorithm selection task. Our interest in the concept of decomposition stems from the fact that it provides a state-of-the-art framework, represented by the Moea/d algorithm [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF], by simply decomposing the multi-objective problem into a number of single-objective sub-problems. We view the decomposition paradigm as an opportunity to leverage existing single-objective features for multi-objective landscape analysis. Our contributions can be summarized as follows.

(i) We propose a new set of high-level multi-objective landscape features based on decomposition. Intuitively, we attempt to capture the multi-objective landscape by aggregating the characteristics from the single-objective subproblem landscapes obtained by decomposition. As such, the proposed features are obtained by first defining single-objective features for each subproblem, and then aggregating them by means of descriptive statistics. (ii) We consider the task of predicting performances of three variants of the stateof-the-art Moea/d algorithm using a tree-based regression model to study the effectiveness of a decomposition-based landscape-aware methodology for automatically selecting the best performing algorithm for a given instance. (iii) Throughout an extensive set of experiments using ρmnk-landscapes with two objectives as a case study, we conduct a systematic analysis on the association between the designed features and the benchmarked landscapes, as well as the association between features and algorithm performance. Our findings reveal that the designed features are able to capture the benchmark parameters (ρ and k), and to substantially improve the so-called single best solver when integrated into a landscape-aware algorithm selection approach.

Outline. In Section 2, we describe the proposed multi-objective landscape features. In Section 3, we study the association of features with benchmark parameters. In Section 4, we investigate the integration of the proposed features into an automated landscape-aware approach. In Section 5, we conclude the paper.

2 From Single-to Multi-objective Features based on Decomposition

Multi-objective Optimization

A multi-objective combinatorial optimization problem can be defined by a set of m objective functions f = (f 1 , f 2 , . . . , f m ), and a discrete set X of feasible solutions in the decision space. Let Z = f (X) ⊆ IR m be the set of feasible outcome vectors in the objective space. To each solution x ∈ X is assigned an objective vector z ∈ Z, on the basis of the vector function f : X → Z. In a maximization context, an objective vector z ∈ Z is dominated by a vector z ∈ Z iff ∀m ∈ {1, . . . , m}, z m z m and ∃m ∈ {1, . . . , m} s.t. z m < z m . A solution x ∈ X is dominated by a solution x ∈ X iff f (x) is dominated by f (x ). A solution x ∈ X is Pareto optimal if there does not exist any other solution x ∈ X such that x is dominated by x. The set of all Pareto optimal solutions is the Pareto set. Its mapping in the objective space is the Pareto front. Our goal is to identify a good Pareto set approximation, for which multi-objective evolutionary algorithms (MOEAs) constitute a popular effective option [START_REF] Deb | Multi-Objective Optimization using Evolutionary Algorithms[END_REF].

Rationale, Methodology and Features Overview

In this work, we get inspiration from the so-called MOEAs based on decomposition [START_REF] Trivedi | A survey of multiobjective evolutionary algorithms based on decomposition[END_REF] in order to design new high-level multi-objective features. This algorithm class is based on searching for good-performing solutions in multiple regions of the Pareto front by decomposing the original multi-objective problem into a number of scalarized single-objective sub-problems. Each sub-problem is obtained by a different parameterization of the same underlying scalarizing function. This is typically what the state-of-the-art Moea/d algorithm [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF] performs, while introducing a cooperation among sub-problem solving. In particular, this offers much flexibility for integrating existing single-objective search operators and solvers, which is actually one of the main reasons for the success of decomposition-based MOEAs. Let us however recall that our main goal is not to design a new multi-objective algorithm, but to design new multi-objective features that can feed the design of a general-purpose landscape-aware approach.

Therefore, we propose to rely on the simple observation that each of the so-defined sub-problems also implies a single-objective landscape that we can attempt to analyze and characterize. In other words, by studying the singleobjective landscape implied by the sub-problems, we should be able to extract some knowledge about the original multi-objective problem. More precisely, the methodology that we adopt in the reminder of this paper consists in: (i) defining a number of single-objective landscapes using decomposition, (ii) extracting single-objective features for each sub-problem landscape, and (iii) aggregating those single-objective features into new multi-objective features. These steps are detailed below.

Sub-problem Landscape Definition. Firstly, we define µ scalarized singleobjective sub-problems, where both the scalarizing function and the µ value are user-defined parameters. Among the different scalarizing functions that may be used, the Chebyshev function is one of the most effective since it can be shown that any Pareto optimal solution can be achieved by solving a well-parameterized Chebyshev sub-problem. In the rest of this paper, we should hence use the Chebyshev scalarizing function:

g(x|ω) = max i∈{1,...,m} ω i • z i -f i (x) , such that x ∈ X, ω = (ω 1 , . . . , ω m
) is a positive weight vector, and z = (z 1 , . . . , z m ) is a reference point such that z i > f i (x) ∀x ∈ X, i ∈ {1, . . . , m}. It should be clear that µ different sub-problems can be obtained by choosing µ different weight vectors, denoted by ω j , j ∈ {1, . . . , µ}.

It is worth noticing that we do not make any assumption about the original (black-box) multi-objective problem, so that we have no information about what region every sub-problem is actually mapping to. Hence, the value of µ as well as the procedure to generate the weight vector can be a critical issue. This is studied in more details later when reporting our empirical results.

Single-objective Landscape Features. Next, we define a landscape for every single-objective sub-problem j ∈ {1, . . . , µ}, for which we compute a number of underlying high-level single-objective features. Following the standard literature on single-objective landscape analysis [START_REF]Recent Advances in the Theory and Application of Fitness Landscapes[END_REF], the landscape of sub-problem j can be defined as a triplet (X, N , g(•|ω j )), such that N : X -→ 2 X is a neighborhood relation defined among solutions for the considered problem; e.g., 1-bit-flips for binary strings, or swaps for permutations.

The considered sub-problem features are based on sampling the so-defined landscape to compute some statistics. Following the standard literature [START_REF]Recent Advances in the Theory and Application of Fitness Landscapes[END_REF][START_REF] Kauffman | The Origins of Order[END_REF], we consider two simple sampling strategies, namely random walks (rws) and adaptive walks (aws). Generally speaking, a walk is an ordered sequence of solutions (x 0 , x 1 , . . . , x ) such that x 0 ∈ X, and x t ∈ N (x t-1 ) ∀t ∈ {1, . . . , }. In a random walk, x t being the current solution, the next solution x t+1 is simply a random neighbor under N . The length of a random walk is a user-defined parameter. In an adaptive walk, the next solution x t+1 is selected to be an improving neighbor with respect to the single-objective scalarizing function g(•|ω j ). Consequently, the length of an adaptive walk is the number of steps performed until no further improvement is possible, x is then a local optimum. Notice that the reference point z * required for computing the scalar fitness values is updated on the basis of the best objective values seen so far during the walk.

Given a sub-problem j ∈ {1, . . . , µ} and a walk (x 0 , x 1 , . . . , x ), we consider the following four classes of single-objective features, as summarized in Table 1: -Fitness value (fv *). In the first class, we compute some statistics informing about the distribution of fitness values observed along the walk. More precisely, we consider the mean (avg) and standard deviation (sd) of the fitness values of collected solutions. We also consider three additional statistics, namely the first auto-correlation coefficient (r1), the kurtosis (kr), and the skewness (sk) of fitness values. The kurtosis and the skewness are standard measures in statistical analysis, while the first auto-correlation coefficient is mostly used in the landscape analysis literature. Denoting by ḡ(•|ω j ) the average fitness value of solutions in the walk, the first auto-correlation coefficient is defined as follows [START_REF] Kauffman | The Origins of Order[END_REF]: -Fitness difference (fd *). In the second class, we compute the average fitness difference with the neighboring solutions for every x i , i ∈ {1, . . . , }:

r1 = -1 t=0 g(x t |ω j ) -ḡ(•|ω j ) • g(x t+1 |ω j ) -ḡ(•|ω j ) -1 t=0 g(x t |ω j ) -ḡ(•|ω j ) 2
1 |N (xi)| x∈N (xi) g(x i |ω j ) -g(x|ω j
) . Similarly, we consider the mean (avg), standard deviation (sd), minimum (min) and maximum (max) fitness difference over solutions from the walk.

-Improving neighbors (in *). In the third class, we compute the proportional number of improving neighbors for each solution x i , i ∈ {0, . . . , }.

Then, we consider the mean, standard deviation, minimum, and maximum of this measure. It is worth noticing that the second and third classes of features require to evaluate the fitness value of neighbors from all solutions from the walk. -Length of the adaptive walk (law *). The fourth class only contains features extracted from the adaptive walk. In particular, we consider the length of the adaptive walk as a feature to characterize the sub-problem landscape. This length was shown to provide an estimation of the number of local optima in single-objective landscape analysis [START_REF] Kauffman | The Origins of Order[END_REF].

Aggregated Multi-objective Features. The features described above are computed for each sub-problem j ∈ {1, . . . , µ}, and then have a dimension µ.

For µ > 1, we need to aggregate these µ-dimensional single-objective features into 1-dimensional multi-objective features. To do so, we use two standard statistics, namely the mean (avg) and the standard deviation (sd). In addition, we use a polynomial regression in order to fit each single-objective feature as a function of the weight vector ω j of sub-problem j. The coefficient of the polynomial model are then used as additional aggregated features. In this study, since we experiment bi-objective optimization problems, we consider a second order polynomial regression and propose to use the first (c1) and the second (c2) coefficients as additional multi-objective features. As summarized in Table 1, we end up with 4 aggregation statistics over respectively 5, 4, and 4 fv, fd, and in features, in addition to 4 aggregated features on the length of adaptive walk. This amounts to a total of 108 decomposition-based multi-objective features.

A Preliminary Exploratory Analysis

As a first step, we analyze the relevance of the proposed features in capturing the characteristics of multi-objective optimization problems regardless of any particular evolutionary search algorithm. Therefore, we conduct a preliminary exploratory analysis in order to highlight the association between the designed features and the properties of well-established benchmark landscapes.

Experimental Setup

Following previous work [START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multi-objective optimization[END_REF], we consider ρmnk-landscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF] as a problemindependent model used for constructing multi-objective multimodal landscapes with objective correlation. Candidate solutions are binary strings of size n. The objective function vector f = (f 1 , . . . , f i , . . . , f m ) is defined as f : {0, 1} n → [0, 1] m such that each objective f i is to be maximized. The objective value f i (x) of a solution x = (x 1 , . . . , x j , . . . , x n ) is an average value of the individual contributions associated with each variable x j . Given objective f i , and each variable x j , a component function f ij : {0, 1} k+1 → [0, 1] assigns a real-valued contribution for every combination of x j and its k epistatic interactions {x j1 , . . . , x jk }. These f ij -values are uniformly distributed in [0, 1]. The objective functions to be maximized can written as:

f i (x) = 1 n n j=1 f ij (x j , x j1 , .
. . , x jk ), ∀i ∈ {1, . . . , m}. In this work, the k epistatic interactions are set uniformly at random among the (n -1) variables other than x j . By increasing the value of k from 0 to (n -1), problem instances can be gradually tuned from smooth to rugged. The f ij -values additionally follow a multivariate uniform distribution of dimension m, defined by an m × m positive-definite symmetric covariance matrix (c pq ) s.t. c pp = 1 and c pq = ρ for all p = q where ρ > -1 m-1 defines the correlation among the objectives. In our work, we focus on bi-objective landscapes, i.e., m = 2. We use a latin hypercube sampling to generate a set of 1 000 balanced instances spanning parameters ranges: n ∈ {50, 51, . . . , 200}, k ∈ {0, 1, 2, . . . , 8} and ρ ∈] -1, 1]. The random walk length is set to = 1 000 across all problem sizes. A unique random walk is performed for all sub-problems, whereas one adaptive walk is performed for each sub-problem. The number of sub-problems is set to µ = 20 and weight vectors are distributed uniformly; i.e., ω j = ((j -1)/(µ -1), 1 -(j -1)/(µ -1)). The neighborhood relation is the standard 1-bit-flip.

For our analysis, we first conduct an exploratory analysis to better visualize and understand the proposed features. Next, we construct a regression model to study the accuracy of features in grasping the global properties of ρmnk-landscapes, and we analyze the correlation among features.

Visual Analysis of Single-objective Features

In Fig. 1, we report the values of single-objective features as a function of subproblems. Due to space restriction, we report a single representative feature for each of the four classes. The blue curves correspond to ρmnk-landscapes with k = 0, the green ones to k = 2 and the red ones to k = 4. The color scales from red to orange, green to cyan, and blue to purple, respectively, and correspond to the objective correlation parameter ρ varying from 1 to -1; i.e., from highly correlated to highly conflicting objectives. For example, the standard deviation of fitness values from a random walk (fv rws sd, bottom left), the average fitness difference from a random walk (fd rws avg, second line, second column) and the standard deviation of improving neighbors (in rdw sd, second line, third column) gives a clear differentiation between landscapes with different k-values. The lower the benchmark parameter k, the lower the standard deviation of fitness values and the average fitness difference. Similarly, the standard deviation of fitness values from an adaptive walk (vf aws sd, top left) and the length of an adaptive walk (law aws avg, top right) seems to be clearly associated with parameter ρ. The flatter the curve rendering the evolution of these two features as a function of weight vectors, the higher the objective correlation ρ.

From our visual inspection, we can conclude that the landscape features are representative of the different global benchmark parameters, which are unknown in a black-box optimization scenario. However, this first analysis considers the single-objective features and not the aggregated 1-dimensional multi-objective features, which are discussed next.

Correlation Analysis of Features and Landscape Parameters

Investigating the accuracy of the designed multi-objective features, we consider a typical machine learning task consisting in predicting the value of the (unknown) global benchmark parameters k and ρ. We respectively construct a random forest classification model and a random forest regression model [START_REF] Breiman | Random forests[END_REF], using the whole set of multi-objective features computed over all considered ρmnk-landscapes. Random forest has the nice property of providing a measure of feature importance for model fitting. In Fig. 2, we report the relative importance of each feature extracted from the random forest models, using the Gini impurity as a measure of quality. Values are averaged over 10 independent repetitions of model fitting.

The first notable observation is that feature importance is different depending on whether we aim at predicting benchmark parameter ρ or k. The objective correlation ρ appears mostly related to a single feature: the standard deviation of the number of improving neighbors (in rws sd). By contrast, deciding on parameter k is related to multiple different features, that mostly correspond to the fitness difference computed from adaptive walks, in particular the maximum fitness difference (fd aws max). To push the analysis further, we show in Fig. 3 the Spearman correlation matrix among a subset of features as well as benchmark parameters ρ and k. It is worth noticing that we do not show all 108 features, but only a subset of 60 features with representative correlation values due to space restriction. Interestingly, we can distinguish four clusters, denoted C 1 , C 2 , C 3 and C 4 in the figure. The largest cluster C 1 contains more than 30 features. There is a positive correlation between the first auto-correlation coefficient of fitness values vf r1 and k, of at least 0.25. This cluster also contains a small subset of 2 features related to the number of improving neighbors in *, with a particularly high correlation value (> 0.8) (the 3 × 3 red points at the center of the figure). In the second cluster C 2 , we observe a large subset of features being negatively correlated with k (the dark blue columns/lines intersecting in a completely red square in the middle-right of the figure). For both clusters related to k (C 1 and C 2 ), features are not extracted from a unique class. This means that features from different classes can be used to characterize k.

The remaining 2 clusters (C 3 and C 4 ) are associated with the second benchmark parameter ρ. In particular, besides the average number of improving neighbors (in rws avg), cluster C 3 also contains the standard deviation (in rws sd), both computed from a random walk. This latter feature, that was found to be the most important to predict ρ in Fig. 2, has a nearly perfect negative correlation with ρ, equals to -1. The last cluster C 4 contains features based on the fitness values computed over adaptive walks (fv aws), which appear to be the most positively correlated with ρ (3 × 3 red points at the bottom right of the figure), with a correlation higher than 0.7.

To summarize, we find that all feature classes are useful to characterize the (unknown) global benchmark parameter k, rendering the degree of non-linearity of the problem. However, the fitness distribution (fv) and the number of improving neighbors (in) classes have more importance than the fitness difference (fd) and the length of adaptive walks (law) for characterizing the (unknown) global benchmark parameter ρ, relating to objective correlation.

Landscape-aware MOEA/D Selection

In this section, we conduct a second set of experiments in order to study the accuracy of the designed features when integrated into an automated algorithm selection approach. We consider the more sophisticated task of selecting the best performing algorithm among an algorithm portfolio. More precisely, we consider three variants of the well-established Moea/d algorithm [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF][START_REF] Marquet | Shake them all! -rethinking selection and replacement in MOEA/D[END_REF] as a case study.

In the following, we start by briefly describing the considered portfolio before addressing our main target task.

Algorithm Portfolio

As mentioned earlier, the Moea/d algorithm is based on a flexible decompositionbased framework that can be configured in different manners. In its baseline variant [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF], Moea/d first decomposes the problem into a number of scalarized sub-problem, as discussed previously. Then, a solution is evolved for each subproblem in a cooperative way. The algorithm iterates over the sub-problems and, at each iteration, an offspring is generated by means of crossover and mutation on the basis of parent solutions selected from the so-called T -neighborhood; i.e., the sub-problems corresponding to the T closest weights in the objective space. The new offspring can then replace any of the sub-problem solutions in the Tneighborhood of the current sub-problem. This corresponds to a standard evolutionary process, where selection and replacement are performed iteratively over sub-problems. In [START_REF] Marquet | Shake them all! -rethinking selection and replacement in MOEA/D[END_REF], it is shown that the selection and replacement underlying the standard Moea/d framework are key algorithm components that highly impact performance. Several generational variants are proposed therein, allowing to tune the selection and replacement underlying the Moea/d framework from fully cooperative (i.e., among all sub-problems) to fully selfish (i.e., independently of any other sub-problem).

Interestingly, it was found that no variant outperforms the other independently of the global benchmark parameter values ρ and k for the considered ρmnk-landscapes. Since such parameters are unknown in a black-box optimization scenario, the study presented in [START_REF] Marquet | Shake them all! -rethinking selection and replacement in MOEA/D[END_REF] leaves open the challenging question of which variant to choose in an automated manner. In addition, this constitutes a perfect and typical setting for the main automated algorithm selection task addressed in this paper. We consider, in the following, three representative variants of the Moea/d framework, exposing different degrees of cooperation among sub-problem solving. For reproducibility, and in order to be consistent with the notations from [START_REF] Marquet | Shake them all! -rethinking selection and replacement in MOEA/D[END_REF], these variants are denoted as follows: (i) Moea/d, referring to the standard variant [START_REF] Zhang | MOEA/D: A multiobjective evolutionary algorithm based on decomposition[END_REF], (ii) Moea/d-sc, a generational variant where selection is performed in a selfish manner for every sub-problem whereas replacement is performed in a cooperative manner, and (iii) Moea/d-ss, a (selfish) generational variant exposing the lower degree of cooperation among sub-problems. Besides population size, the three variants have the same set of parameters: δ = 1, nr = 2, p mut = 1 n and p cr = 1. Due to space restriction, we refer to [START_REF] Marquet | Shake them all! -rethinking selection and replacement in MOEA/D[END_REF] for a full description of these Moea/d variants.

In order to highlight the relevance of this portfolio in studying the accuracy of our features when integrated into a landscape-aware algorithm selection approach, we briefly report their relative performance using exactly the same set of ρmnk-landscapes as in the previous section. Every algorithm is executed 20 times on each instance, using a population size equals to n, and a budget of 10 6 evaluations. The performance of an algorithm is computed as its hypervolume relative deviation w.r.t. the best-found approximation set for each instance. The Table 2. Performance matrix of the three Moea/d variants. The diagonal reports the number of times where the corresponding algorithm is statistically outperformed by another one (the lower the better). The other cells report how many instances (out of 1000) the algorithm in the corresponding line is statistically better than the algorithm in the corresponding column (the higher the better). hypervolume measures the area covered by an approximation set and enclosed by a reference point [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. For a given instance, the reference point is set to the best value seen across all runs for each objective. We then count the number of times an algorithm is statistically outperformed using a two-sided Mann-Whitney test with a p-value of 0.05 and a Bonferroni correction. Results are reported in Table 2.

Moea

We clearly see that each algorithm is outperformed by another one on a subset of instances. The basic Moea/d variant seems to have a reasonably good behavior, since it is less-often outperformed than the two other variants overall (see diagonal). A more detailed analysis, omitted due to space restriction, shows that there is a complex interaction between algorithm performance and the benchmark parameters ρ and kwhich can be summarized as follows: (i) the smaller k and ρ, the better Moea/d and Moea/d-sc against Moea/d-ss, (ii) the larger ρ (highly correlated), the better Moea/d-ss, and (iii) the larger k and the smaller ρ, the better Moea/d-sc. Of course, this general trend has some exceptions, but it shows the impact of the unknown benchmark parameters on the relative performance of algorithms.

Automated Algorithm Selection

Task and Experimental Methodology. We study the accuracy of the proposed features by investigating the selection of the best performing Moea/d variant. For this purpose, we adopt the following standard supervised-learning approach. We first train three models in order to predict the performance of every considered Moea/d variant. We use the average hypervolume deviation as defined in the previous section as a measure of algorithm performance on a given instance, which then corresponds to our output prediction variable. Considering an unseen test instance, the landscape features are first computed, the performance of each algorithm is then predicted on the basis of the trained models, and the algorithm having the best prediction is selected as the recommended one. We consider the same set of ρmnk-landscapes described in the previous section. We adopt a standard validation methodology where an instance is selected for training with probability 0.9 and for testing with probability 0.1. We use a set of 100 random regression trees to learn and predict the expected relative hypervolume deviation. Reported values are computed over 50 independent runs. We recall that the proposed multi-objective features rely on some weight vectors µ. We consider a variable number of weight vectors in the range µ ∈ {1, 2, 3, 4, 5, 6, 10, 20}. Additionally, we consider two alternatives for generating weight vectors, namely uniform or random. In a random setting, a weight vector is generated uniformly at random. In a uniform setting, the weight vector are evenly distributed in the objective space. In particular, for µ = 1, the weight vector selected in the uniform setting is (0.5, 0.5); i.e., the "middle" of the objective space. In this case, the multi-objective features are simply the same than the corresponding single-objective features from the single scalarized sub-problem. For µ = 2, our uniform setting corresponds to weight vector (1, 0) and (0, 1). This means that our features are obtained by aggregating the single-objective features computed independently for each objective of the original multi-objective problem. The impact of this setting is carefully analyzed in our experiments.

Prediction Accuracy. In order to assess the prediction accuracy, we compute three complementary measures. The first two directly relate to the prediction error: the percentage of times the selected algorithm does not have the best hypervolume deviation in average, and the percentage of times the selected algorithm is statistically outperformed by at least one other algorithm. The third indicator, which is a straightforward adaptation from [START_REF] Lindauer | The algorithm selection competition series 2015-17[END_REF], measures the gap between: (i) the performance of the single best solver (SBS) having the best performance in average over the training set (without model training), and (ii) the performance of the virtual best solver (VBS), obtained by a model that would make perfect predictions. More precisely, let I train and I test be the set of training and testing instances, respectively, and let rhf(A, i) be the average relative hypervolume deviation of a given algorithm A ∈ A = {Moea/d,Moea/d-ss,Moea/d-sc} on instance i ∈ I train ∪ I test . For every algorithm A ∈ A and instance subset J, let rhf(A, J) = 1 |J| i∈J rhf(A, i). We define SBS as the algorithm having the best rhf value on the training set I train , i.e., SBS = arg min A∈A rhf(A, I train ) . We define VBS as the virtual 'algorithm' obtained by a perfect prediction model (an oracle); i.e., the algorithm with the best rhf(•, i) value for each i ∈ I test . Finally, let Recommended Solver (RS) be the algorithm predicted by the actual trained model. The merit indicator is:

M = rhf(RS, I test ) -rhf(V BS, I test ) rhf(SBS, I test ) -rhf(V BS, I test )
It should be clear that: (i) A merit of 0 indicates that the model does not make any error, (ii) A merit in the range [0, 1[ indicates that the model is more efficient than the SBS but worse than the VBS, (iii) a merit greater than 1 indicates that the model is worse than the SBS. Achieving a merit value of 0 is clearly a very challenging task and one seeks for a merit value below 1 (better than the SBS) and as close as possible to 0 (the VBS).

Experimental Results and Discussion. Our main results are summarized in Fig. 4 showing the accuracy indicators as a function of the number of weight vectors µ and their type (random or uniform). For completeness, we also show the R 2 coefficient obtained by the training models. Different observations can be extracted from Fig. 4. We first clearly see that the choice of the weight vector distribution is of critical importance. In fact, a random choice does not obtain a good accuracy, except when the number of weights µ is substantially large. By contrast, a uniform strategy appears to perform reasonably well, even when the number of weights is low. Interestingly, for uniform weights, the worst accuracy is obtained with µ = 2. Notice that such a setting is even substantially outperformed by a random choice of weight vectors. This indicates that computing single-objective features independently for each objective is not a recommended strategy. By contrast, computing features for decomposed sub-problems is effective even when using a very low number of weights. This indicates that a decomposition-based approach for multi-objective landscape analysis contains a valuable information about algorithm performance. Surprisingly, we found that a uniform choice of few weight vectors with µ ∈ {1, 3} performs reasonably well, although increasing µ > 3 allows to obtain a better accuracy. The relatively good results achieved with µ = 1 are however to be interpreted very carefully, taking into account that the shape of the Pareto front for ρmnk-landscapes, although having different magnitude in the objective space, is convex, symmetric and centered in the middle of the objective space, regardless of the values of ρ and k [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. Although this is a recurrent observation for many multi-objective combinatorial optimization problems, one might need to carefully choose µ when tackling problems with different Pareto front shapes. Such considerations are left for future investigations.

At last, in order to further show the accuracy of the proposed multi-objective landscape features, we experiment a baseline random forest model using the (unknown) global benchmark parameters ρ and k as input variables to predict algorithm performance. Contrary to a black-box scenario where the knowledge about ρ and k is not available, the accuracy of such a 'white-box' model should highlight the relevance and reliability of the proposed black-box features. We found that such a model trained with ρ and k obtains an average merit of 0.41. Comparatively, black-box features obtain an average merit of 0.31, 0.37 and 0.29 respectively, for µ ∈ {1, 3, 20} uniform weight vectors. This once again indicates that the proposed approach is very effective, and that the designed high-level black-box features seem to provide more accurate prediction models than the global benchmark parameters, hence allowing to substantially improve over the single best solver, and also to get closer to the ideal virtual best solver.

Conclusion and Open Issues

In this paper, we push a step towards the development of automated landscapeaware selection and configuration approaches by proposing a set of multi-objective landscape features and analyzing their effectiveness in grasping the global properties of black-box multi-objective combinatorial optimization problems, together with their efficiency in selecting the best performing algorithm among three variants of the Moea/d state-of-the art algorithm. The proposed features are based on the simple idea of aggregating the single-objective features extracted from a number of sub-problems obtained by decomposition. Our empirical analysis on a wide range of ρmnk-landscapes provides insights into the accuracy of the proposed approach. However, it also raises some interesting research questions.

For instance, problems with more than two objectives, where the distribution of weight vector is expected to play an even more important role, as well as other real-like problems, are to be studied. Moreover, we excluded the cost of feature computation from our analysis, which can be an issue when the total affordable budget for the optimization task is restricted. An interesting observation is that we only need to extract features from a relatively low number of sub-problems, as considering a single sub-problem is already shown to provide a reasonably good performance. Additionally, solely one random walk is required to compute the considered single-objective features, so that one can eventually end up with negligible cost features. In this respect, a more fine-grained analysis of the cost-vs-importance of features is to be conducted in our future investigations. Finally, an interesting question is to conduct a systematic analy-sis of the proposed decomposition-based features w.r.t. existing dominance-and indicator-based features, and to analyze their relative cost and accuracy. Given that decomposition-based features were shown to be effective and does not need any (costly) dominance-and indicator-based computations, it is our hope that unifying the two classes of features would allow us to end up with efficient highlevel state-of-the-art features for landscape-aware multi-objective optimization.

Fig. 1 .

 1 Fig. 1. Feature values of ρmnk-landscapes decomposed into 20 sub-problems (n = 25, each color correspond to a particular configuration of ρ and k).

Fig. 2 .

 2 Fig. 2. Relative importance of features to predict ρ (left) and k (right).

Fig. 3 .

 3 Fig. 3. Pairwise correlation among a subset of features.

Fig. 4 .

 4 Fig. 4. R 2 (top left), merit (top right) and error rates (bottom) according to the number of weights µ and their distribution.

Table 1 .

 1 A summary of the proposed landscape features.

	description	sub-problem features random walk adaptive walk	MO features
	Fitness values	fv rws s s ∈ {avg, sd, r1, kr, sk} fv aws s	fv rws s r ; fv aws s r r ∈ {avg, sd, c1, c2}
	Fitness differences	fd rws s s ∈ {avg, sd, min, max} fd aws s	fd rws s r ; fd aws s r r ∈ {avg, sd, c1, c2}
	Improving neighbors	in rws s s ∈ {avg, sd, min, max} in aws s	in rws s r ; in aws s r r ∈ {avg, sd, c1, c2}
	Walk length	-	law	law r r ∈ {avg, sd, c1, c2}