Estefania Yap
email: estefaniay@student.unimelb.edu.au

Mario A Muñoz
email: munoz.m@unimelb.edu.au

Kate Smith-Miles
email: smith-miles@unimelb.edu.au

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

Instance Space Analysis of Combinatorial Multi-objective Optimization Problems

Keywords: Multi-objective optimization, black-box combinatorial optimization, landscape analysis, feature-based performance prediction

In recent years, there has been a continuous stream of development in evolutionary multi-objective optimization (EMO) algorithms. The large quantity of existing algorithms introduces difficulty in selecting suitable algorithms for a given problem instance. In this paper, we perform instance space analysis on discrete multi-objective optimization problems (MOPs) for the first time under three different conditions. We create visualizations of the relationship between problem instances and algorithm performance for instance features previously identified using decision trees, as well an independent feature selection. The suitability of these features in discriminating between algorithm performance and understanding strengths and weaknesses is investigated. Furthermore, we explore the impact of various definitions of "good" performance. The visualization of the instance space provides an alternative method of algorithm discrimination by showing clusters of instances where algorithms perform well across the instance space. We validate the suitability of existing features and identify opportunities for future development.

I. INTRODUCTION

Many real world problems consist of two or more, often conflicting, objectives. These are known as multi-objective optimization problems (MOPs), whose solution consists of a set of optimal trade-offs between the objectives. Evolutionary Multi-objective Optimization (EMO) algorithms are used to solve MOPs, due to their ability to maintain a diverse population of solutions spanning the trade-offs in a single simulation run.

Since the pioneering work of Schaffer [START_REF] Schaffer | Some experiments in machine learning using vector evaluated genetic algorithms (artificial intelligence, optimization, adaptation, pattern recognition)[END_REF], there has been substantial interest in the development of EMO algorithms. The now large number of available algorithms leads to the question of which one to choose, since the quality of the solution is affected by the choice of algorithm. This challenge is known as the algorithm selection problem [START_REF] Rice | The algorithm selection problem[END_REF] which aims to predict which algorithm will perform well with minimized cost to solve a given problem, where cost is measured in the resources used such as time. However, understanding the relationship between algorithm performance and problem instances is not a simple task. This is because problem difficulty is affected by many factors, such as multi-modality or ruggedness [START_REF] Malan | A survey of techniques for characterising fitness landscapes and some possible ways forward[END_REF]. If the optimization problem can be stated analytically, Funding was provided by the Australian Research Council through the Australian Laureate Fellowship FL140100012. such factors can be measured, but typical MOPs are black-box problems described only by sample input (decision variables) and outputs (objective function values). Given that no a priori information is available for black-box problems, how can we identify the relationship between instance features and algorithm performance?

The fitness landscape of black-box optimization problem instances offers an analysis framework, extracting features that characterize instance topology. These features quantify the challenges encountered by algorithms, independent of problem classes. Landscape analysis in single-objective optimization has receiving growing interest [START_REF] Mersmann | Benchmarking evolutionary algorithms: Towards exploratory landscape analysis[END_REF], with new features being developed. However, while they can be extended to multiobjective optimization, there may be features that are not relevant due to additional challenges incurred by the conflicting nature of the objectives. To counter this, features that are multi-objective in nature are necessary to characterize these problems. Such research exists [START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF]- [START_REF] Liefooghe | What makes an instance difficult for black-box 0-1 evolutionary multiobjective optimizers?[END_REF] but many require enumeration of all Pareto optimal solutions and are therefore impractical for the algorithm selection problem which is attempting to predict performance based only on instance data, not algorithm performance data.

In recent work, Liefooghe et al. [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] introduce several new cheap features using short random and adaptive walks. This work showed significant features that helped discriminate between the performance of three algorithms (IBEA [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], MOEA/D [START_REF] Zhang | Moea/d: A multiobjective evolutionary algorithm based on decomposition[END_REF], NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]) on a range of problem instances using ⇢mnk landscapes. Liefooghe et al. used a CART decision tree to provide an easily interpretative model to distinguish which algorithm is likely to be best for a given instance. A more complex random forest model [START_REF] Geurts | Extremely randomized trees[END_REF] was also explored.

In this paper, we build upon this work by carrying out the first Instance Space Analysis (ISA) [START_REF] Smith-Miles | Towards objective measures of algorithm performance across instance space[END_REF] of combinatorial multi-objective problems using the same instances, features and algorithms as [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF]. Our aim is twofold: to create a visual complement using the CART decision tree features, and explore how the strengths and weaknesses of algorithms depend on the chosen criteria for good performance (i.e. best performing algorithm, or close to best). ISA also provides insights into how features change across the space and their correlation with algorithm performance. As such, we are able to generate a mapping of each algorithm's "footprint," i.e., the regions of the space where it has strengths or weaknesses. Support vector machines (SVMs) are then used to decide which algorithm is recommended across various regions of the instance space, based on a set of calculated features. Through this analysis, we also assess the diversity of the ⇢mnk problem instances, and identify sparse areas in the instance space. Our analysis is comprised of three stages:

1) visually supporting conclusions of [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] -performing the first instance space analysis on the CART subset of features to see if similar rules can be visualized, 2) extending explanations -analysis using feature selection methods to consider additional features that may provide further insights, and 3) relaxing performance -exploring the impact when the definition of "good" is relaxed from best hyper-volume to be within 1% of best hyper-volume. The remainder of this paper is organized as follows. In Section II, we introduce the general multi-objective combinatorial optimization problem, and the benchmark instances used for evaluation. We then describe the features that are utilized in this paper. In Section III we describe the instance space methodology for generating a 2D visualization of algorithm performance across instances. Results are presented in Section IV where we validate the relationship between features and algorithm performance in the previous paper [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF]. Alternative instance spaces generated by different features are also presented, and we investigate the impacts of relaxed conditions of "good" performance to obtain stronger insights into algorithm strengths and weaknesses. Section V concludes the paper with a discussion of future research directions.

II. BACKGROUND

A. Definition

Without loss of generality, consider a multi-objective minimization problem with m objectives and n decision variables:

min y = f (x) = (f 1 (x) , ..., f m (x)) where x = (x 1 , ..., x n ) 2 X y = (y 1 , ..., y m ) 2 Y
where x is known as the decision vector, X is the search space, y the objective vector and Y the objective or fitness space, and f (x) the fitness function. In multi-objective combinatorial optimization, X is a discrete set, i.e. X := {0, 1}

n . Given two solutions x (1) and x (2) , x (1) dominates x (2) (denoted as x (1) x (2) ) if it is better in at least one objective, and at least as good in all others. Solution x (1) is said to be non-dominated with respect to x (2) , and x (2) is dominated with respect to x (1) . This relation then holds for the objective vectors, with y (1) is non-dominated with respect to y (2) . A solution x ? is Pareto optimal (PO) if there is no solution x 2 X that dominates it. Similarly, y ? 2 Y is non-dominated if there does not exist any y 2 Y such that y ? is dominated by y. The set of PO solutions forms the Pareto set, and their corresponding mapping into the objective space is known as the Pareto front (PF). Generating a good approximation of the PF is the main goal in solving MOPs.

B. Benchmark Instances

The selection of a representative subset of MOP instances is required in order to study the problem domain. For the purposes of the study, ⇢mnk-landscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space: Mnk-landscapes with correlated objectives[END_REF] are used to construct the test instances for combinatorial MOPs. This is due in part to their ability to control properties such as multimodality and objective correlation. A ⇢mnk-landscape can be formally defined as:

max f i (x) = 1 n n X j=1 f ij (x j , x j1 , . . . x j k ) i 2 {1, . . . , m} s.t. x j 2 {0, 1} j 2 {1, . . . , n}
where n is the number of bits in the decision space, m the number of objectives and k the number of epistatic interactions influencing the contribution of each x j . The fitness f i (x) of a solution x 2 X corresponds to the average value of the contributions in its n fitness components. Problem instances can be tuned from smooth to rugged by increasing k from 0 to (n 1). The parameter ⇢ > 1 m 1 controls the degree of correlation between the objectives.

C. Features

The purpose of features is to quantify structural characteristics about a problem instance. They are used to summarize instances, and to identify the cause of difficulty. Features are used for gaining insights into the relationship between algorithms and instances, and consequently enabling suitable algorithm selection using performance prediction.

Multi-objective landscapes require the inclusion of multiple PO solutions, due to the trade-offs resulting from conflicting objectives. Similarly to the single-objective case, characteristics such as multimodality and ruggedness are empirically linked to instance difficulty and algorithm performance [START_REF] Weinberger | Correlated and uncorrelated fitness landscapes and how to tell the difference[END_REF], [START_REF] Kerschke | Towards analyzing multimodality of continuous multiobjective landscapes[END_REF]. Therefore, it is desirable to identify features that correlate well with these characteristics. In particular, useful features should be captured using the fewest function evaluations.

Liefooghe et al. [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] introduced local features which compute landscape information by using random or adaptive walks to sample the neighbourhood of solutions. The neighbourhood relation is defined as N : X 7 ! 2 X , and is here based on the 1-bit-flip operator. In a random walk, a random neighbour is selected at each step. During an adaptive walk, an improving (i.e., dominating) neighbour is selected at each step using a Pareto hill climber (PHC) [START_REF] Verel | On the structure of multiobjective combinatorial search space: Mnk-landscapes with correlated objectives[END_REF].

