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

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

Luís Paquete
email: paquete@dei.uc.pt

Jin-Kao Hao
email: jin-kao.hao@univ-angers.fr

Experiments on Local Search for Bi-objective Unconstrained Binary Quadratic Programming

This article reports an experimental analysis on stochastic local search for approximating the Pareto set of bi-objective unconstrained binary quadratic programming problems. First, we investigate two scalarizing strategies that iteratively identify a high-quality solution for a sequence of sub-problems. Each sub-problem is based on a static or adaptive definition of weighted-sum aggregation coefficients, and is addressed by means of a state-of-the-art single-objective tabu search procedure. Next, we design a Pareto local search that iteratively improves a set of solutions based on a neighborhood structure and on the Pareto dominance relation. At last, we hybridize both classes of algorithms by combining a scalarizing and a Pareto local search in a sequential way. A comprehensive experimental analysis reveals the high performance of the proposed approaches, which substantially improve upon previous best-known solutions. Moreover, the obtained results show the superiority of the hybrid algorithm over non-hybrid ones in terms of solution quality, while requiring a competitive computational cost. In addition, a number of structural properties of the problem instances allow us to explain the main difficulties that the different classes of local search algorithms have to face.

Introduction

The unconstrained binary quadratic programming (UBQP) problem is one of the most challenging problem from single-objective combinatorial optimization [START_REF] Kochenberger | The unconstrained binary quadratic programming problem: a survey[END_REF]. Given a collection of n items such that each pair of items is associated with a profit value that can be positive, negative or zero, the UBQP problem seeks a subset of items that maximizes the sum of their paired values. The value of a pair is summed up only if the two corresponding items are selected. From a computational point-of-view, a feasible solution to a UBQP instance can be represented as a binary string of size n. Each position from the binary string maps to a particular variable that indicates whether the corresponding item is included in the subset of selected items or not. Beyond its theoretical significance [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF], the utility of UBQP has been demonstrated on a wide variety of application fields [START_REF] Kochenberger | The unconstrained binary quadratic programming problem: a survey[END_REF]. Furthermore, a number of NP-hard combinatorial optimization problems can be recast as UBQP problems, such as graph coloring, max-cut, set packing, set partitioning, or maximum clique, among others [START_REF] Kochenberger | The unconstrained binary quadratic programming problem: a survey[END_REF]. The singleobjective UBQP problem has received a growing interest in recent years [START_REF] Glover | Efficient evaluations for solving large 0-1 unconstrained quadratic optimisation problems[END_REF][START_REF] Kochenberger | The unconstrained binary quadratic programming problem: a survey[END_REF], and a multi-objective extension of UBQP has been proposed recently [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF].

In this paper, we focus on bi-objective UBQP, where two profit values are associated with each pair of items. By optimizing both sums of profit values simultaneously, we can improve the descriptive power of the conventional singleobjective UBQP problem, and provide a more general formulation. However, as for many problems from multi-objective combinatorial optimization, the biobjective UBQP problem raises several difficulties for heuristics design. In particular, the number of optimal solutions can be very large [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF], and determining whether a candidate solution is optimal is NP-complete, even in the singleobjective case [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. For these reasons, we design and experiment with multiobjective stochastic local search algorithms, and measure their efficiency and their effectiveness on instances with different dimensions and correlation degrees between the objective function values. Furthermore, we analyze the problem structure to learn more about those difficulties, and to improve the design of algorithms. The contributions of the paper are two-fold.

(i) We characterize the features of small-size, enumerable bi-objective UBQP instances. More particularly, we analyze the number of global and local optimal solutions, based on scalarizing functions and on the Pareto dominance relation; and we examine the connectedness between optimal solutions. (ii) We design and analyze local search algorithms for bi-objective UBQP, including two scalarizing approaches, a Pareto-based approach, and a hybrid approach combining these two complementary search strategies. The designed algorithms substantially improve over the previous attempts in solving largesize bi-objective UBQP instances [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF]. More importantly, our experimental analysis allows us to better understand how the performance of these classes of algorithms relates to the structural properties of the search space, explaining the high efficiency and effectiveness of the proposed approaches.

The remainder of the paper is organized as follows. In Section 2, we present the bi-objective UBQP problem. In Section 3, we study the characteristics of small-size instances. In Section 4, we introduce four stochastic local search algorithms for bi-objective UBQP. In Section 5, we analyze the performance of these algorithms on a set of large-size bi-objective UBQP instances. In Section 6, we finally conclude the paper and discuss further research directions.

Bi-objective UBQP

This section presents the problem formulation, some definitions related to multiobjective combinatorial optimization, and the problem instances that are investigated in the paper.

Problem Formulation

The bi-objective UBQP (bUBQP) problem can be formalized as follows [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF].

max f 1 (x) = n i=1 n j=1 q ij 1 x i x j max f 2 (x) = n i=1 n j=1 q ij 2 x i x j subject to x ∈ {0, 1} n (1)
where (f 1 , f 2 ) is the pair of objective functions to be maximized, n is the number of items, Q 1 = (q ij 1 ) and Q 2 = (q ij 2 ) are both an n × n matrix of constant values, either positive, negative or zero. As in the single-objective case, the solution space X = {0, 1}

n is defined on binary strings of size n; its size is then 2 n .

Definitions

We denote by Z ⊆ IR 2 the feasible region in the objective space, i.e. the image of feasible solutions when using the maximizing function vector

f = (f 1 , f 2 ) such that Z = f (X). The Pareto dominance relation is defined as follows. A solution x ∈ X is dominated by a solution x ′ ∈ X, denoted as x ≺ x ′ , if f k (x) ≤ f k (x ′ )
for all k ∈ {1, 2}, with at least one strict inequality. If neither x ≺ x ′ nor x ′ ≺ x holds, then both solutions are mutually non-dominated. A solution x ∈ X is Pareto optimal if there does not exist any other solution

x ′ ∈ X such that x ≺ x ′ .
The set of all Pareto optimal solutions is the Pareto set, and its mapping in the objective space is the Pareto front. One of the most challenging issues in multiobjective optimization is to identify a minimal complete Pareto set, i.e. one Pareto optimal solution mapping to each point from the Pareto front. Since the bUBQP problem is both NP-hard and intractable [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF], approximate algorithms like stochastic local search are well suited to identify a Pareto set approximation.

Problem Instances

Following [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF], the definition of each bUBQP objective function is based on a matrix Q k , k ∈ {1, 2}. As in the single-objective UBQP instances available in the OR-lib [START_REF] Beasley | OR-library: Distributing test problems by electronic mail[END_REF], non-zero matrix integer values are randomly generated following a uniform distribution in [-100, +100]. The density d ∈ [0, 1] gives the expected proportion of non-zero entries in the matrix. Following a Bernoulli distribution of parameter d, a given entry at position (i, j) is set to zero on both matrices, i.e. q ij 1 = q ij 2 = 0. Moreover, we define a correlation coefficient ρ between the data contained in the two matrices. The positive (respectively negative) data correlation decreases (respectively increases) the degree of conflict between the objective function values. The generation of correlated data follows a multivariate uniform distribution of dimension 2 [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF]. As reported in Fig. 1(a), the coefficient ρ allows to tune the correlation between objective function values with a high accuracy. The considered bUBQP problem instances as well as an instance generator are available at the following URL: http://mocolib.sf.net/.

Characteristics of Small-Size bUBQP Instances

Thereafter, we study the impact of the density d and of the objective correlation ρ on the number of Pareto optimal solutions, Pareto local optimal solutions, supported solutions, and on the connectedness property of bUBQP instances. More particularly, we consider a density d ∈ {0.2, 0.4, 0.6, 0.8, 1.0} and an objective correlation ρ ∈ {-0.9, -0.7, -0.4, -0.2, 0.0, +0.2, +0.4, +0.7, +0.9}. The problem size is set to n = 18 in order to enumerate the solution space exhaustively. For each parameter setting, 30 independently generated random instances are considered. Experimental results are given in Fig. 1(b-f ). In the following, we provide a detailed analysis of these statistics.

Pareto Optimal Solutions

Fig. 1(b) shows the proportion of Pareto optimal solutions in the solution space. Interestingly, the density d does not affect the size of the Pareto set. However, the objective correlation ρ modifies the number of Pareto optimal solutions to several orders of magnitude. Indeed, almost 0.05% of the solution space correspond to non-dominated solutions for conflicting objectives (ρ = -0.9), whereas this number drops to less than 0.003% for correlated objective (ρ = +0.9). As a consequence, the larger the objective correlation ρ, the lower the cardinality of the Pareto set. This means that an algorithm is expected to take more time to identify the whole Pareto set when the objectives are in conflict.

Supported Solutions

In multi-objective optimization, scalarizing approaches consist in transforming the original problem into a single-objective one by means of an aggregation of the objective function values. A typical example is the weighted-sum scalarizing function [START_REF] Ehrgott | Multicriteria optimization[END_REF] that can be defined as follows.

g λ (x) = λ 1 • f 1 (x) + λ 2 • f 2 (x) (2)
where x ∈ X is a candidate solution, and λ = (λ 1 , λ 2 ), such that λ 1 , λ 2 ≥ 0, is a weighting coefficient vector. Supported solutions are non-dominated solutions which are optimal with respect to a weighted-sum aggregation of the objective functions. Their corresponding objective vectors are located on the boundary of the convex hull of the Pareto front [START_REF] Ehrgott | Multicriteria optimization[END_REF]. On the contrary, non-supported solutions are not optimal for any setting of the weighting coefficient vector λ. In order to explain the ability of scalarizing multi-objective optimization approaches to identify a large portion of Pareto optimal solutions, we should put the problemrelated properties in relation with the proportion of supported solutions. Fig. 1(c) shows the proportion of supported solutions in the Pareto set. Once again, the matrix density d has a very small influence. However, when the objective correlation increases, and despite the absolute number of supported solution actually gets lower, the proportion of supported solutions on the Pareto set increases. For highly correlated objectives (ρ = +0.9), nearly all Pareto optimal solutions are supported (this is even the case for some of the instances). On the contrary, for conflicting objectives (ρ = -0.9), only 15% of Pareto optimal solutions are supported. By putting this property in relation with algorithm design, we can assume that scalarizing approaches should be more suited to approximate the Pareto set of bUBQP instances with correlated objectives.

Connectedness

In the following, we describe some properties related to the connectedness of the Pareto set [START_REF] Ehrgott | Connectedness of efficient solutions in multiple criteria combinatorial optimization[END_REF]. We follow the definition of k-Pareto graph from [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF]. The k-Pareto graph is a graph P G k = (V, E), where each vertex in V corresponds to a Pareto optimal solution, and there is an edge e ij ∈ E between two nodes i and j only if the shortest distance between solutions x i and x j ∈ X, with respect to a given neighborhood, is below a bound k. For bUBQP, we adopt the Hamming distance on binary strings. This corresponds to the number of moves performed with the bit-flip neighborhood operator. Fig. 1(d) shows the ratio between the size of the largest connected component in the 1-Pareto graph (P G k=1 ) and the size of the Pareto set. The objective correlation ρ has a clear impact on this feature. Indeed, the proportion of Pareto optimal solutions in the largest connected component decreases from ρ = -0.9 to ρ = 0.4, and then slightly increases from ρ = 0.4 to ρ = 0.9. Overall, we can expect to reach 50% to 95% of the whole Pareto set by iteratively exploring the neighborhood of an approximation set initialized with at least one non-dominated solution. However, when there are several connected components in the 1-Pareto graph, it may happen that the distance between those components is small. Fig. 1(e) reports the smallest distance k such that the k-Pareto graph becomes connected, i.e. for all pairs of vertices x i , x j ∈ V in P G k , there is a path between x i and x j . When this minimal distance k is around 9, which is the average distance between random solutions for n = 18, we can conclude that the distance between Pareto optimal solutions is large. Actually, for bUBQP instances, this minimal distance is clearly smaller (between 2 and 3 in average). This means that finding a subset of non-dominated solutions can actually help to identify additional ones, which then may constitute a valuable asset for initializing local search algorithms.

Pareto Local Optimal Solutions

In Fig. 1(f ), we report the proportion of Pareto Local Optimal (PLO) solutions [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF] in the solutions space. A solution x ∈ X is a PLO with respect to a neighborhood structure N if there does not exist any neighboring solution x ′ ∈ N (x) such that x ≺ x ′ . As above, the neighborhood structure is taken as the 1-bit-flip, which is directly related to a Hamming distance 1. Once again, the distribution d does not seem to affect the number of PLO. However, similar to the trend observed on the Pareto set cardinality, the objective correlation ρ modifies the number of PLO to several orders of magnitude, from 20% of the solution space for ρ = -0.9 to less than 0.02% for ρ = +0.9. Therefore, by assuming that the difficulty for Pareto-based local search gets higher when the number of PLO is large, the difficulty of bUBQP instances might increase with the degree of conflict between the objectives.

Local Search for bUBQP

In this section, we give the working principles of four stochastic local search algorithms for identifying a Pareto set approximation to bUBQP instances. We start by introducing the algorithmic components shared by the different approaches. Then, we present the search strategies of two scalarizing local search algorithms, one Pareto-based local search algorithm, as well as a hybrid approach where a scalarizing and a Pareto local search phases are sequentially applied.

Main Ingredients

Neighborhood Relation. Similarly to the previous analysis, the neighborhood structure of the proposed local search algorithms is based on the 1-bit-flip operator: Two feasible solutions are neighbors if they differ exactly on one variable. In other words, a given neighbor can be reached by changing the value of a binary variable to its complement from the current solution. The size of the 1-bit-flip neighborhood structure is equal to the problem size n. As in the single-objective UBQP, each bUBQP objective function can be evaluated incrementally. We follow the fast incremental evaluation procedure proposed in [START_REF] Glover | Efficient evaluations for solving large 0-1 unconstrained quadratic optimisation problems[END_REF] to calculate the move gain of a given neighbor. For each objective, the whole set of neighbors is evaluated in linear time. As a consequence, the objective values of all neighboring solutions are evaluated in O(n) in the two-objective case.

Tabu Search. The tabu search algorithm proposed in [START_REF] Glover | Diversification-driven tabu search for unconstrained binary quadratic problems. 4OR-Q[END_REF] is reported to be one of the best-performing approaches for single-objective UBQP. In order to extend it to the multi-objective case, we consider a simple weighted-sum aggregation, as presented in Section 3.2, so that the initial objective vector values are (temporarily) transformed into a single scalar fitness value. Once the objective values of a given neighboring solution have been (incrementally) evaluated, we compute its scalar fitness value with respect to the weighted-sum problem (Eq. 2) for a given definition of the weighing coefficient vector. As a short-term memory, we maintain the tabu list as follows: Revisiting solutions is avoided within a certain number of iterations, called the tabu tenure. The tabu tenure of a given variable x i is denoted by tenure(i). Hence, variable x i will not be flipped again for a number of tenure(i) iterations. Following [START_REF] Glover | Efficient evaluations for solving large 0-1 unconstrained quadratic optimisation problems[END_REF], we set the tabu tenure of a given variable x i after it has been flipped as tenure(i) = tt + rand [START_REF] Glover | Diversification-driven tabu search for unconstrained binary quadratic problems. 4OR-Q[END_REF], where tt is a user-given parameter and rand [START_REF] Glover | Diversification-driven tabu search for unconstrained binary quadratic problems. 4OR-Q[END_REF] gives a random integer value in [START_REF] Aneja | Bicriteria transportation problem[END_REF][START_REF] Glover | Diversification-driven tabu search for unconstrained binary quadratic problems. 4OR-Q[END_REF]. From the set of neighbors produced by all non-tabu moves, we select the one with the best (highest) fitness value. However, all the neighbors are always evaluated, and a tabu move can still be selected if it produces a better solution than the current global best; this is called an aspiration criterion in tabu search [START_REF] Glover | Diversification-driven tabu search for unconstrained binary quadratic problems. 4OR-Q[END_REF]. The stopping condition is satisfied when no improvement has been performed within a given number of moves α, called the improvement cutoff. For more details on the tabu search algorithm for single-objective UBQP, the reader is referred to [START_REF] Glover | Diversification-driven tabu search for unconstrained binary quadratic problems. 4OR-Q[END_REF].

Scalarizing Local Search with Uniform Weights (SLS unif )

The first approach consists in solving different settings of the weighted-sum problem (Eq. 2) by means of multiple weighting coefficient vectors defined in a way that the whole region of the Pareto front is covered in the objective space.

For solving each scalarizing sub-problem, any algorithm for the resulting singleobjective problem version can potentially be applied. In our case, we use the tabu search algorithm detailed above as a (single-objective) solver. Let us consider a set of µ uniformly defined weighting coefficient vectors (λ 0 , . . . , λ i , . . . , λ µ-1 ), such that λ i 1 = i/(µ -1) and λ i 2 = 1λ i 1 . Each weighting coefficient vector λ i corresponds to a scalarizing sub-problem WS i . We start by identifying a highquality solution with respect to the first objective function, corresponding to the scalarizing sub-problem WS 0 , associated with the weighting coefficient vector λ 0 = (0, 1). The final solution is then used as a seeding solution for solving the next sub-problem WS 1 . We iterate this principle, each time the initial solution for sub-problem WS i being the one that is returned by the tabu search algorithm for the previous sub-problem WS i-1 . At last, in order to avoid a bias in the search process towards one objective, we re-run the same strategy by considering the reversed sequence of weighting coefficient vectors (λ µ-1 , λ µ-2 , . . . , λ 0 ). The algorithm outputs the union of non-dominated solutions generated during these two phases. The resulting scalarizing local search with uniform weights (SLS unif ) adapts the "double two-phase local search" from [START_REF] Paquete | A study of stochastic local search algorithms for the biobjective QAP with correlated flow matrices[END_REF] to bUBQP by using the single-objective tabu search procedure for solving scalarizing sub-problems.

Dichotomic Scalarizing Local Search (SLS dicho )

Similarly, the second approach is based on solving a sequence of scalarizing sub-problems, by means of a weighted-sum aggregation function, with the singleobjective tabu search algorithm. However, unlike SLS unif that defines them a priori, the weighting coefficient vectors are now iteratively determined based on the solutions identified at previous steps. The resulting SLS dicho approach follows the principles of dichotomic search from exact bi-objective optimization [START_REF] Aneja | Bicriteria transportation problem[END_REF], and adapt them to a local search engine strategy, similarly to [START_REF] Dubois-Lacoste | Improving the anytime behavior of two-phase local search[END_REF]. Notice that, by using any exact algorithm instead of the tabu search procedure, such a dichotomic search would output the (exact) set of supported solutions [START_REF] Ehrgott | Multicriteria optimization[END_REF]. Unfortunately, this would require to solve an NP-hard problem for each scalarizing sub-problem. Indeed, each sub-problem corresponds to a single-objective UBQP instance.

We start by identifying a high-quality solution for each separate objective. Let x 1 (resp. x 2 ) be the approximate solution found by tabu search for objective f 1 (resp. f 2 ). Both solutions are then added to a sequence U F = x 1 , x 2 , . . . , arranged in the decreasing order of f 1 -values. Next, at each step of the algorithm, we define a weighting coefficient vector

λ = f 2 (x 2 ) -f 2 (x 1 ), f 1 (x 1 ) -f 1 (x 2
) , corresponding to the sub-problem to be solved in the current iteration. It gives a search direction that is perpendicular to the segment defined by f (x 1 ) and f (x 2 ) in the objective space. Let x be the solution identified by tabu search for this definition of λ.

If f 1 (x 1 ) > f 1 (x) > f 1 (x 2 ) and f 2 (x 2 ) > f 2 (x) > f 2 (x 1
), then x is added to the sequence U F . Otherwise, we remove x 1 from U F and add it to an external set U T . Following [START_REF] Dubois-Lacoste | Improving the anytime behavior of two-phase local search[END_REF], for each scalarizing sub-problem, we use the solutions found in previous iterations to seed the search process of the current iteration. Based on preliminary experiments, both x 1 and x 2 are here used as an initial solution for two independent runs of the tabu search procedure based on λ. The SLS dicho algorithm iterates this principle until U F contains less than two elements, and returns the non-dominated solutions from U F ∪ U T .

Pareto Local Search (PLS)

Let us now consider a Pareto approach based on a set of solutions and local search principles. In contrast to scalarizing approaches, the selection process is here directly based on the Pareto dominance relation. A typical example is the Pareto Local Search (PLS) algorithm [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF]. An archive of mutually non-dominated solutions found so far is maintained in two different sets: V F for non-dominated solutions whose neighborhood has not yet been explored, and V T for solutions whose neighborhood has already been explored. These two sets are used in order to avoid a useless re-evaluation of a solution's neighborhood. The algorithm starts with a set of mutually non-dominated solutions to initialize V F , typically a single random solution. At each iteration, one unvisited solution is chosen at random from V F . All its neighboring solutions are (incrementally) evaluated and checked for insertion in the archive. The current solution is then discarded from V F and added to V T , and dominated solutions are removed from V F ∪ V T . The algorithm stops once V F is empty, i.e. all solutions from the archive are visited. PLS always terminates and returns a maximal Pareto local optimum set [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF].

Two-Phase Local Search (TP-LS)

The final algorithm consists in a hybrid two-phase approach, where SLS dicho and PLS are applied in a sequential way. It combines two fundamentally different and complementary search strategies: a scalarizing and a Pareto-based approach. In the first phase, SLS dicho is applied to identify a set of approximate supported solutions, as described in Section 4.3. This set of mutually non-dominated solutions is then used to initialize the archive V F of the PLS algorithm, and is further improved by exploring the neighborhood of its own content until no improvement is possible. Hence, contrary to the conventional PLS, the search process does not start with a single random solution, but with a set of good-quality solutions identified by a scalarizing approach. The performance of the designed two-phase local search algorithm (TP-LS) should be impacted by the connectedness property for the problem under consideration; the more connected the Pareto optimal solutions, the easier to identify new non-dominated solutions from identified ones. Notice that TP-LS shares similar principles with existing approaches proposed for other problem classes [START_REF] Dubois-Lacoste | A hybrid TP+PLS algorithm for bi-objective flow-shop scheduling problems[END_REF][START_REF] Lust | Two-phase Pareto local search for the biobjective traveling salesman problem[END_REF][START_REF] Paquete | A study of stochastic local search algorithms for the biobjective QAP with correlated flow matrices[END_REF]].