During a walk of size `, the neighbourhood N of each sample is explored. The proportion of dominated (#inf), dominating (#sup), incomparable (#inc), non-dominated (#lnd), and supported (#lsupp) solutions in the neighbourhood are measured. Features relating to the hyper-volume are also calculated. These include the hyper-volume covered by each For samples collected during random and adaptive walks, features are measured and averaged. Additionally, the first autocorrelation coefficient of random walks is known to characterize the ruggedness of the single-objective landscapes [START_REF] Weinberger | Correlated and uncorrelated fitness landscapes and how to tell the difference[END_REF]. Therefore, to accommodate this in the multi-objective case, the first autocorrelation coefficient of the above features during random walks is also measured. Lastly, the correlation between objectives is measured during the random walk. The full list of 26 features obtained during random and adaptive walks can be found in Table I.

D. Prior Work

Several new local features defined above were introduced in the work of Liefooghe et al. [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] with the intention of characterizing the structure of combinatorial MOPs. Tree-based predictive models were used to highlight the differences in important features between local and global search heuristics on small instances.

Large scale experiments were also performed for featurebased performance prediction. We utilize the same experimental setup in this paper. A total of 1,000 ⇢mnk-landscapes are created using latin hypercube sampling. The parameters of these instances are generated within the domains of n 2 {64, . . . , 256}, k 2 {0, . . . , 8, }, m 2 {2, . . . , 5, } and

⇢ 2 [ 1, 1] with ⇢ > 1 m 1 .
The local features in Section II-C are measured, as well as the values of m and n, as they are known parameters in black-box settings.

Three state-of-the-art EMO algorithms were chosen as representatives of different evolutionary approaches. The indicator based method is represented by IBEA [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF], scalarization by MOEA/D [START_REF] Zhang | Moea/d: A multiobjective evolutionary algorithm based on decomposition[END_REF], and dominance by NSGA-II [START_REF] Deb | A fast and elitist multiobjective genetic algorithm: NSGA-II[END_REF]. These are run using default parameters in jMetal 4.5 [START_REF] Durillo | jmetal: A java framework for multi-objective optimization[END_REF]. All algorithms use a population size of 100, 2-point crossover with a rate of 0.9, and a bit-flip mutation rate of 1/n. The budget is set at 10 6 evaluations. Performance is measured by the hypervolume metric, with "good" initially defined as the best hypervolume. Later we relax this goodness definition to be within 1% of best hyper-volume.

A CART decision tree was used in [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF] to distinguish between the preferred algorithm for an instance, based on measured features. The features selected in the decision tree included lnd avg aws, hv r1 rws, hvd avg rws, lnd avg aws and hv avg aws. A randomized trees model was used to generate the final output, with an error rate of 0.014 for the best statistical rank.

III. METHODOLOGY

Instance Space Analysis (ISA) is a methodology that allows the visualization of a group of problem instances in a 2D space. Within the instance space it is possible to observe trends in algorithm performance; hence, we are able to identify and objectively measure the regions of strength (and conversely, weaknesses) for an algorithm, which are known as footprints. Moreover, ISA gives us information about the distribution of the instances, allowing us to identify sparse regions where new instances may be needed. First described by Smith-Miles et al. [START_REF] Smith-Miles | Towards objective measures of algorithm performance across instance space[END_REF], ISA can be automatically performed through publicly available webtools [START_REF] Smith-Miles | Melbourne algorithm test instance library with data analytics (MATILDA)[END_REF] or using the MATLAB toolkit [START_REF] Muñoz | Instance space analysis: A toolkit for the assessment of algorithmic power[END_REF]. In general, ISA involves the following steps:

1) collecting the meta-data, which corresponds to a table

where each row represents a test instance and each column is a feature or the performance of an algorithm; 2) selecting a subset of features which best distinguish between different algorithm performances, and variations in problem instances; 3) projecting onto a 2D space for visualization; 4) measuring the algorithm footprints; and 5) generating new test instances in sparse regions.

In this study, we will not be generating new test instances, and instead focus on the analysis of the algorithm footprints. Details on all the steps of the ISA methodology can be found in previous work [START_REF] Kandanaarachchi | Instance space analysis for unsupervised outlier detection[END_REF], [START_REF] Muñoz | Instance spaces for machine learning classification[END_REF]. However, we briefly summarise them here for reference. Once the meta-data has been collected, a binary performance measure is calculated that defines when an algorithm is "good." As mentioned above, "good" is here initially defined as largest hyper-volume or, later, as being within 1%. Then, the features are pre-processed by bounding the outliers within the mean plus or minus five times the interquartile range. Next, both the features and performance measures are normalized and standardized using Box-Cox and Z transformations.

We filter out those features that are not highly predictive of algorithm performance, by calculating the absolute value of the Pearson correlation between them. Then, we cluster these meta-features using as similarity metric 1 |⇢ x,y |, where ⇢ x,y is the correlation between two features. Silhouette analysis is used to determine the number of clusters. Retaining one meta-feature per cluster to reduce redundancy, we consider all valid combinations. Using PCA with two components as dimensionality reduction, and Random Forests (RF) to predict the binary measure, we determine which projection has the highest out-of-the-bag average accuracy. We use PCA and RFs as they are cheaper alternatives than the algorithms used in the following stages. This approach ensures that our selected meta-features are predictive of algorithm performance across the portfolio.

Next, we use the Prediction Based Linear Dimensionality Reduction (PBLDR) method [START_REF] Muñoz | Instance spaces for machine learning classification[END_REF] to find the final 2D projection, which creates the most linear trends of algorithm performance and feature values across the instance space, to assist visualization of directions of hardness and feature correlations. BFGS is used to solve numerically the underlying optimization problem in PBLDR, which is known to have infinite solutions. As such, we calculate 30 different projections and select the one with the highest topological preservation, defined as the correlation between high-and low-dimensional distances.

The instance space analysis assists in guiding automatic algorithm selection by highlighting which algorithms are compatible within regions of the space. For new instances, features are measured and projected onto the instance space. The location of this instance and whether it lies within any algorithm's footprint can be used to guide selection. We repeat this analysis for each of the three stages mentioned in Section I: (A) supporting conclusions drawn using the CART subset, (B) extending explanations through the use of the feature selection method used above and (C) relaxing performance by defining "good" performance as within 1% of best hyper-volume.

IV. RESULTS

A. Supporting Conclusions

The final projection matrix of the generated instance space using the CART features, shown in Fig. 2, is defined by the following linear projection to 2D. Each instance is a point in the 2D space shown in Fig 2,defined

While MOEA/D is the algorithm that performs best on average, good performance is almost evenly split between NSGA-II, as seen in Fig. 1b and1c. The footprints for each algorithm are mostly clustered. Fig. 2 shows hv avg aws has a good correlation with the splitting of spaces in which each algorithm is superior. NSGA-II performs well on instances with a larger value, while the opposite is true for MOEA/D. This is a similar recommendation given by the CART model. This feature is known to strongly negatively correlate with the number of objectives m [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF], suggesting that MOEA/D is preferred for more objectives. Thus, our approach shows an alternative approach to the CART model and generates a visual representation of performance. In all generated instance spaces, it is evident that IBEA only performs well on a subset of the instances that NSGA-II performs well on. Fig. 7 shows that it is never recommended by the SVM as the preferred algorithm. The SVM in Fig. 7a shows that MOEA/D and NSGA-II are each selected as the preferred algorithm in approximately half of the instances.

B. Extending Explanations

When using feature selection on the complete set, only two features were shared with those in the CART subset (hv avg aws, lnd avg rws) shown in Fig. 4. The final projection matrix of the generated instance space in Fig. 4 is defined by: We observe the negative correlation between m and hv avg aws, with MOEA/D's footprint occupying instances with larger m. The performance of MOEA/D on larger values of sup avg rws is also observed. This feature is correlated with the objective correlation ⇢ [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF], suggesting that MOEA/D performs well when ⇢ is larger, i.e. when the degree of conflict between the objectives is small.

Z f eatselect = 2
The subset of features selected shows a similar partitioning of the instance space between MOEA/D and NSGA-II's performance, as shown in Fig. 3b and3c. However, the projection seems to favour NSGA-II slightly, showing a clearer partition of best performance. Here, the feature hv avg aws correlates with algorithm performance, as do inc avg rws and lnd avg rws. As a result, these features are able to generate a useful instance space for inferring algorithm performance.

The projection of the instance space here is more advan- tageous for NSGA-II, as seen earlier in its footprint. This is visible in Fig. 7b for the SVM, showing preference on NSGA-II for more than half of the instances.

C. Relaxing Performance

If we now relax the definition of good performance to include algorithms whose hyper-volume is within 1% of the best, a different set of features are selected to explain performance. The final projection matrix of the generated instance space in Fig. 6, can be defined by:

Z relaxedperf = 2
The feature hv avg aws is once again selected as a good predictor for algorithm performance.

In relaxed conditions of "good", Fig. 5c shows that the space where NSGA-II displays good performance is more cleanly partitioned. However, MOEA/D performs so strongly that its footprint overlaps with that of NSGA-II's, without any clear separation. While the earlier analyses showed the two algorithms complementing each other, here the performance of MOEA/D is superior. Due to the shared good performance on much of the space, the SVM has difficulties in identifying clustered regions within the instance space where MOEA/D performs well. Therefore, the features currently used may

D. Algorithm Selection Results

The test instances for all analyses cover a large portion of the instance spaces generated with minor empty pockets closer to the origin, as seen in Fig. 7. We conclude that the set of test instances studied are quite diverse and span the instance space sufficiently to infer performance behaviours across the instance space.

V. CONCLUSIONS AND FUTURE WORK

In this paper, we carried out the first Instance Space Analysis (ISA) of combinatorial MOPs, using as a basis the work by Liefooghe et al. [START_REF] Liefooghe | Landscape-aware performance prediction for evolutionary multi-objective optimization[END_REF]. ISA complements that of the decision tree analysis, offering additional insights into algorithm strengths and weaknesses through visualizations. The development of such insights into the transparency of algorithm performance allows for better automated algorithm selection for new, untested instances. When an untested instance is encountered, features are evaluated mapped into their instance space co-ordinates. Given these co-ordinates, the suitable algorithms can be inferred by selecting those which show dominant performance over the region. Moving forward, there is the opportunity to extend the approach to consider a broader range of algorithms and instance features to explore a more comprehensive understanding of algorithm behaviours. Investigation of the algorithm footprints showed that MOEA/D displays strong performance under relaxed conditions of goodness, as seen in Fig. 5. However, its weaknesses are sparsely located. Therefore, despite MOEA/D being the higher performing algorithm, the SVM selected NSGA-II more often due to NSGA-II's consistent region of good performance. While these features may be sufficient in selecting the best algorithm by hyper-volume, when the condition of goodness is relaxed, they may not provide enough information to separate algorithm performance. We therefore focus our attention toward developing features that hold more discrimination between algorithms in future, as well as expanding our algorithm portfolio. Furthermore, we are interested in developing features that can be calculated cheaply to minimize the total cost of function evaluations. The amalgamation of this work will provide further insights, assisting with better informed algorithm selection in future in the multi-objective space.

Fig. 1 :

 1 Fig.1: Algorithm footprints in the projected instance space of (1) using the CART subset of features.

Fig. 2 :

 2 Fig. 2: Distribution of normalized CART subset of features in the projected instance space of (1). The color scale corresponds to normalized feature values.

Fig. 3 :

 3 Fig. 3: Algorithm footprints in the projected instance space of (2) using feature selection.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Distribution of normalized subset of features in the projected instance space of (2).

Fig. 6 :

 6 Fig. 6: Distribution of normalized subset of features in the projected instance space of (3) when performance is relaxed.

TABLE I :

 I Features collected via sampling methods, used in this paper.

	Description	Random walk	Adaptive walk
	estimated correlation between objectives		f cor rws
	length of walk			length aws
		First autocorrelation Average	Average
	proportion of dominated solutions	#inf r1 rws	#inf avg rws	#inf avg aws
	proportion of dominating solutions	#sup r1 rws	#sup avg rws	#sup avg aws
	proportion of incomparable solutions	#inc r1 rws	#inc avg rws	#inc avg aws
	proportion of non-dominated solutions	#lnd r1 rws	#lnd avg rws	#lnd avg aws
	proportion of supported solutions	#lsupp r1 rws	#lsupp avg rws #lsupp avg aws
	hyper-volume covered by a single solution	hv r1 rws	hv avg rws	hv avg aws
	difference of hyper-volume covered by a single solution and its neighbourhood	hvd r1 rws	hvd avg rws	hvd avg aws
	hyper-volume covered by the entire neighbourhood	nhv r1 rws	nhv avg rws	nhv avg aws
	solution (hv), the difference between the hyper-volume cov-
	ered by the neighbour and current solution (hvd), and the
	hyper-volume covered by the whole neighbourhood (nhv).

  by:

		2	0.1628	0.6232	3	> 2	lnd avg rws	3
	Z CART =	6 6 6 4	0.4698 0.4443 0.3854	0.1364 0.079 0.048	7 7 7 5	6 6 6 4	hv r1 rws hvd avg rws lnd avg aws	7 7 7 5
			0.249	0.7237			hv avg aws


Omar Abdelkafi
email: omar.abdelkafi@univ-lille.fr

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

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

Darrell Whitley
email: whitley@cs.colostate.edu

On the Design of a Partition Crossover for the Quadratic Assignment Problem

We conduct a study on the design of a partition crossover for the QAP. On the basis of a bipartite graph representation, we propose to recombine the unshared components from parents, while enabling their fast evaluation using a preprocessing step for objective function decomposition. Besides a formal description and complexity analysis of the proposed crossover, we conduct an empirical analysis on its relative behavior using a number of large-size QAP instances, and a number of baseline crossovers. The proposed operator is shown to have a relatively high intensification ability, while keeping execution time relatively low.

Introduction

One of the key ingredients in the success of evolutionary algorithms is the design of effective and efficient crossover operators. In the context of gray-box optimization, problem-specific properties can help in designing dedicated operators. This is the case of partition crossovers, developed for a number of combinatorial problems (e.g., TSP [START_REF] Tinós | Generalized asymmetric partition crossover (gapx) for the asymmetric tsp[END_REF], SAT [START_REF] Chicano | Enhancing partition crossover with articulation points analysis[END_REF], NK-landscapes [START_REF] Tinós | Partition crossover for pseudo-boolean optimization[END_REF]), and allowing to efficiently explore large search spaces. The term partition crossover is here used in a general sense, to render the idea of decomposing the variables according to their values in the parents, and then recombining them in such a way that the best improving offspring can be computed efficiently. In other words, a partition crossover is based on the idea of optimal recombination of two parents, which requires to find the best possible offspring among all possible ones, while fulfilling the genotype inheritance principle. As such, two issues must be considered.

Firstly, one should specify the recombination mechanism allowing gene transmission while fully exploring the improvement potential of parents. In particular, the shared genes from parents are kept identical and the other genes are to be properly mixed. Secondly, since the set of possible offspring underlying such a process is typically huge, one must rely on some properties to compute the best offspring in the most efficient manner. For instance, the so-called k-bounded binary problems [START_REF] Tinós | Partition crossover for pseudo-boolean optimization[END_REF] can be decomposed as a linear combination of sub-functions of at most k variables. This guarantees that the contributions of non interacting variables to the global fitness are additive, and hence, the choice of the optimal gene sequence when performing recombination can be performed greedily in an efficient manner. Another example is TSP [START_REF] Tinós | Generalized asymmetric partition crossover (gapx) for the asymmetric tsp[END_REF], where the cost of the global tour is additive with respect to the length of different sub-paths sharing the same cities, but in different order, in the parent solutions.

In this paper, we are interested in the Quadratic Assignment Problem (QAP) [START_REF] Koopmans | Assignment problems and the location of economic activities[END_REF][START_REF] Loiola | A survey for the quadratic assignment problem[END_REF][START_REF] Misevicius | Comparison of crossover operators for the quadratic assignment problem[END_REF][START_REF] Sahni | P-complete approximation problems[END_REF], for which no partition crossover has been developed so far, despite the broad range of dedicated algorithms. A main difficulty comes from the definition of its objective function, where the contribution of every single variable is sensitive to all other variables. Hence, it becomes a challenging issue to design a decomposition process which enables both parent recombination and fast offspring evaluation. This is precisely the aim of our work, and our contribution is to be viewed as a first step towards the design of an efficient and effective partition crossover operator for the QAP. More specifically, we rely on an intuitive bipartite graph representation for decomposition. We then show how recombination and offspring evaluation can be performed on that basis. Moreover, we conduct an empirical study rendering the performance of the designed crossover comparatively to existing ones, either by itself, or when combined with a fast local search process. Our analysis shows that the proposed crossover has a high intensification power while keeping execution time relatively low.

The paper is organized as follows. In Section 2, we introduce the QAP and review some existing crossovers. In Section 3, we describe our main contribution towards the design a partition crossover for the QAP. In Section 4, we report our experimental findings. In Section 5, we conclude the paper.

Background

Problem Definition

The Quadratic Assignment Problem (QAP) [START_REF] Koopmans | Assignment problems and the location of economic activities[END_REF][START_REF] Sahni | P-complete approximation problems[END_REF] aims at assigning n facilities I, to n locations J. Let f hi be the flow between facilities h and i, and d sj be the distance between locations s and j. The objective is to minimize the sum of the products between flows and distances, such that each facility is assigned to exactly one location. The solution space can be defined as the set Π of permutations of {1, . . . , n}. Given a permutation solution π, the i th element π(i) corresponds to assigning facility i to location π(i). The QAP is then stated as:

arg min π∈Π h∈I i∈I d π(h)π(i) f hi (1)
The QAP is NP-hard [START_REF] Sahni | P-complete approximation problems[END_REF], and is considered as one of the most difficult problems from combinatorial optimization. As such, we have to rely on heuristic approaches such as stochastic local search and evolutionary algorithms [START_REF] Loiola | A survey for the quadratic assignment problem[END_REF].

Representative Crossover Operators

In this paper, we are interested in designing efficient and effective crossover operators. Among the large number of hybrid genetic algorithms for the QAP; see, e.g., [START_REF] Lim | Efficient genetic algorithms using simple genes exchange localsearch policy for the quadratic assignment problem[END_REF][START_REF] Merz | Fitness landscape analysis and memetic algorithms for the quadratic assignment problem[END_REF][START_REF] Misevicius | An improved hybrid genetic algorithm: new results for the quadratic assignment problem[END_REF][START_REF] Vázquez | A hybrid genetic algorithm for the quadratic assignment problem[END_REF]], crossover appears to be a crucial component. In this respect, we can distinguish two families of crossovers [START_REF] Misevicius | Comparison of crossover operators for the quadratic assignment problem[END_REF]. In half-structured crossovers, the offspring preserve only a part of the parent genes, whereas the remaining part is typically generated at random. In fully-structured crossovers, the offspring genes are obtained by preserving the ones from parents. In our work, we consider four usual and representative baseline crossovers, two from the first family, and two from the second one. We first describe the half-structured crossovers, namely Opx and Ux, and then the fully structured ones, namely Cx and Spx. Notice that these crossovers can also find applications in other problems such as TSP [START_REF] Larrañaga | Genetic algorithms for the travelling salesman problem: A review of representations and operators[END_REF].

The Opx crossover. This is a standard one-point crossover, e.g., [START_REF] Lim | Efficient genetic algorithms using simple genes exchange localsearch policy for the quadratic assignment problem[END_REF]. Given two parent permutations, a random point is selected to define two parts for each parent. A new offspring is obtained by first preserving the locations from the first parent up to the chosen random point. The remaining part is filled by copying the elements from the second part of the second parent, excluding those that were already copied from the first parent. This may lead to the situation where some facilities are not assigned to any locations. The offspring is hence complemented at random, using the remaining locations that were not yet included.

The Ux crossover. This is a standard uniform crossover operator, e.g., [START_REF] Tate | A genetic approach to the quadratic assignment problem[END_REF]. Some locations are selected at random following a Bernoulli distribution with parameter 1/2. The locations occupied by the selected facilities in the first parent are copied. The locations occupied by the remaining facilities in the second parent are also copied unless they were already included from the first parent. The offspring is complemented by randomly assigning the missing locations.

The Cx crossover. This is the so-called cycle crossover [START_REF] Merz | A comparison of memetic algorithms, tabu search, and ant colonies for the quadratic assignment problem[END_REF][START_REF] Merz | Fitness landscape analysis and memetic algorithms for the quadratic assignment problem[END_REF][START_REF] Oliver | A study of permutation crossover operators on the traveling salesman problem[END_REF]. For clarity, we consider the example of Fig. 1a. All shared assignments are copied, i.e., facilities 1 and 7 assigned to locations 5 and 9. The crossover starts iterating from the first different facility assignment, which is facility 2 assigned to location 3 in the first parent. Looking at the second parent, facility 2 is assigned to location 8. This location is occupied by facility 3 in the first parent. Similarly, facility 3 is assigned to location 4 in the second parent, and so on. A cycle -alternating between the same set of locations in a different order on both parents -is then detected when arriving to facility 8 assigned to location 3 in the second parent. Hence, one parent is selected at random and the so-computed locations are preserved in the offspring. This procedure is repeated until all facilities are assigned.

The Spx crossover. This is the so-called Swap Path Crossover [START_REF] Ahuja | A greedy genetic algorithm for the quadratic assignment problem[END_REF][START_REF] Glover | Genetic algorithms and scatter search: unsuspected potentials[END_REF]. The crossover is based on iteratively swapping unshared locations. In Fig. 1b, it starts with facility 2 assigned to location 3 (resp. 6) in the first (resp. second) parent. In the first (resp. second) parent, location 6 (resp. 3) is occupied by facility 5 (resp. 8). Hence, a swap is performed between locations 2 and 5 in the first parent, as well as between locations 2 and 8 in the second parent. Two offspring are obtained, respectively to the first and second parent. Then, the best of the two replaces its corresponding parent, say the first one as in Fig 1b . Notice that after this iteration, the two new parents have one more location in common.

The same is then repeated iteratively until both parents become the same. The output offspring is the best ever created during all iterations.

A Partition Crossover for the QAP

A partition crossover is based on the idea of decomposing the evaluation function as well as the set of variables during recombination. This decomposition should enable to construct and to evaluate offspring using partial evaluations.

The fastest the evaluation step, the fastest the exploration of a large number of offspring, eventually leading to an improving one, hence exploring the socalled dynastic potential of parents at best. It should be clear that none of the previously-described crossovers fulfill this requirement, although Cx and Spx attempt to construct an offspring based on the idea of gene transmission.

In the following, we aim at designing a new partition crossover for the QAP. We start by introducing some notations, then we describe a decomposition process for the QAP, and the underlying recombination (Proposition 1) and function evaluation process (Corollary 1), in a formal, but intuitive, manner.

Recombination based on a Bipartite Graph Representation

Following standard notations from graph theory, a bipartite graph G is a graph whose nodes can be divided into two disjoint sets U (G) and V (G), such that an edge in E(G) can only connect a node in U to one in V . Given a subset of nodes

Q ⊆ U and R ⊆ V , we denote by G[Q, R] the subgraph induced by Q ∪ R in G.
Let π be a permutation solution for the QAP. Let us define the bipartite graph whose nodes sets are respectively the set of facilities I and locations J, and where every facility i is connected to its unique location j according to permutation π. This corresponds in fact to a very natural representation of a feasible assignment of the facilities to the locations; see Fig. 2.

Definition 1. Let G π = (I, J, E) be the bipartite graph such that E(G π ) = {(i, j) | π(i) = j}.
Let us now consider two permutation solutions π 1 and π 2 which will play the role of parents for our target crossover. Let I (resp., I) be the set of facilities that are assigned (resp., not assigned) to the same locations, denoted J (resp., J).

Definition 2. Let I = {i ∈ I | π 1 (i) = π 2 (i)} and I = I \ I. Then, let J = {π 1 (i) | i ∈ I}, and J = J \ J.
Notice that for any facility i ∈ I, the corresponding locations π 1 (i) and π 2 (i) are different, and both belong to J. Let us now define the bipartite graph G π1π2 obtained by merging the edges of G π1 and G π2 ; see Fig. 2.

C 3 in G π1π2 [I, J], with U (C 1 ) = {1, 3, 4}, U (C 2 ) = {5, 6}, and U (C 3 ) = {7, 9}. Definition 3. Let G π1π2 = (I, J, E(G π1 ) ∪ E(G π2 ))
Let us focus on the connected components of the so-obtained graph. First, we have that every facility i ∈ I is connected, by exactly two edges in G π1π2 , to a unique location j = π 1 (i) = π 2 (i) ∈ J, and j is not connected to any other facility. Hence, this implies exactly one connected component with these two nodes connected by two parallel edges. Apart from such components (in I ∪ J), the other components of interest connect nodes in I to nodes in J. Then, Definition 4. Let k be the number of connected components in G π1π2 [I, J], and let

C = {C 1 , C 2 , • • • , C k } be the set of these connected components.
Notice that k = 0 iff π 1 = π 2 . In the following, we assume that π 1 = π 2 , and hence k ≥ 1. By definition, we also have that for every ∈ {1, . . . , k}, C is a bipartite graph whose edges form a cycle that are alternatively in E(G π1 ) and in E(G π2 ). An offspring can hence be constructed by: (i) preserving the shared assignments of facilities in I, and (ii) choosing for every set of facilities U (C ), ∈ {1, . . . , k}, the locations they are connected to either in G π1 (first parent) or in G π2 (second parent). More formally, let m : {1, . . . , k} → {1, 2} be an arbitrary mapping function. Such a mapping can be used to decide which of π 1 or π 2 to consider when choosing the locations of U (C ) for every , i.e., if m( ) = 1 then use π 1 , otherwise use π 2 . This is stated in the following proposition summarizing the proposed crossover recombination mechanism. Proposition 1. Given two permutations π 1 and π 2 and an arbitraty mapping function m : {1, . . . , k} → {1, 2}, ⊗ π as defined in the following is a feasible permutation offspring.

• ∀i ∈ I, ⊗ π(i) = π 1 (i) = π 2 (i) • ∀i ∈ I, ⊗ π(i) = π m( i) (i), where i ∈ {1, . . . , k} is such that i ∈ U (C i )
The previous proposition is to recall the Cx crossover, where a cycle corresponds to a connected component in our formalism. However, we consider to explore not solely one random offspring, but the whole set of offspring solutions that can be constructed by Proposition 1 to find the best possible one. Since, the recombination process is fully and uniquely determined by the choice of the mapping function m, we have to consider all possible mappings. Obviously, there exist 2 k possibilities (including parents) leading to as much possible offspring. Since evaluating one offspring from scratch costs Θ(n 2 ), a naive approach to compute the fitness values of these possible offspring has a complexity of Θ(2 k • n 2 ), which can be prohibitive. In the next section, we show how to reduce it.

Decomposition of the Evaluation Function

Following the previous notations, let m : {1, . . . , k} → {1, 2} be an arbitrary mapping function, and let us define for every and in {1, . . . , k}:

Q = h∈U (C ) i∈U (C ) d π m( ) (h)π m( ) (i) f hi ; Q = h∈ I i∈U (C ) d π1(h)π m( ) (i) f hi R = h∈ I i∈ I d π1(h)π1(i) f hi ; Q = i∈ I h∈U (C ) d π m( ) (h)π1(i) f hi
where Q represents the contribution of facilities in C w.r.t. the facilities in another connected component C . Similarly, Q (resp. Q ) represents the contribution of facilities in C w.r.t. facilities having the same locations in both parents π 1 and π 2 . Finally, R represents the pairwise contribution of the facilities that have the same locations in π 1 and π 2 . The following proposition shows that the QAP objective function can be decomposed using these contributions.

Proposition 2. Let ⊗ π be a permutation offspring as defined in Proposition 1. Then,

f ( ⊗ π) = R + k =1 Q + Q + k =1 Q Proof.
By decomposing the facilities w.r.t. I, I, we get:

f ( ⊗ π) = h∈I i∈I d⊗ π(h) ⊗ π(i) f hi = h∈I i∈I d⊗ π(h) ⊗ π(i) f hi + h∈I i∈ I d⊗ π(h) ⊗ π(i) f hi = h∈I k =1 i∈U (C ) d⊗ π(h) ⊗ π(i) f hi + h∈ I k =1 i∈U (C ) d⊗ π(h) ⊗ π(i) f hi + h∈I i∈ I d⊗ π(h) ⊗ π(i) f hi + h∈ I i∈ I d⊗ π(h) ⊗ π(i) f hi = k =1 h∈U (C ) k =1 i∈U (C ) d⊗ π(h) ⊗ π(i) f hi + h∈ I k =1 i∈U (C ) d⊗ π(h) ⊗ π(i) f hi + k =1 h∈U (C ) i∈ I d⊗ π(h) ⊗ π(i) f hi + h∈ I i∈ I d π1(h)π1(i) f hi = k =1 k =1 h∈U (C ) i∈U (C ) d π m( ) (h)π m(l) (i) f hi + k =1 h∈ I i∈U (C ) d π1(h)π m( ) (i) f hi + k =1 i∈ I h∈U (C ) d π m( ) (h)π1(i) f hi + R = k =1 k =1 Q + k =1 Q + k =1 Q + R
As a result, we obtain the following corollary which follows from the fact that the contributions appearing in the decomposition of Proposition 2 can only have a constant number of values, for all possible choices of the mapping m.

Experimental Analysis

In the rest of the paper, we provide an empirical analysis of the designed partition crossover, denoted by Px.

Experimental Setup

We consider the following 8 QAP instances from the literature4 [14,15]: tai343e0i with i ∈ {0, . . . , 7}. They have been selected due to their challenging size of n = 343 facilities and locations, which is rarely addressed in the literature. We consider the following scenarios. Scenario #1. The goal of this scenario is the study the relative ability of the Px crossover to find an improving offspring. As depicted in the high-level template of Experiment 1, we consider two settings where the initial parents are either (i) random solutions, or (ii) local optima. For the latter case, we run a basic hill climbing local search using the standard swap neighborhood to construct the initial local optima. Starting from a random permutation, the best swap move is performed until no improvement is possible. The local search is executed as much times as needed to find as much different local optima as needed. All competing crossovers are applied 100 times using 100 pairs of different initial parents. Once an offspring has been generated by crossover, we also consider to check if it is a local optimum w.r.t. the swap neighborhood by running the local search again. Scenario #2. The goal of this scenario is to study the relative performance of the Px crossover when plugged into a simple evolutionary algorithm. As depicted in the high level template of Experiment 2, we consider a hybrid evolutionary algorithm embedding the swap-based local search discussed previously. In each generation, a new offspring is generated by performing crossover followed by a local search with probability p. We consider a simple random parent selection and a non-elitist replacement where the newly generated offspring replaces the oldest individual. The value of p is chosen in the set {0, 0.05, 0.2}. This allows us to study the relative behavior of crossover using a variable amount of local search, ranging from no local search at all (p = 0), to a small (p = 0.05) and a high (p = 0.2) amount. For each configuration, 10 independent runs are performed with a maximum number of generations G = 1 000 and a population size of n.

CPU running time.

Let us finally notice that we manage to analyze the CPU execution time. It is hence important to recall that the implementations of all algorithms were optimized as much as possible. In particular, computing the best move in a swap-based local search can be performed in a very efficient manner for the QAP. In fact, this can be done in an incremental manner on the basis of the current solution [START_REF] Taillard | Robust taboo search for the quadratic assignment problem[END_REF][START_REF] Taillard | Comparison of iterative searches for the quadratic assignment problem[END_REF]. For the sake of fairness when analyzing execution time, this well-established and important consideration is carefully implemented (a) Number of improving offspring q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 5 10 in all experiments. Moreover, when running the Px crossover, we restrict the maximum number of explored offspring to at most 2 15 , which means that when the number of connected components k > 15, not all possible offspring from Proposition 1 are explored. All algorithms are implemented in C++ and run on an Intel Xeon(R) CPU E3-1505M v6 3.00GHz. Let us finally notice that more advanced settings including other QAP instances and algorithms, as well as finely tuned operators/components using automated algorithm configuration methods, etc, are left for future investigations.

Experimental Analysis and Results

We start our analysis by reporting our findings from Experiment 1 in Fig. 3. Crossover improvement ratio. In Fig. 3a, we can see that the Px crossover has a significantly higher improvement rate. Using random parents, it can produce an improving offspring in almost 100% of the cases, whereas all other considered crossovers have a ratio of around 28%. Using local optimal parents, the Px is still able to find improving offspring in around 30% of the cases. This is to contrast with the other crossovers that fail in almost all cases. This first set of observations shows that the Px crossover has a relatively high intensification ability since it is even able to improve over local optima. Local optimality. In Fig. 3b, we further show the number of moves performed by the local search initialized with the generated offspring (in line 3), i.e., 0 moves means that the offspring is a local optima w.r.t. the swap neighborhood. For all crossovers except Px, the local search improves the constructed offspring independently of using random or local optima parents. This means that these crossovers are more diversification-oriented, since even when using local optimal parents, they are likely to produce an inferior offspring that can be improved by local search. Hence, we can think about these crossovers as acting in a perturbative manner. The situation is completely different for the Px crossover. On one side, with random parents, where the offspring improvement ratio was found to be almost 100%, the local search can still find improvements. On the other side, with local optimal parents, where the offspring improvement ratio is about 30%, Table 2: Ranks of crossovers (a lower rank is better). The first (resp. second) part is with respect to solution quality (resp. execution time). For each value of parameter p (local search ratio), a rank c indicates that the corresponding crossover was found to be significantly outperformed by c other ones w.r. the local search cannot find improvements, which indicates that the produced offspring is also likely to be a local optimum. This means that the Px crossover has some ability to act as a tunneling operator allowing to jump from two local optima to a new improving local optimum (with around 30% success rate).

Solution quality vs execution time. In the following, we report our findings from Experiment 2. In Table 2, we rank the different algorithms according to: (i) the quality of the best solution that the evolutionary algorithm is able to find during its execution, and (ii) the total CPU execution time. The situation is clearly different depending on whether a local search is used or not. When performing only crossover with no local search (p = 0), the Px crossover is able to find substantially better solutions, followed by Spx, and then by Opx, Ux and Cx, which do not show any significant difference statistically. Without surprise, these crossovers run however faster than Px and Spx. Notice here the extremely high cost of Spx. Interestingly, the proposed Px crossover implies a relatively reasonable increase in terms of running time (about 3.5 seconds) compared against Spx (about 287 seconds).

When performing crossover with local search (p ∈ {0.05, 0.2}), we can first see that, compared to not using local search at all, all variants can find much better solutions, while having a higher execution time. Interestingly, the Px and Spx crossovers (which were previously found to perform better than Opx, Ux, and Cx) now provide slightly worst solutions. On the other side, when analyzing the execution time, we found that using the Px crossover is significantly faster than all other variants, with at least 25% CPU time gain; e.g., for p = 0.2, the average execution time over all configurations is of 506 seconds with Px, against 677, 706, 714, and 871 with Px, Cx, Opx, Ux, and Spx, respectively. As commented before, knowing that our swap-based hill climbing local search was carefully implemented using a state-of-the-art fast incremental evaluation procedure for finding the best move [START_REF] Taillard | Robust taboo search for the quadratic assignment problem[END_REF][START_REF] Taillard | Comparison of iterative searches for the quadratic assignment problem[END_REF], such an observation might be surprising at first sight. However, it can be explained from two perspectives. Firstly, due to Corollary 1, performing the Px crossover is reasonably fast. Secondly, the Px crossover was previously found to have a relatively high intensification power. Thus, it is likely to produce a high-quality offspring, eventually being a local optima. Hence, it is more likely that the local search stops more quickly when attempting to improve the produced offspring. In this case, the cost of the local search is also reduced, hence leading to a decrease in the overall CPU time.

Finally, the previous results about solution quality and execution time are found to hold at any generation, independently of the configuration. This is illustrated in Fig. 4 rendering the convergence profile and the CPU execution time as a function of generations for the first QAP instance. This also confirms that the Px crossover is more intensification-oriented, and should be complemented by other diversification mechanisms when effectively integrated into more advanced evolutionary search processes.

QAP connected components. To complement our analysis, we provide further observations on the characteristics of the QAP and the relevance of the Px q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 1 2 3 4 5 crossover. More precisely, remember that given k connected components implied by the unshared facility assignments I, the Px crossover is able to provide the best over 2 k possible offspring. In Fig. 5, we report the average value of k (Left), as well the percentage of unshared facility assignments (Right), i.e., I /n in %, over all generations and all considered Px runs. We can see that for all considered QAP instances, solutions contain very few shared facility assignments, i.e., less than 1.5%, and the value of k stays relatively low, i.e., 6 in average with the exception of a few outliers exceeding 10. This suggests that the exploration power of the Px crossover can be improved in different ways, since its time complexity guarantees (Corollary1) should be able to support the fast evaluation of much more offspring solutions. For instance, it would be interesting to investigate the splitting of existing connected components into smaller ones at the aim of processing, and hopefully finding, more improving offspring.

Conclusion

In this paper, we presented our first investigations on the design of a partition crossover for the QAP. The proposed recombination and evaluation process is proved to provide a reasonable trade-off between the intensification power and the running time complexity. Our empirical study provides first insights towards the design and integration of more powerful partition crossovers for the QAP.

In particular, a future challenging issue is to investigate complementary decomposition techniques allowing to break the QAP bipartite graph into further smaller connected components, while maintaining a reasonable recombination and evaluation cost. In this respect, a promising idea would be to investigate the knowledge about the flow and distance values in order to identify the critical facilities and locations when performing decomposition. The challenge is then to guide the recombination process by identifying the most promising components to consider while keeping evaluation cost as low as possible.

  (a) The Cx crossover. (b) First steps of the Spx crossover.

Fig. 1 :

 1 Fig. 1: Illustration of the Cx and Spx crossovers.

Fig. 2 :

 2 Fig. 2: Illustration of a bipartite graph representation. Solid (resp. dashed) edges are with respect to G π1 (resp. G π2 ). We have I = {2, 8}, J = {3, 8} and the corresponding edges in gray. There are k = 3 connected components C 1 , C 2 and C 3 in G π1π2 [I, J], with U (C 1 ) = {1, 3, 4}, U (C 2 ) = {5, 6}, and U (C 3 ) = {7, 9}.

Corollary 1 .

 1 The whole set of offspring that can be generated by Proposition 1 can be explored and evaluated in O(n 2 + k 2 • 2 k ) time.Proof. Let us first notice that computing the components of the bipartite graph can be done in Θ(n) time. Over all the possible mapping functions m, (m( ), m( )) can only take 4 different values, namely, (1, 1), (1, 2), (2, 1) and (2, 2). Hence, for every fixed values of and , Q can only take 4 possible values. These four values depend solely on π 1 and π 2 . Therefore, there can only be Θ(k 2 ) possible values for Q over all and , and all possible choices of m. All of these Θ(k 2 ) values can be precomputed by a simple preprocessing step. Since by definition C and C do not share any nodes for every = , this preprocessing step takes obviously O(n 2 ) time. Similarly, there are only Θ(k) possible values for Q and Q over all possible values of , and all possible choices of m. Hence, they can also be precomputed in Θ(n 2 ) time. Finally, R does not depend on m, and can also be precomputed in O(n 2 ) time. To summarize, all possible values taken by Q , Q , Q , and R can be precomputed in Θ(n 2 ) time and stored in Θ(k 2 ) memory space before even any specific choice of the mapping function m is made. Now, let us consider a specific choice for the mapping function m, which fully determines an offspring permutation ⊗ π according to Proposition 1. Then, according to Proposition 2, and given the contributions Q , Q , Q , and R were already precomputed by the previous discussion, it takes Θ(k 2 ) time to compute f ( ⊗ π). The corollary follows since there are 2 k possible mapping functions m.

Experiment 1 :23Experiment 2 : 6

 126 Pseudo-code of the first experimental scenario 1 Let π1 and π2 be either (i) two random solutions or (ii) two local optima; Apply crossover on π1 and π2 to obtain an offspring ⊗ π; Apply a local search with ⊗ π as initial solution to obtain ⊗ π ; Pseudo-code of the second experimental scenario Input: G: maximum number of generations; p ∈ [0, 1]: local search ratio; 1 Generate an initial random population P ; 2 for g = 1 to G do 3 Apply crossover on two randomly selected parents π1 and π2 to obtain Replace the oldest individual of the population with ⊗ π;

Fig. 3 :

 3 Fig. 3: Results from Experiment 1. Results are over all eight instances.

Fig. 4 :

 4 Fig. 4: Experiment 2 on first instance. Top: convergence profile. Bottom: execution time. Left: p = 0 (no local search). Right: p = 0.2.

Fig. 5 :

 5 Fig. 5: Left: Number of connected components (k). Right: I /n in %.

  t. a Wilcoxon statistical test at a significance level of 0.05. Rows are w.r.t. instances.

		p = 0	p = 0.05	p = 0.2
	#Ins Opx Ux Cx Spx Px Opx Ux Cx Spx Px Opx Ux Cx Spx Px
		Mean average deviation to the best (in subscript)
	1	2810 2809 2808 1796 0710 00.38 00.39 00.37 00.37 00.37 00.26 00.24 00.30 10.32 10.34
	2	2105 2105 2105 1103 094 00.04 00.04 00.03 00.04 10.05 00.02 00.03 00.02 00.03 00.03
	3	2115 2115 2115 1113 0102 00.04 00.03 00.03 00.04 10.05 00.02 00.02 00.02 30.03 30.03
	4	296 296 296 195 084 00.04 00.04 00.04 00.04 10.05 00.03 00.02 00.02 10.03 20.03
	5	2107 2107 2108 1106 096 00.04 00.03 00.03 00.04 10.04 00.02 00.02 00.02 10.03 20.03
	6	2123 2123 2124 1121 0110 00.06 00.05 00.05 00.06 10.06 00.04 00.04 00.03 10.05 10.05
	7	2115 2115 2115 1113 0101 00.04 00.04 00.04 00.05 00.04 00.03 00.03 00.04 00.04 00.03
	8	2110 2110 2110 1108 0 96 00.04 00.03 00.03 00.03 10.04 00.02 00.02 30.03 30.03
			Mean CPU execution time (in subscript)
	1	00.61 10.62 20.75 4286 33.39 0194 0198 0196 4478 0185 1705 1714 1676 4870 0502
	2	00.61 10.62 20.72 4287 33.44 0202 0202 0199 4475 0187 2698 1697 1662 4857 0494
	3	00.60 10.63 20.70 4286 33.38 0210 0208 0206 4481 0195 1712 1730 1682 4876 0510
	4	00.60 10.62 20.75 4287 33.41 0193 1198 0194 4468 0175 1724 1728 1693 4880 0513
	5	00.60 10.62 20.74 4287 33.39 0189 0187 0183 4463 0172 1681 2684 1643 4842 0488
	6	00.61 10.62 20.72 4286 33.39 0209 0212 0206 4481 0193 1714 1725 1690 4876 0515
	7	00.61 00.62 20.69 4286 33.37 0207 0206 0205 4477 0193 1697 2701 1670 4875 0505
	8	00.61 10.62 20.70 4287 33.40 0221 0223 0219 4494 0200 1715 1734 1693 4890 0521

http://mistic.heig-vd.ch/taillard/problemes.dir/qap.dir/qap.html

Acknowledgments. The three first authors are supported by the French national research agency (ANR-16-CE23-0013-01) and the Research Grants Council of Hong Kong (RGC Project No. A-CityU101/16).


Martin Drozdik
email: martin@iplab.shinshu-u.ac.jp

Kiyoshi Tanaka

Hernan Aguirre

Sébastien Verel

Arnaud Liefooghe

Bilel Derbel

Sebastien Verel

An Analysis of Differential Evolution Parameters on Rotated Bi-objective Optimization Functions

Keywords: differential evolution, rotational invariance, multi-objective optimization, parameter analysis

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

Introduction

Differential evolution [START_REF] Price | Differential Evolution: A Practical Approach to Global Optimization[END_REF] started as a simple single-objective continuous optimization heuristic. The need for a versatile multi-objective optimizer has motivated researchers to generalize the basic algorithm for multi-objective problems. Now we have a great number of multi-objective DE variants. Many of them use the same mechanism to generate new individuals. In a problem with n variables a new individual is created using a crossover variation operator which randomly selects k; k ≤ n variables which are perturbed. The magnitude of the mutation is generated by scaling a difference of randomly chosen individuals.

Many research papers on DE such as [START_REF] Denysiuk | Many-objective Optimization Using Differential Evolution with Variable-wise Mutation Restriction[END_REF] or [START_REF] Tea | DEMO: Differential Evolution for Multiobjective Optimization[END_REF] provide little insight into how the authors chose the parameters for their benchmarking. We find this striking since many authors choose their parameters such that the crossover operator perturbs only a small number of variables in an existing individual. In other words, the search for the Pareto optimal set proceeds along the coordinate axes. Since these algorithms perform very well [START_REF] Denysiuk | Many-objective Optimization Using Differential Evolution with Variable-wise Mutation Restriction[END_REF][9], we have a suspicion that this may be due to some characteristic of the problem, such as separability, that makes it easy to optimize along the axes. This would mean that if the axes are transformed, the algorithm should lose some performance.

Very strict warning against the practice of perturbing a small number of variables at a time has been raised as soon as 1996 by Salomon [START_REF] Salomon | Re-evaluating Genetic Algorithm Performance Under Coordinate Rotation of Benchmark Functions[END_REF]. Salomon empirically demonstrated that the stellar performance of many popular singleobjective genetic algorithms owes to the fact that most of the benchmark functions were separable and that the low mutation rate caused them to be optimized one component at a time. Once Salomon stripped the separability by rotating the principal axes of the benchmark functions, many algorithms were significantly slowed down, while some failed to converge completely. Salomon's theoretical results state that, in some cases, the probability of finding the global optimum can drop below that of random search. We are concerned that the same is true for the multi-objective realm since many authors perform their experiments with separable test functions.

In DE the number of variables that are perturbed is controlled by a parameter. If all variables are perturbed, the algorithm has the same performance regardless of rotation. Let us have a parameter setting, that perturbs only a small proportion of the variables, which outperforms a setting that perturbs all variables. In this work we attempt to answer this question: Is this exceptionally good performance on a problem with a particular alignment of the coordinate axes balanced by exceptionally bad performance on a different alignment?

We do this empirically by observing the performance of a simple multiobjective algorithm DEMO (Differential evolution for multi-objective optimization) [START_REF] Tea | DEMO: Differential Evolution for Multiobjective Optimization[END_REF] on a bi-objective subset of the WFG (Walking Fish Group) test suite [START_REF] Huband | A Review of Multiobjective Test Problems and a Scalable Test Problem Toolkit[END_REF]. We run all our experiments with a fixed population size and a fixed number of variables, while varying the parameters. Then we gradually rotate the problems in a controlled manner and observe the new behavior.

The answer to our question is, unexpectedly, negative. We find a statistically significant difference between the performance on the rotated problems and the original ones. Closer inspection reveals that a systematic performance loss happens when we rotate the separable problems, but the performance is still significantly better than for a rotationally invariant algorithm. We find that this happens for multi-modal problems, while single-modal problems exhibit the behavior we would expect from the work of Salomon.

In the following section we provide background information on DE and on the previous related work on DE parameters. In Section 3 we introduce the experimental design, where we explain which problems are used and why were they chosen. In addition we introduce a new performance metric called the relative hypervolume, and explain the controlled manner in which the rotations are generated. In Section 4 we present our data along with a discussion. Finally, in Section 5, we present the conclusion. Trim the P to size N using non-dominated sorting [START_REF] Tea | DEMO: Differential Evolution for Multiobjective Optimization[END_REF] and MNN diversity [START_REF] Kukkonen | A Fast and Effective Method for Pruning of Nondominated Solutions in Many-Objective Problems[END_REF] 24 end 2 Background

Differential Evolution

In this section, we describe the variant of DE which we use in this work. It is a slightly modified version of the DEMO algorithm [START_REF] Tea | DEMO: Differential Evolution for Multiobjective Optimization[END_REF] described in Algorithm 1. The modified parts are highlighted with yellow color in lines 14 and 23. Let us look at the algorithm in detail. First, the population of the algorithm is randomly initialized (line 1). Then, the algorithm runs for a fixed number of generations (evolutionary loop). In each generation, DE iterates through the entire population generating a trial individual, which is compared to an existing target individual. The trial is generated utilizing the traditional method in lines 4 to 13. Here we introduce the parameters of DE.

To generate a new individual, three distinct individuals are selected from the population. By forming a difference between two of them, scaling it by a fixed parameter F, and adding to a third individual we obtain a so called mutant individual (line 5). The trial individual is created by crossover between the mutant and the target. First, a randomly chosen variable from the mutant is inherited (line 6). Next, each other variable is inherited from the mutant with a fixed crossover probability of Cr. Otherwise it comes from the target.

After the individual is generated we project it to the decision space. This is a modification of the original algorithm which did not explicitly deal with domain issues. The purpose is to keep the algorithm as simple as possible while being able to optimize problems with simple constraints. Next the trial is compared to the target. If one of them is dominated by the other, we discard the dominated one. If they are mutually non-dominated, we keep them both. At the end of the generation loop, the population is trimmed to size N using non-dominated sorting and the M nearest neighbor diversity estimation procedure [START_REF] Kukkonen | A Fast and Effective Method for Pruning of Nondominated Solutions in Many-Objective Problems[END_REF]. We chose this procedure because it achieves a better distribution along the Pareto front than the original crowding distance computation.

Note that Cr = 1 is the only value of Cr for which the DE algorithm is rotationally invariant with probability 1. Rotational invariance does not by itself imply good performance. Its merit is that it allows us to generalize a single observation to an entire invariance class [START_REF] Hansen | Impacts of Invariance in Search: When CMA-ES and PSO Face Ill-Conditioned and Non-Separable Problems[END_REF].

Crossover Probability and Separability

In this section we summarize what we know about the relationship between the separability of the test functions and the good choice of Cr parameter.

There are many different types of separability. One of the simplest is additive

separability. A function f : D ⊆ R n → R is called additively separable if: ∃f 1 , . . . , f n such that f (x 1 , . . . , x n ) = n i=1 f i (x i )
The most important consequence of additive separability is that the n-dimensional problem can be optimized sequentially one variable at a time. Therefore separable problems are not subject to the curse of dimensionality [START_REF] Salomon | Re-evaluating Genetic Algorithm Performance Under Coordinate Rotation of Benchmark Functions[END_REF].

Salomon [START_REF] Salomon | Re-evaluating Genetic Algorithm Performance Under Coordinate Rotation of Benchmark Functions[END_REF] illustrates the problems of algorithms which vary the individuals one variable at a time on a quadratic function of two variables in Figure 1. The ellipses in the left part are contours of a separable quadratic function. We can see two individuals on one of the contours. The blue individual represents an individual in a randomized algorithm. If we mutate one variable of this individual, the probability to get an improvement in the objective function is relatively high, since the improvement intervals d 1 , d 2 are long. If we rotate the coordinate axes, thus rendering the function non-separable, the improvement intervals shrink.

One more illustration of problems which arise is using a sequential deterministic algorithm which finds the optimum with respect to one variable at a time. The red individual illustrates the path of one such an algorithm. When the function is aligned with the axes, this algorithm achieves optimum in just two iterations, while in the rotated case the algorithm not only progresses slower, but never actually reaches the optimum.

Huband et al. from the Walking Fish Group (WFG) define separability from the optimizational standpoint [START_REF] Huband | A Review of Multiobjective Test Problems and a Scalable Test Problem Toolkit[END_REF] The multi-objective model is fundamentally different from the single-objective model because all objectives are being optimized simultaneously. The global optima of each optimized function constitute only a relatively small subset of the Pareto optimal set. Therefore, it is appropriate to ask if the problems of sequential algorithms which are illustrated in Figure 1 persist in multi-objective optimization. Also, while additively separable unimodal functions are inherently similar to the quadratic function in Figure 1, it is not clear if the intuition holds for multi-modal functions or for functions which are WFG-separable but not additively separable.

Variance as a Common Currency

Probably the most significant work on the theoretical properties of DE has been written by Zaharie [START_REF] Zaharie | Critical Values for the Control Parameters of Differential Evolution Algorithm[END_REF]. Let us collect all the trial vectors that are generated in the course of one generational loop of Algorithm 1 into a set P trial . Then the relationship between the variance in decision space of P and P trial is given by the simple equation E[V ar(P trial )] = cE[V ar(P )] where:

c = 2F 2 Cr + Cr 2 -2Cr N + 1 (1)
Zaharie omits the fact that in most DE variants the individuals which generate the trial individual are chosen distinct from the target individual (Algorithm 1 line 4). However her results hold unchanged also after adding this assumption.

The work of Zaharie is important since it transforms the two parameters into a single number c (common currency) which has a very intuitive interpretation. If c < 1 we see that the algorithm tends to contract the population while if c > 1 it expands the population. Based on empirical data Kukkonen concluded in [START_REF] Kukkonen | An Empirical Study of Control Parameters for The Third Version of Generalized Differential Evolution (GDE3)[END_REF] [7] that a good choice of parameters is one that satisfies c ∈ [1.0; 1.5] with the upper bound not very strict.

Experimental Design

In this section, we describe which test problems we chose and why. We explain what we mean by rotating the problem and we propose a new performance metric which we use.

WFG Problems

In order to explore the relationship between the control parameters of DE and the characteristics of the problem, we chose 4 problems from the WFG test suite [START_REF] Huband | A Review of Multiobjective Test Problems and a Scalable Test Problem Toolkit[END_REF]. These problems have been chosen since they have the same Pareto front and contain all possible combinations of the WFG-separability and modality characteristics. They are summarized in Table 1. We chose the number of variables to be 10 of which one is a positional variable.

Rotations in R n

As humans we have a very good intuitive understanding of rotation in 2 or 3 dimensional space. However in higher dimensions things are not as intuitive as they might seem. An elementary rotation by the angle φ is characterized by the matrix:

R e = cos(φ) -sin(φ) sin(φ) cos(φ)

We can generalize this rotation to n-dimensional space by taking an n-dimensional identity matrix I and replacing I i,i , I i,j , I j,i , I j,j by R e 1,1 , R e 1,2 , R e 2,1 , R e 2,2 respectively. We can see that the rotation is not executed around an axis as we might intuitively feel, but around an n-2 dimensional subspace which is coincidentally a 1-dimensional axis in the intuitive 3-dimensional case. For our experiments, we generate the rotation matrix R by applying a rotation to each n -2 dimensional subspace in sequence, one rotation after the other.

We rotate the entire decision space (DS). This way the entire Pareto optimal set is always attainable since the entire decision space rotates along. In the case of WFG problems this means rotating a n-dimensional hyper-box. For example, in order to initialize the population in Algorithm 1 in the rotated DS (line 1), we first initialize the population in the original DS and then multiply by R -1 . Similar process is used to project the individual to the rotated DS on line 14. To evaluate the objective value of an individual we first multiply the decision vector by R and evaluate the original objective functions.

Relative Hypervolume

In our experiments, we use only one performance metric, the hypervolume (HV) [START_REF] Zitzler | Evolutionary Algorithms for Multiobjective Optimization: Methods and Applications[END_REF], since it includes information on both convergence and spread of the individuals. With WFG problems, it is not easy to choose the reference point for the HV. Even if we choose the point as tight as possible, there are some individuals after the initialization of DE which dominate the reference point. Therefore the HV at the start is not zero and it is hard to say if a certain attained HV is good or bad. Moreover, it is hard to make quantitative comparisons based on HV. If some algorithm achieves HV of 100 and another one achieves a HV of 99.98, it may seem that the difference is not very big, but it all depends on the HV at initialization. If the algorithms started with HV = 0, the interpretation of the results would be quite different from one where HV = 99.99 at the start.

We attempt to mitigate this problem by subtracting the HV at initialization (HV init ) and normalizing the result using the maximal attainable hypervolume (HV max ). We define the relative hypervolume (RHV) in the following equation:

RHV := HV -HV init HV max -HV init (2)
We compute HV max deterministically by integrating the space between the true Pareto front (PF) and the reference point. From (2), we have RHV ∈ [1; -∞). RHV = 1 implies convergence, RHV at initialization is 0 and RHV < 0 indicates an algorithm which is receding from the Pareto front. We use RHV since its normalized nature is more intuitive and it is more robust with respect to the selection of the reference point. It may be more meaningful to compare two algorithm runs in terms of RHV. If we have two algorithm runs starting from the same randomly initialized population then the ratio of their RHVs is independent of the choice of the reference point. 4 On the other hand, two independent runs which produce the same final population may yield different relative hypervolume.

Results and Discussion

In our experiments we varied the parameters F ∈ [0.05; 1.5], Cr ∈ [0; 1] equidistantly with a resolution of 0.05. For each combination we performed 10 runs of Algorithm 1. To simplify the setup, the population size was kept constant at 100 individuals and the length of each run was 250 generations. We explored the rotations from 0 to 90 degrees with a resolution of 5 degrees. In the following we discuss our results on a subset of the experimental data. To simplify the analysis, in each section we keep either F, Cr or the rotation angle fixed.

Fixed Rotation Angle

Figure 2 shows the average RHV on non-rotated problems. For illustration, we plot the combinations of F and Cr which result in c = 1.0, 1.5 and 3.0 according to [START_REF] Denysiuk | Many-objective Optimization Using Differential Evolution with Variable-wise Mutation Restriction[END_REF]. The circle marks the combination of parameters with the best RHV. The data in all our figures from now on is presented from left to right in the same order as in Table 1. The separable problems (S) are on the left, the non-separable (NS) on the right, while the unimodal (UM) are on the inside and multi-modal (MM) ones are near the page margins. For each problem, an L-shaped favorable region containing RHV of 0.8 and higher, roughly corresponds to c ∈ [1; 1.5]. Low value of Cr is more robust, since it allows for a wider interval of F values. Unexpectedly, this holds also for non-separable problems WFG6 and WFG9.

The effect of introducing a rotation by 5 degrees is shown in Figure 3. The two figures seem identical, but the ratio of these averages in Figure 4 reveals a difference. A value of less than 1 indicates that the rotation caused the performance to decrease. We highlighted the contour at level 1 and marked the maximal and minimal value by circles. In order to make the results most readable we chose a color scale of [0.5; 1.2] for separable problems and [0.6; 1.7] for non-separable problems. The separable problems on the left half exhibit a performance loss consistent with Salomon's single-objective results. Performance dropped for almost all Cr smaller than 1. Non-separable WFG6 and WFG9 do not show such a systematic decrease. In some areas we even see an increase of performance. It seems that there is relatively little difference between the rotated and nonrotated data. These result may seem not as significant as Salomon's. However, there is an important methodological difference. When he mentions that the performance on the rotated benchmark is six orders of magnitude worse than the performance on the non-rotated benchmark ([8, p.273]), he means that the minimal attained value 2.65 • 105 is six orders of magnitude worse in absolute numbers. But the value at initialization was three orders of magnitude greater yet. This means that both algorithms started somewhere near 2.65 • 10 8 and the non-rotated one progressed to 2.65 • 10 -1 while the rotated one progressed to 2.65•10 5 . In terms of relative hypervolume, this would be a very small difference 5 . In order to provide a scale-independent comparison, we compared all data using a two-tailed Wilcoxon signed rank test at a significance level of 0.05. For separable problems in Figures 4 and5 we separate the parameter space with a dashed line into two areas. The area on the right is such that the rotated and non-rotated data is not significantly different, while on the left there is a significant decrease in performance. The data for non-separable problems contains areas of both significant decrease and significant decrease, as well as areas with no significant difference so in this case the separation cannot be plotted so compendiously. The effects are more visible with 45 degree rotation in Figure 5. Again, there is a systematic decrease in performance for the separable problems for Cr < 1. However, this decrease does not imply that Cr = 1 is a good choice. Looking at Figures 2 and3, we see that Cr = 1 is a consistently bad choice for the multi-modal problems WFG4 and WFG9.

Fixed F

In Figures 2 and3 we see that F = 0.5 is compatible with many different values of Cr and achieves consistently good performance. The average RHV for F = 0.5 is shown in Figure 6. For multi-modal problems WFG4 and WFG9, very low values of Cr are consistently good for all studied rotations, while for uni-modal problems WFG6 and WFG7 big values of Cr yield a consistently good performance. On the other hand, poor performance is achieved with big values of Cr for multi-modal problems and small values for uni-modal problems. The data for WFG4 and WFG9 suggests that the exceptionally good performance of a small Cr setting does not have to be balanced by an exceptionally bad performance after the problem is rotated. Based on the observation from Figure 6 we see that for each problem either Cr = 0.1 or Cr = 0.9 perform well through the observed spectrum of rotations.

Fixed Cr

Figures 7 and8 we see data with a fixed value of Cr = 0.1 and Cr = 0.9 respectively. For Cr = 0.1 the regions with the best performance are for rotations which are either close to 0 or 90 degrees. This is true also for non-separable problems, but it is more visible for separable problems. The data for Cr = 0.9 seems different. The choice of Cr close to 1 means that the algorithm is nearly rotationally invariant. The gained robustness with respect to coordinate rotation is balanced by lost robustness in the choice of F. Almost in all cases the interval with favorable values of F became shorter.

We see that for values of Cr < 1 there is a performance loss when the coordinate axes are rotated, but does the performance drop bellow that of a rotationally invariant choice of Cr = 1? The data supporting a negative answer is presented in Figure 9. Here we divided the average RHV with Cr = 0.1 by the average RHV attained with a rotationally invariant Cr = 1. The interpretation of the dashed and full contour lines is the same as for Figures 4 and5. For WFG4, the setting of Cr = 0.1 statistically significantly outperformed Cr = 1 for all rotations and all values of F. This means a definitive negative answer to our main question. The results are similar for the second multi-modal problem WFG9. Here we see a small region in which the data for Cr = 0.1 and Cr = 1 are not significantly different and Cr = 1 is significantly better in a few isolated cases. The unimodal problems on the other hand show that Cr = 1 is significantly better for most rotations and for the best performing values of F.

Conclusion

In this work we showed how the behavior of the differential evolution algorithm on bi-objective problems changes when the coordinate axes of the decision space are rotated. Our findings show that the change is significant even for small rotations. There is a consistent drop in performance on separable problems while the qualitative properties of the change for non-separable problems are much less predictable. Unexpectedly, for multi-modal problems, low values of crossover probability perform better through the observed spectrum of rotations. As a future work we propose to see if this holds for problems other than the ones we studied and if this is the case, to find the cause of this behavior.

Algorithm 1 : 9 X 11 X

 1911 Modified DEMO [9] algorithm 1 initialize P = {X1, ..., XN } uniformly randomly in the decision space 2 for generation := 1 to Gmax do Evolutionary loop 3 for target := 1 to N do Generational loop 4 randomly generate mutually distinct r1, r2, r3 = target 5 Xmutant := Xr 1 + F(Xr 2 -Xr 3 ) 6 randomly generate inv ∈ {1, . . . , n} 7 for i := 1 to n do 8 if rand(0.0; 1.0) < Cr or i = inv then trial,i := Xmutant,i 10 else trial,i := Xtarget,i 12 end 13 end 14 project X trial to decision space 15 if Xtarget dominates X trial then 16 discard X trial 17 else if X trial dominates Xtarget then 18 replace Xtarget with X trial 19 else if Xtarget and X trial are mutually non-dominated then 20 add X trial to the end of the population 21 end 22 end 23

Fig. 2 :

 2 Fig. 2: Average RHV without rotation

Fig. 3 :

 3 Fig. 3: Average RHV with rotation angle of 5 degrees

Fig. 4 :Fig. 5 :

 45 Fig. 4: Average RHV with a rotation of 5 degrees Average RHV without a rotation

Fig. 6 :Fig. 7 :

 67 Fig. 6: Average RHV for F = 0.5

Fig. 8 :Fig. 9 :

 89 Fig.8: Average RHV for Cr = 0.9

  . A variable x i is separable if the set of global optima of a problem: argmin Fig. 1: Illustration of variable-wise optimization on a rotated quadratic function is the same for any choice of the other variables x 1 , . . . , x i-1 , x i+1 , . . . , x n . For example, an additively separable function is WFG-separable, hence WFG-separability is a generalization of additive separability. The authors define a separable multiobjective problem as one where each objective is separable. The majority of the frequently used DTLZ and ZDT problems are WFG-separable [3], while their objective functions are not additively separable.

	x 2	Improvement intervals Separable (aligned) quadratic function	x 2	Non-separable (rotated) quadratic function
		Global optimum
		Sequential optimizer
	d 2	x 1	d 2	x 1
		d 1		d 1

xi f (x 1 , . . . , x n )

Table 1 :

 1 Characteristics of the selected WFG problems

	WFG4 WFG7 WFG6 WFG9
	separable yes	yes	no	no
	unimodal no	yes	yes	no

Given that the reference point is dominated by all individuals in the population.

Assuming that the minimum of the given function is 0, the difference would be on the order of 10 -3 .