Experimental Analysis

Experimental Design

We conduct an experimental study on the influence of the problem size (n) and of the objective correlation (ρ) over the performance of the proposed local search algorithms for approximating the Pareto set of bUBQP problem instances. In addition, we consider the best-known approximation sets (best-known) identified by multiple variants of evolutionary and memetic algorithms proposed in [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF].

We investigate the following instance parameter setting: a problem dimension n ∈ {1000, 2000, 3000, 4000, 5000} and a correlation between the objective function values ρ ∈ {-0.5, -0.2, 0.0, +0.2, +0.5}. The density of the matrices is set to d = 0.8. One instance, generated at random, is considered per parameter combination. This leads to a total of 25 problem instances. A set of 30 runs per instance is performed for each algorithm. All the algorithms start with a random solution. The tabu tenure tt is set to n/150 and the improvement cutoff α is set to n. At last, for each phase of SLS unif , µ = 101 weighting coefficient vectors (λ 0 , . . . , λ i , . . . , λ 100 ) are uniformly defined as λ i 1 = i/100 and λ

i 2 = 1 -λ i 1 .
Since all the algorithms have a natural stopping condition, we measure their performance in terms of approximation set quality and computational cost. For each instance, we examine the quality of the Pareto set approximations identified by the competing algorithms in terms of hypervolume and epsilon indicators [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]. First, we compute the hypervolume relative deviation (hypervolume) as (hv(R)hv(A))/hv(R), where A ⊆ Z is an approximation set and R is a reference set. The reference set is the best-found approximation over all tested configurations for the instance under consideration. Let z - k (resp. z + k ) be the worst (resp. best) value obtained over all approximation sets for objective f k , the reference point z = (z 1 , z2 ) for the hypervolume calculation is set to zk

= z - k -(z + k -z - k ) • 10 -2 , k ∈ {1, 2}.
Additionally, the epsilon indicator (epsilon) gives the minimum multiplicative factor by which an approximation set has to be shifted in the objective space in order to weakly dominate the reference set. In both cases, a lower indicator-value is better.

Experimental Results

A summary of our computational results is presented in Table 1, following the presentation from [START_REF] Bader | HypE: An algorithm for fast hypervolume-based manyobjective optimization[END_REF]. The first line corresponds to the bUBQP instance with ρ = -0.5 and n = 1000, and reports the quality of the Pareto set approximation obtained by the different algorithms with respect to hypervolume. The average hypervolume relative deviation obtained by SLS unif , SLS dicho , PLS, TP-LS and best-known over the 30 executions is respectively 0.009, 0.006, 0.002, 0.000 and 0.031. The ranking obtained by means of a pairwise Wilcoxon signed-rank nonparametric statistical test gives the following order for this particular setting:

(1) TP-LS, (2) PLS, (3) SLS dicho , (4) SLS unif , and ( 5) best-known. Complementarily, Fig. 2 shows the average indicator-values for a subset of instances (the error bars indicate the confidence interval within a significance level of 10 -2 ). The results from best-known are omitted for a better readability.

Clearly, all the local search algorithms investigated in the paper largely improve over the previous best-known approximation sets from [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF]. Indeed, for the 25 bUBQP instances under investigation, best-known obtains the lowest rank for 23 of them in terms of hypervolume and epsilon. A simple approach like SLS unif is able to obtain better best-known results in all the instances but one. Among the algorithms proposed in the paper, SLS unif is repeatedly dominated by the others with respect to both indicators. The only notable exceptions are for instances with correlated objectives where SLS unif performs better than PLS in terms of hypervolume (while it is slightly worse in terms of epsilon), and for large-size instances with conflicting objectives where PLS encounters more difficulties compared with other approaches. In both cases, the reason seems to be that the approximation sets identified by PLS badly covers the lexicographically optimal regions of the Pareto front. This is also the reason why SLS dicho outperforms PLS on more than half of the instances with respect to hypervolume, whereas the same only happens four times with respect to epsilon.

Interestingly, the quality of the approximation sets identified by scalarizing approches (SLS unif , SLS dicho ) slightly increases with the objective correla-Table 1. Comparison of the competing local search algorithms and of the previous best-known approximation [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF] with respect to the hypervolume relative deviation (hypervolume) and to the unary multiplicative epsilon indicator (epsilon). The first value stands for the number of algorithms that statistically outperform the one under consideration with respect to a pairwise Wilcoxon signed-rank non-parametric statistical test with a p-value of 10 -2 by using a Bonferroni correction (lower is better). The number in brackets stands for the average indicator-value, rounded to 10 -3 (lower is better). Bold ranking values correspond to the best-performing algorithm for the instance and the indicator under consideration.

4 (0.037) 0.5 3000 2 (0.003) 1 (0.002) 3 (0.010) 0 (0.001) 3 (0.010) -0.5 4000 3 (0.007) 1 (0.003) 2 (0.006) 0 (0.000) 4 (0.092) -0.2 4000 3 (0.007) 1 (0.003) 2 (0.004) 0 (0.001) 4 (0.077) 0.0 4000 3 (0.007) 1 (0.003) 2 (0.005) 0 (0.001) 4 (0.093) 0.2 4000 2 (0.004) 1 (0.002) 3 (0.006) 0 (0.001) 4 (0.047) 0.5 4000 2 (0.003) 1 (0.002) 3 (0.008) 0 (0.001) 4 (0.014) -0.5 5000 2 (0.007) 1 (0.002) 3 (0.020) 0 (0.000) 4 (0.141) -0.2 5000 2 (0.007) 1 (0.003) 3 (0.008) 0 (0.001) 4 (0.130) 0.0 5000 3 (0.006) 1 (0.003) 2 (0.006) 0 (0.001) 4 (0.130

2

(1.008) +0.5 1000 3

(1.011)

3

(1.012)

2

(1.008) 0

(1.002)

1 (1.005) -0.5 2000 3 (1.009) 2 (1.007) 1 (1.003) 0 (1.001) 4 (1.026)
-0.2 2000 2 (1.008) 2 (1.008) 1 (1.003) 0 (1.001) 4 (1.027) 0.0 2000 3 (1.009) 2 (1.007) 1 (1.005) 0 (1.001) 4 (1.025) +0.2 2000 3 (1.009) 2 (1.008) 1 (1.004) 0 (1.001) 4 (1.019) +0.5 2000 2 (1.011) 1 (1.009) 1 (1.008) 0 (1.002) 4 (1.014)
-0.5

-0.2 tion ρ, as the proportion of supported solutions; see Fig. 1(c). By comparing SLS unif with SLS dicho , the later is always at least as good as the former for all the instances we investigated but one (ρ = +0.5, n = 1000). The reason is that the SLS unif algorithm is limited on the approximation set size that it is able to identify (a fixed number of 101 weighting coefficient vectors times two phases, i.e. 202 solutions at most); see Fig. 3. On the contrary, SLS dicho adaptively determines a number of weighting coefficient vectors based on the solutions it iteratively identifies. As a consequence, it takes advantage of manipulating an unbounded approximation set, that allows SLS dicho to obtain better indicatorvalues overall. Still, as shown in Fig. 3, the number of solutions identified by both scalarizing approaches, which only seek for supported solutions, are lower than Pareto-based approaches by several orders of magnitude. However, the number of solutions found by all approaches reduces with the objective correlation ρ, as the number of Pareto optimal solutions reported in Section 3.1.

Overall, in terms of approximation quality, there is a clear advantage to TP-LS. Actually, hybridizing SLS dicho and PLS allows to obtain statistically better approximation sets, in terms of hypervolume and epsilon, than all the other competing algorithms, for all the instances. In particular, there is a substantial improvement of the indicator-values, showing that TP-LS is consistently able to identify a high-quality approximation set, which is very close to the reference set, especially for instances with conflicting objectives. Indeed, more than 99.67% of the best-found hypervolume is covered by the approximation set identified by TP-LS in the worst case. The epsilon indicator-value is always less than 1.003. The weakness of PLS in identifying good-quality lexicographical solutions seems to be overcome by initializing the archive with high-quality scalarizing solutions. This means that finding some non-dominated solutions can actually help to identify additional ones, as conjectured in Section 3.3. More surprisingly, the hybrid TP-LS approach also allows to improve the performance of PLS in terms of computational time. As reported in Fig. 4, the running time of TP-LS is below the one of PLS for most of the instances. By analyzing this more carefully, TP-LS actually allows to generate much more candidate solutions, and as many potentially non-dominated solutions as PLS, by performing much less pairwise comparisons, particularly for large-size instances. As well, the overhead of TP-LS compared to only performing the first phase as in SLS dicho is almost insignificant. At last, Fig. 4 also reveals that the loss of SLS unif in terms of quality is in fact counter-balanced by a very short computing time, which is lower than all other algorithms by several orders of magnitude. Actually, one single run of SLS unif allows to improve the aggregated best-known results from [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF] by running faster than each single run of the algorithms from [START_REF] Liefooghe | A hybrid metaheuristic for multiobjective unconstrained binary quadratic programming[END_REF].

Conclusions

In this paper, we designed and analyzed stochastic local search algorithms for identifying a Pareto set approximation in bi-objective unconstrained binary quadratic programming. First, we designed a local search approach that iteratively identifies an approximate solution for multiple scalarizing sub-problems, such that the whole region of the Pareto front is covered in the objective space. The resulting SLS unif algorithm is based on a weighted-sum aggregation function, on a set of uniformly defined weighting coefficient vectors, and on a state-of-theart tabu search procedure for the single-objective version of the problem under consideration. The scalarizing sub-problems are solved in sequence, such that the solution found at a given iteration is used to initialize the search process of the subsequent iteration. SLS unif allows to obtain a substantial improvement over the best-know solutions from previous studies with much less computations. Furthermore, the computational ressources required for instances with thousands of variables is less than a few minutes. However, the given number of weighting coefficient vectors severely limits the cardinality of the approximation set identified by SLS unif . This user-defined parameter cannot be easily set a priori without performing an expensive experimental campaign. For this reason, we considered an improved scalarizing approach, based on the dichotomic search principles, that defines the weighting coefficient vector to be used in the current iteration based on solutions identified in previous iterations. As a consequence, there is no user-defined limit on the number of solutions manipulated by this search strategy. The designed SLS dicho algorithm allows to improve the results of SLS unif in terms of quality, while inducing an extra cost in terms of computing time. Still, for these two scalarizing approaches, that both seek supported optimal solutions only, the cardinality of the obtained approximation set is less than those of Pareto-based approaches. In addition, we highlighted that the performance of scalarizing local search decreases with the degree of conflict between the objectives, following the proportional number of supported solutions. Next, a Pareto local search algorithm, that iteratively explores the neighborhood of an archive of mutually non-dominated solutions, obtains similar results in terms of quality indicator-values, while requiring even more computational resources. Indeed, the PLS strategy obtains larger approximation sets while exploring much less candidate solutions. The bottleneck of its effectiveness is the number of comparisons required to maintain the archive, while the bottleneck of its efficiency is a poor quality in finding strong solutions at the extremes of the Pareto front. At last, we designed a two-phase algorithm that applies SLS dicho and PLS in a sequential way, the output of the former being the input of the latter. The TP-LS approach significantly surpasses the other algorithms over all the configurations we experimented. The computational overhead is negligible compared with the stand-alone SLS dicho approach, TP-LS being able to generate as many solutions as PLS while performing much less comparisons. This hybrid algorithm profits from the closeness that exists between optimal solutions. Starting with a subset of approximate (supported) non-dominated solutions, TP-LS is able to identify additional ones, then improving the overall approximation set quality.

Extending the stochastic local search approaches investigated in the paper to unconstrained binary quadratic programming with more than two objectives would improve even more the expressive ability of the problem formulation. This would provide a more general unifying modeling and solution framework for multi-objective optimization that could potentially enable an efficient reformulation and resolution of a wide class of large-scale and NP-hard multi-objective combinatorial optimization problems.

Fig. 1 .

 1 Fig. 1. (a) Spearman correlation coefficient between the objective function values, (b) ratio of the number of Pareto optimal solutions to the solution space size, (c) ratio of the number of supported solutions to the Pareto set size, (d) ratio of the size of the largest connected component of the Pareto graph for Hamming distance 1 to the Pareto set size, (e) minimal Hamming distance to connect the Pareto graph, and (f ) ratio of Pareto local optimal solutions to the solution space size, with respect to the objective correlation ρ, for d = 0.2 (•), d = 0.4 ( ), d = 0.6 (⋄), d = 0.8 (△) and d = 1.0 (▽). For each parameter setting, average values and confidence intervals (with a significance level of 10 -2 ) are reported over 30 independently generated random instances. The problem size is n = 18. Notice the log-scale on the y-axis for (b) and (f ).

Fig. 2 .

 2 Fig. 2. Comparison of SLSunif (•), SLSdicho ( ), PLS (⋄) and TP-LS (△) with respect to hypervolume (top) and epsilon (bottom) for ρ = -0.5 (left), ρ = 0.0 (center), and ρ = +0.5 (right). A lower value is better.

Fig. 3 .

 3 Fig. 3. Comparison of SLSunif (•), SLSdicho ( ), PLS (⋄) and TP-LS (△) with respect to the size of the approximation set found for ρ = -0.5 (left), ρ = 0.0 (center), and ρ = +0.5 (right). Notice the log-scale on the y-axis.

Fig. 4 .

 4 Fig. 4. Comparison of SLSunif (•), SLSdicho ( ), PLS (⋄) and TP-LS (△) with respect to to the CPU time (top), the number of comparisons (middle), and the number of visited solutions (bottom), for ρ = -0.5 (left), ρ = 0.0 (center), and ρ = +0.5 (right). Notice the log-scale on the y-axis.


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

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

Fabio Daolio
email: fdaolio@shinshu-u.ac.jp

Hernan Aguirre

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

Hernán Aguirre

A feature-based performance analysis in evolutionary multiobjective optimization

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Black-box multiobjective combinatorial optimization problems are characterized by a discrete solution space and by multiple objective functions, such as cost, profit, or risk, that are ill-defined, computationally expensive, or for which an analytical form is not available. Due to the black-box nature of the objective functions, problem-specific algorithms are usually excluded to identify or approximate the Pareto set, so that an increasing number of general-purpose evolutionary multiobjective optimization (EMO) algorithms and other randomized search heuristics have been proposed in recent years [START_REF] Coello Coello | Evolutionary Algorithms for Solving Multi-Objective Problems[END_REF]. However, the overall amount of understanding about the pros and cons of different EMO algorithm designs and configurations with respect to a given problem structure is rather scarce. Due to the increasing number and complexity of black-box multiobjective optimization problems and algorithms, one of the most difficult challenges is to devise and exhibit a number of general-purpose problem characteristics and statistical methodologies allowing to explain the dynamics and the performance of EMO algorithms.

Recently, a few attempts to explain the performance of randomized search heuristics based on relevant fitness landscape features have been proposed for single-objective optimization problems of continuous and combinatorial nature; see e.g. [START_REF] Daolio | Local optima networks and the performance of iterated local search[END_REF][START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF][START_REF] Mersmann | A novel feature-based approach to characterize algorithm performance for the traveling salesperson problem[END_REF]. In this paper, we address the issue of feature-based performance analysis for EMO algorithms according to the main characteristics of 0-1 multiobjective optimization problems. We first extend our previous works by summarizing a number of problem properties and fitness landscape features for blackbox 0-1 multiobjective optimization [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF]. They include features extracted from the problem input data, like variable correlation, objective correlation, and objective space dimension [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF], as well as features from the Pareto set [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF][START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF], the Pareto graph [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF] and the ruggedness and multimodality of the fitness landscape [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. Then, we analyze the correlation between those features and the performance of an EMO algorithm. More particularly, we investigate the expected running time of the global SEMO algorithm [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF] with restart to identify a (1+ε)-approximation of the Pareto set on a large number of small-size enumerable multiobjective NKlandscapes with objective correlation, i.e. ρMNK-landscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. Our analysis shows the relative influence of each individual problem feature on the algorithm performance. In particular, the running time of global SEMO appears to be predominantly impacted by the ruggedness of the fitness landscape, more than other features like the number of Pareto optimal solutions. Additionally, we investigate different formulations of a multiple linear regression model. This allows us to discuss the joint effect of different subsets of features in capturing the dynamics of the algorithm. The ruggedness and the multimodality, but also the number of objectives, the correlation between them, and the hypervolume of the Pareto front turn out to be the most impactful characteristics that allow to explain the performance of global SEMO for ρMNK-landscapes.

The remainder of the paper is organized as follows. Section 2 details the problem and algorithm settings of our analysis. Section 3 summarizes the problem features under consideration in the paper. Section 4 introduces different regression models to explain the performance of global SEMO for enumerable ρMNK-landscapes. Section 5 concludes the paper and suggests further research.

Problem and Algorithm Settings

In this paper, we are interested in the ability of evolutionary multiobjective optimization (EMO) algorithms to identify a Pareto set approximation for black-box multiobjective combinatorial optimization problems. In particular, we investigate the (estimated) running time of global SEMO [START_REF] Laumanns | Running time analysis of evolutionary algorithms on a simplified multiobjective knapsack problem[END_REF] with restart to identify a (1 + ε)-approximation of the Pareto set on a large bench of enumerable ρMNKlandscapes with different structural properties. We consider the maximization of an objective function vector f = (f 1 , . . . , f m ) over the discrete set of solutions X = {0, 1}

n , where m is the number of objectives, and n is the problem size. X is the solution space, and Z = f (X) ⊆ IR m is the objective space. A solution x ∈ X is dominated by a solution x ′ ∈ X if ∀i ∈ {1, . . . , m}, f i (x) f i (x ′ ) and ∃i ∈ {1, . . . , m} such that f i (x) < f i (x ′ ). The set of solutions that are not dominated by any other is the Pareto set, and its image in the objective space is the Pareto front.

ρ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]. They extend single-objective NK-landscapes [START_REF] Kauffman | The Origins of Order[END_REF] and multiobjective NK-landscapes 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 solution space is X = {0, 1}

n . The parameter k refers to the number of variables that influence a particular position from the bit-string. 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 function f i is to be maximized. The problem can be formalized as follows.

max f i (x) = 1 n n j=1 f ij (x j , x j1 , . . . , x j k ) , i ∈ {1, . . . , m} s.t. x j ∈ {0, 1} , j ∈ {1, . . . , n} (1)
As in the single-objective case, each separate objective function 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 . Indeed, for each objective f i , i ∈ {1, . . . , m}, and each variable x j , j ∈ {1, . . . , n}, 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 the range [0, 1]. As a consequence, the individual contribution of a variable x j depends on the value of x j , as well as on the values of k < n other variables {x j1 , . . . , x j k }. In this work, the epistatic interactions, i.e. the k variables that influence the contribution of x j , are set uniformly at random among the (n -1) variables other than x j , following the random neighborhood model from [START_REF] Kauffman | The Origins of Order[END_REF]. By increasing the number of epistatic interactions k from 0 to (n -1), problem instances can be gradually tuned from smooth to rugged. In ρMNK-landscapes, f ij -values additionally follow a multivariate uniform distribution of dimension m, defined by an m× m positive-definite symmetric covariance matrix (c pq ) such that c pp = 1 and c pq = ρ for all p, q ∈ {1, . . . , m} with p = q, were ρ > -1 m-1 defines the objective correlation degree; see [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF] for details. The positive (respectively negative) data correlation ρ allows to decrease (respectively increases) the degree of conflict between the objective function values. The same correlation coefficient ρ is then defined between all pairs of objectives, and the same epistatic degree k and epistatic interactions are set for all the objectives.

Global SEMO

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

n → Z such that Z ⊆ IR m , like in ρMNK-landscapes. It maintains an unbounded archive A of non-dominated solutions found so far. The archive is initialized with one random solution from the solution space. At each iteration, one solution is chosen at random from the archive x ∈ A. Each binary variable from x is independently flipped with a rate of 1 n in order to produce an offspring solution x ′ . This mutation operator is ergodic, meaning that there is a non-zero probability of jumping from any point to any other point in the solution space. The archive is then updated by keeping the non-dominated solutions from A ∪ {x ′ }. 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 budget of function evaluations.

Performance Measure

Let ε be a constant value such that ε 0. The (multiplicative) ε-dominance relation ( ε ) can be defined as follows. For all x, x ′ ∈ X,

x ε x ′ if f i (x) (1 + ε) • f i (x ′ ), ∀i ∈ {1, . . . , m}. A set X ε ⊆ X is an (1 + ε)-
approximation of the Pareto set if for any solution x ∈ X, there is one solution x ′ ∈ X ε such that x ε x ′ . This is equivalent to finding an approximation set whose multiplicative epsilon quality indicator value with respect to the (exact) Pareto set is lower than (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 here 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 running time 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 running time" [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 with restarts. Indeed, after (t -1) failures, each one requiring T f evaluations, and the final successful run with T s evaluations, the total running time is T = t-1 i=1 T f + T s . By taking the expectation value and by considering that the probability of success after (t -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 ] (2)
In our case, the success rate p s is estimated with the ratio of successful runs over the total number of executions (p s ), the expected running time for unsuccessful runs E[T f ] is set to a constant limit on the number of function evaluation calls T max , and the expected running time 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 t s ts i=1 T i (3)
where t 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 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 solution space exhaustively. The solution space size is then 2 18 . A set of 30 different landscapes, independently generated at random, are considered for each parameter combination: ρ, m, and k. They are made available at the following URL: http://mocobench.sf.net.

We set a target ε = 0.1. The time limit is set to T max = 2 n • 10 -1 < 26 215 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 of CPU resources available: m = 5 and ρ ∈ {-0.2, 0.0}. Those experiments have then been discarded due to missingness. Overall, this represents a total amount of 2 980 instances times 100 executions, that is 298 000 simulation runs.

Features to Characterize Problem Difficulty

In Table 1, we give a number of general-purpose problem features, either directly extracted from the problem instance (low-level features), or computed from the enumerated Pareto set and solution space (high-level features). Obviously, since the features require the solution space to be completely enumerated, they are not practical for performance prediction purposes. However, we still include them in order to examine their impact on the algorithm performance. For the case of ρMNK-landscapes, the neighborhood is induced by the bit-flip operator, which is directly related to the Hamming distance between solutions. For the computation of the hypervolume, the reference point is set to the origin z ⋆ = (0.0, . . . , 0.0). For a more comprehensive explanation of those features and a correlation analysis between them, we refer to [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF]. In the next section, we relate the value of those features for enumerable ρMNK-landscapes to the performance of G-SEMO.

Problem Features vs. Algorithm Performance

In this section, we conduct a linear regression analysis on the correlation between the problem features presented in the previous section and the performance of G-SEMO. The algorithm performance is defined as the expected running time ert, in terms of the number of evaluation function calls, required by the algorithm to identify a (1 + ε)-approximation of the Pareto set. We first detail our methodological setup. Then, we analyze the individual as well as the joint impact of problem features on the algorithm performance. At last, we compare the accuracy of regression models for different objective space dimensions.

Methodological Setup

Linear Regression. In order to provide an explanatory model for the algorithm performance, we perform a linear regression, whose general model can be formalized as follows:

y = β 0 + β 1 • v 1 + β 2 • v 2 + . . . + β p • v p + e ( 4
)
where y is the response variable, (v 1 , v 2 , . . . , v p ) are the explanatory variables, and e is the usual error term. In our case, the response variable to be explained is the expected running time of G-SEMO: log(ert), and the p explanatory variables correspond to selected problem features as detailed in Section 3. The response of the linear model is here log-transformed in order to better approach linearity; see Fig. 1 (top-left). Also, an order of magnitude is in general sufficiently relevant for the running time of an EMO algorithm. Such linear regression models are usually fitted using an ordinary least-square minimization. A least-square estimator βi is then produced for each regression coefficient β i , i ∈ {1, . . . , p}, by minimizing the sum of squared residuals between the fitted and the observed values. The case of a single explanatory variable (p = 1) is known as simple linear regression, whereas the extension to more than one explanatory variables (p > 1) known as multiple linear regression. Notice that, although a linear regression model might not be able to reliably catch the existing correlation between explanatory variables as well as non-linear dependencies with the response variable, it has the advantage of being simple and easily interpretable.

Regression Accuracy. In the following, when measuring the accuracy of a linear regression model, we will be interested in the following statistics:

-The absolute correlation coefficient (r) measures the linear association between the predicted and the actually observed values (the conventional Pearson correlation coefficient is here used). Its absolute value ranges from 0 to 1. The closer r to 1.0, the better the fitting. Actually, an r-value of 1.0 indicates that the linear regression line perfectly fits the data. -The mean absolute error (MAE) measures the average value of the absolute difference between the values predicted by the regression model and the values actually observed (the residuals). It aggregates the magnitudes of the prediction errors into a single predictive power measure to compare different models. Clearly, the lower the MAE, the better the regression model. -The root mean-square error (RMSE) measures the square root of the average squared difference between the values predicted by the regression model and the values actually observed. Similarly, the lower the RMSE, the better the regression model. Notice that the RMSE tends to favor a regression model that avoids large errors even though it produces a less satisfactory fit overall, whereas the MAE tends to favor a regression model that produces occasional large errors while being reasonably good on average. -The relative absolute error (%RAE) corresponds to the MAE relative to the basic model that always predicts the mean, i.e. when no explanatory variables are used in the regression model (p = 0). As a consequence, smaller values are better, and a RAE higher than 100% indicates that the corresponding model is actually worse than this basic model in terms of MAE. -The root relative squared error (%RRSE) corresponds to the RMSE relative to the basic model that always predicts the mean.

For each of those statistics, we report the values observed on the training set, i.e. the set of data used to build the model. In addition, we perform a 10-fold cross-validation in order to assess how the results of the regression model generalize to an independent data set. The original data set is randomly partitioned into 10 samples of equal size. The cross-validation process is repeated 10 times, with each of the samples being used exactly once as the validation data. For each sample, the above statistics are computed, and then averaged in order to produce a cross-validated r, MAE, RMSE, %RAE, and %RRSE value. Data Preparation. Although features might usually have to be normalized appropriately in order to get rid of scaling issues and for a fair comparison between them, we here chose not to normalize them in order to ease the interpretation of the different models. Notice that normalizing the features would only result in a change on the value of the estimated regression coefficient βi , i ∈ {1, . . . , p}. Moreover, this is not an issue within our experiments since we explicitly generate instances of the same size (n = 18) that takes their objective values in the same hyper-box [0, 1] m .

Individual Impact of Problem Features

The scatter plots between each feature and the log-transformed estimated running time of G-SEMO log(ert) is reported in Fig. 1. Notice that some features (npo, nconnec, nplo) have been log-transformed in order to better approach linearity. Additionally, the statistics of all possible simple linear regression models, one for each feature, are reported in Table 2, from the lowest to the highest absolute correlation value. The individual impact of each feature is analyzed below. First of all, four features are not directly linearly correlated to the expected running time of G-SEMO: the number of objective functions m, the objective correlation ρ, the cardinality of the Pareto set log(npo), and the hypervolume hv. For each of these features, the absolute correlation coefficient value is under 0.25, and the prediction error is around the one of the most basic model that always predicts the mean. Surprisingly, there is no direct connection with the two main low-level features from problem input data m and ρ. At least, the link between the running time and those features is not a direct linear correlation, but a more complex model will be analyzed in the next section. As well, the cardinality and the hypervolume of the Pareto set, features closely related to the final goal of the search process, do not explain the variance of log(ert) by themselves.

The features related to the connectedness of the Pareto set are all weakly correlated to log(ert). The absolute correlation coefficients of the number of connected components log(nconnec), the proportional size of the largest connected component lconnec, and the minimum distance to be connected kconnec are between 0.37 and 0.49. The more connected the Pareto set, the smaller the running time of G-SEMO. The algorithm performance is also moderately correlated with the average and maximal distance between Pareto optimal solutions, avgd and maxd (the absolute correlation coefficient values are 0.60 and 0.47, respectively). The larger the distance between Pareto optimal solutions in the solution space, the larger the running time of G-SEMO. Interestingly, the cardinality of the Pareto set has a smaller impact on the performance of G-SEMO than the distance between solutions in the Pareto set. Moreover, the multimodality of the landscape, in terms of the number of Pareto local optimal solutions log(nplo), is moderately correlated to the running time of G-SEMO: the more Pareto local optima, the longer the running time (the correlation coefficient is 0.46).

At last, the only strong correlation appears with the feature related to the ruggedness of the landscape. Indeed, the number of epistatic interactions k is highly correlated to the efficiency of G-SEMO (the correlation coefficient is 0.85). The more rugged the landscape, the longer it takes to identify a (1 + ε)-approximation of the Pareto set. In other words, by taking the features individually, the model based on k is the one that gives the highest accuracy. On average, it allows to predict the logarithm of the runtime of G-SEMO within ±0.77 of the observed value. Since the RMSE is much larger (1.00), this sug-gests that the deviation to this average value might be large, as we can see also on Fig. 1 (first line, second column). Such a regression accuracy is around twice better than the most basic model that always predicts the mean (±1.53), and largely better than the second more accurate simple linear regression model based on avgd (±1.15).

Having the ruggedness of the landscape k as a more important (individual) feature than the number of objectives m and the objective correlation ρ might be surprising at first sight. Indeed, the number of Pareto optimal solutions is known to increase exponentially with the number of objectives and the degree of conflict between them; see e.g. [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]. However, let us remind that the algorithm under consideration in the paper (G-SEMO) actually handles an unlimited approximation set size. It is then only slightly affected by the minimum number of solutions required to obtain a (1 + ε)-approximation of the Pareto set [START_REF] Papadimitriou | On the approximability of trade-offs and optimal access of web sources[END_REF]. Actually, depending on ρ, the estimated expected running time of G-SEMO is 23 to 118 times larger for rugged two-objective instances than smoother five-objective instances. For instance, when ρ = 0.2, the average ert-value is equal to 61 836 for k = 10 and m = 2 while it is equal to 2 691 for k = 2 and m = 5. Similarly, for ρ = 0.9, the average ert-value is 33 922 for k = 10 and m = 2, and only 287 for k = 2 and m = 5. In accordance with known results from single-objective optimization [START_REF] Barnett | Ruggedness and neutrality -the NKp family of fitness landscapes[END_REF], the ruggedness of the landscape k seems to largely impact the running time of EMO algorithms.

Overall, analyzing the individual impact of problem features supports the hypothesis that the structural properties identified in the previous section can help to understand the performance of a simple dominance-based EMO algorithm like G-SEMO for ρMNK-landscapes. In the next section, different multiple linear regression models are examined in order to better explain the running time of G-SEMO, based on the joint effect of these problem-related characteristics.

Joint Impact of Problem Features

We start by fitting the response variable log(ert) against all the low-level and high-level features presented in Section 3. The statistics related to the model accuracy are provided in Table 3 (line 1). For this complete multiple linear regression model, the correlation coefficient is over 0.9. Overall, it allows to explain the performance of G-SEMO with a much higher accuracy compared to the simple linear regression model based on k only, and outperforms the basic model that always predicts the mean by around 60%. Moreover, the model has a high degree of generalization. Indeed, the RMSE and the cross-validated RMSE are very close to each other. The same happens with the MAE (the difference between both is always under 10 -2 ). The scatter plot of the actual vs. the predicted performance values is given in Fig. 2 (left). This allows us to visualize how the model accuracy varies depending on the hardness of the problem instance: The model seems to slightly underestimate the runtime for easier and harder instances whereas it rather overestimates it for an intermediate instance difficulty. Although the regression coefficients related to each feature are not interpretable due to the different scaling of the metric values, not all regression coefficients are statistically significant in a general multiple linear regression model. In order to eliminate the influence of the less significant regression coefficients, we proceed by backward elimination. Starting from the inclusion of all features, we iteratively remove the feature that has the lowest impact on the increase of the MAE until no feature remains. This allows us to produce a ranked list of features by traversing the feature space from one side to the other and recording the order that attributes are deleted. Hence, the attributes that are deleted in the last steps have a more meaningful impact on the model. The steps of the backward elimination are sketched in Table 3, where one feature is removed at every line. Notice that a forward selection, that does the opposite procedure of iteratively adding attributes, ends up with a similar ranking on the importance of features, except that maxd is in the latter case more important than log(npo) and log(nconnec) (detailed results are not reported due to space limitation).

This feature selection analysis allows us to gain further insights about which subset of features obtains the highest accuracy. Indeed, the error increase is almost insignificant until the deletion of ρ in the model (line 8 in Table 3), where the correlation coefficient drops from 0.91 to 0.89, the MAE rises from 0.58 to 0.61 and the RMSE rises from 0.77 to 0.80. We can then conclude that ρ and subsequent attributes constitute the most significant subset for explaining the algorithm performance. Actually, a more compact model, with only the five most significant features, constitute an acceptable alternative, and has almost the same accuracy than the full model; see also Fig. 2 (right). Once again, the most important feature seems to be k, now followed by m, log(nplo), hv and ρ, in the order of importance. As a consequence, although they are not able to catch all the variations of log(ert) individually (see Section 4.2), the joint effect of all three low-level features from the problem input data (k, m and ρ) is relevant for explaining the running time of G-SEMO. Moreover, there is one high-level feature related to the hypervolume of the Pareto set hv and to the multimodality log(nplo). As well, the number of Pareto optimal solutions is not a significant addition to the regression model. We attribute this to the fact that the hypervolume incorporates a more relevant information related to the Pareto front for the algorithm behavior.

Explanatory Models vs. Objective Space Dimension

In this section, we build a separate regression model for each objective space dimension m ∈ {2, 3, 5}. The statistics related to the model accuracy are provided in Table 4. The information for the complete model mixing all m-values is also given in the tables in order to facilitate the comparison between the different models. Additionally, the scatter plot of the actual vs. the predicted performance values is given in Fig. 3. First of all, the regression model for a particular number of objectives is always slightly more accurate than the global model for all m-values (whatever the m-value, there is an improvement in terms of correlation, %RAE or %RRSE over the complete model). Actually, the worsening in terms of MAE or RMSE for larger m-values is only an artefact that the expected running time of G-SEMO increases with the number of objectives; see Fig. 1 (first line, third column). This means that constructing a regression model per objective space dimension only allows to reduce the prediction error to a very small extent.

As reported in Table 5, by focusing on the five most relevant features identified in the previous section (k, m, log(nplo), hv and ρ, following the order of importance), we are also able to construct a more compact regression model for each m-value with a satisfactory response compared to the model using all features. The difference in the correlation coefficient and the average error between the models with all features and their compact counterpart is always less than 2 • 10 -2 . In fact, applying a feature selection procedure by backward elimination for each of the models until the MAE increases by more than 10 -2 ends up with following subset feature selection, following the order of importance (detailed results are omitted due to space restriction):

k, hv and log(nplo) for m = 2; -k, log(nplo) and hv for m = 3;

k and ρ for m = 5.

Interestingly, this means that the objective correlation ρ is not a relevant feature for the models with m ∈ {2, 3}. We attribute this to the fact that the number of Pareto local optimal solutions log(nplo) increases with ρ in these cases, whereas the correlation between both is much lower for m = 5. For the same reason, ρ is selected over hv in the latter case. In addition, given that the proportion of Pareto local optimal solutions in the solution space increases exponentially with m (we know, for instance, that more than 95% of the solution space correspond to Pareto local optimal solutions for ρ = -0.2 and m = 5 in average [START_REF] Verel | On the structure of multiobjective combinatorial search space: MNK-landscapes with correlated objectives[END_REF]), the multimodality of the landscape, corresponding to log(nplo), is not relevant anymore for larger objective space dimensions. However, in all cases, the ruggedness of the landscape k is again the most relevant feature.

Conclusions

In this paper, we investigated the impact of problem features on the running time of a simple dominance-based EMO algorithm with restart, that maintains an unbounded archive of non-dominated solutions found so far. The topology of an arbitrary problem instance, in terms of ruggedness, multimodality, objective space dimension, objective correlation, cardinality and hypervolume of the Pareto set, as well as distance and connectedness between non-dominated solutions, has been examined for a large set of enumerable multiobjective NKlandscapes with objective correlation. First, a simple linear regression analysis revealed that the ruggedness of the landscape had the more critical effect on the algorithm performance. Second, a more-advanced multiple linear regression analysis allowed us to highlight the more significant subset of problem features.

As in the single-objective case [START_REF] Barnett | Ruggedness and neutrality -the NKp family of fitness landscapes[END_REF][START_REF] Kauffman | The Origins of Order[END_REF], the ruggedness and the multimodality of the landscape affect the algorithm running time to a large extent. Additionally, the number of objectives, the correlation between them, and the hypervolume of the Pareto front to be covered are all jointly impactful in the multiobjective case. At last, although problem features have a different impact depending on the objective space dimension, the degree of explanation they are able to provide together is always as meaningful for the algorithm performance. Overall, our feature-based analysis was able to highlight the main relationships between the structural properties of the landscape and the performance of the algorithm. This allowed us to better understand the behavior and the performance of this EMO algorithm class. The problem characteristics under analysis in the paper validate the relevance of our methodology for explaining the performance of EMO approaches. However, it remains an open question if there exist supplementary features that could better capture the problem difficulty for different problem and algorithm classes, and if more general regression models would allow to better apprehend the correlations among the features as well as their (non-linear) dependencies with the algorithm performance. Furthermore, the goal of the paper was on understanding the algorithm behavior and performance rather than blindly recommending the best-performing approach, but a natural extension for future research is to investigate the prediction power of the regression models proposed in the paper, based on existing works from single-objective optimization [START_REF] Hutter | Algorithm runtime prediction: Methods & evaluation[END_REF]. Following the algorithm selection problem formulated by Rice in the 1970s [START_REF] Rice | The algorithm selection problem[END_REF], this would allow us to design a portfolio approach for selecting the most appropriate algorithm configuration, based on a relevant structural characterization of the multiobjective problem instance to be solved. For that purpose, extending our paradigm with more-advanced regression models based on problem features that can be estimated inexpensively for large-size instances is currently under investigation.

  low-level features (k) Number of variable (epistatic) interactions (m) Number of objective functions (ρ) Correlation between the objective function values high-level features (npo) Number of Pareto optimal solutions [1, 8] (hv) Hypervolume value [18] of a the Pareto set [1] (avgd) Average distance between Pareto optimal solutions [11] (maxd) Maximum distance between Pareto optimal solutions [11] (nconnec) Number of connected components in the Pareto set [14] (lconnec) Proportion of the largest connected component of the Pareto set [10] (kconnec) Minimal Hamming distance to connect the Pareto set [14] (nplo) Number of Pareto local optimal solutions [15]

Fig. 1 .

 1 Fig. 1. Histogram of the distribution of log(ert)-values over all the instances (top-left), as well as scatter plots and regression lines for each feature vs. log(ert).

Fig. 2 .

 2 Fig. 2. log(ert) vs. fitted values for (left) the model with all features; and (right) the model with a selected subset of features (i.e. k, m, log(nplo), hv, ρ).

Fig. 3 .

 3 Fig. 3. log(ert) vs. fitted values for (left) m = 2, (middle) m = 3, (right) m = 5. All features are part of the models.

Table 1 .

 1 Summary of low-level and high-level features investigated in the paper.

Table 2 .

 2 Summary statistics of simple linear regression models, each one being based on a single problem feature. Values are rounded to 10 -2 .

				training set			10-fold cross validation
	feature	r MAE RMSE %RAE %RRSE	r MAE RMSE %RAE %RRSE
	ρ	0.03	1.52	1.89	99.66	99.94 0.01	1.53	1.88	99.75	100.00
	hv	0.04	1.52	1.89	99.66	99.94 0.02	1.53	1.89	99.73	99.99
	log(npo)	0.08	1.54	1.88 100.61	99.70 0.07	1.54	1.88 100.66	99.74
	m	0.21	1.51	1.85	98.63	97.86 0.20	1.51	1.85	98.69	97.91
	kconnec	0.37	1.38	1.76	89.99	93.03 0.37	1.38	1.76	90.01	93.06
	log(nconnec) 0.40	1.44	1.73	94.21	91.47 0.40	1.44	1.73	94.27	91.52
	log(nplo)	0.46	1.41	1.67	92.53	88.56 0.46	1.42	1.67	92.57	88.59
	maxd	0.47	1.30	1.67	85.07	85.07 0.47	1.30	1.67	85.09	88.36
	lconnec	0.49	1.31	1.65	85.67	87.35 0.49	1.31	1.65	85.69	87.37
	avgd	0.60	1.15	1.51	75.31	79.95 0.60	1.15	1.51	75.33	80.00
	k	0.85	0.77	1.00	50.67	52.95 0.85	0.78	1.00	50.67	52.96
	none	0.00	1.53	1.89 100.00	100.00 0.04	1.53	1.89 100.00	100.00

Table 3 .

 3 Summary statistics of the multiple linear regression model with backwardelimination feature selection. Values are rounded to 10 -2 .

				training set			10-fold cross validation
		r MAE RMSE %RAE %RRSE	r MAE RMSE %RAE %RRSE
	all features	0.91	0.58	0.76	37.74	40.34 0.91	0.58	0.76	37.86	40.48
	\ maxd	0.91	0.58	0.76	37.73	40.38 0.91	0.58	0.76	37.84	40.50
	\ log(nconnec) 0.91	0.58	0.76	37.73	40.49 0.91	0.58	0.77	37.83	40.60
	\ log(npo)	0.91	0.58	0.77	37.72	40.53 0.91	0.58	0.77	37.80	40.62
	\ lconnec	0.91	0.58	0.77	37.74	40.54 0.91	0.58	0.77	37.81	40.61
	\ avgd	0.91	0.58	0.77	37.83	40.58 0.91	0.58	0.77	37.90	40.64
	\ kconnec	0.91	0.58	0.77	37.95	40.66 0.91	0.58	0.77	38.00	40.71
	\ ρ	0.91	0.61	0.80	39.98	42.20 0.91	0.61	0.80	40.04	42.26
	\ hv	0.89	0.67	0.87	43.59	46.29 0.89	0.67	0.87	43.63	46.33
	\ log(nplo)	0.88	0.70	0.91	45.83	48.19 0.88	0.70	0.91	45.86	48.21
	\ m	0.85	0.77	1.00	50.67	52.95 0.85	0.78	1.00	50.67	52.96
	\ k	0.00	1.53	1.89 100.00	100.00 0.04	1.53	1.89 100.00	100.00

Table 4 .

 4 Summary statistics of the multiple linear regression models for all objective space dimensions (∀m) and for each individual objective space dimension m ∈ {2, 3, 5}. All features are part of the models. Values are rounded to 10 -2 .

				training set			10-fold cross validation
		r MAE RMSE %RAE %RRSE	r MAE RMSE %RAE %RRSE #inst
	∀m	0.91	0.58	0.76	37.74	40.34 0.91	0.58	0.76	37.86	40.48	2980
	m = 2 0.92	0.52	0.68	35.86	38.62 0.92	0.53	0.68	36.15	38.95	1350
	m = 3 0.92	0.52	0.69	36.12	39.02 0.92	0.53	0.70	36.60	39.51	1048
	m = 5 0.93	0.64	0.81	36.86	37.72 0.92	0.65	0.83	37.62	38.53	582

Table 5 .

 5 Summary statistics of the multiple linear regression models for all objective space dimensions (∀m) and for each individual objective space dimension m ∈ {2, 3, 5}. Only the subset of significant features (k, m, ρ, hv, log(nplo)) are part of the models. Values are rounded to 10 -2 .

				training set			10-fold cross validation
		r MAE RMSE %RAE %RRSE	r MAE RMSE %RAE %RRSE #inst
	∀m	0.91	0.58	0.77	37.95	40.66 0.91	0.58	0.77	38.00	40.71	2980
	m = 2 0.92	0.53	0.69	36.14	39.17 0.92	0.53	0.69	36.25	39.29	1350
	m = 3 0.92	0.54	0.71	37.35	39.89 0.92	0.55	0.71	37.62	40.14	1048
	m = 5 0.92	0.64	0.82	37.30	38.23 0.92	0.65	0.83	37.54	38.52	582

Acknowledgements. This work was partially supported by the Japanese-French research project "Global Research on the Framework of Evolutionary Solution Search to Accelerate Innovation" (JSPS, 2013-2016).