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

Roberto Santana
email: roberto.santana@ehu.eus

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

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

Multi-objective NK landscapes with heterogeneous objectives

Keywords: Multi-objective optimization, heterogeneous objectives, landscape analysis, NK landscapes

published or not. The documents may come

INTRODUCTION

Multi-objective optimization problems (MOPs) with heterogeneous objectives are those where the objective functions differ in one or several aspects, such as scaling, landscape, evaluation time, or theoretical and practical difficulty [START_REF] Greco | Understanding complexity in multiobjective optimization (Dagstuhl Seminar 15031)[END_REF]. Studies on multi-objective evolutionary algorithms (MOEAs) do not usually pay attention to the heterogeneity of the objectives, and to the way it might influence search difficulty. This is particularly the case for multiobjective benchmark functions where the focus is usually put on global characteristics of the multi-objective landscape, or on the shape of the Pareto front (PF). However, real-world problems might exhibit a significant variability in the objectives' characteristics, and recent research has addressed the question of heterogeneous objectives and proposed multi-objective approaches to deal with them [START_REF] Allmendinger | Heterogeneous objectives: state-of-the-art and future research[END_REF]. Much of this previous research has been, however, focused on problems where the heterogeneity arises in evaluation times or latencies, that is, when each objective takes a different amount of time to be evaluated [START_REF] Allmendinger | Multiobjective optimization: When objectives exhibit non-uniform latencies[END_REF][START_REF] Allmendinger | hang on a minute': Investigations on the effects of delayed objective functions in multiobjective optimization[END_REF][START_REF] Blank | Constrained bi-objective surrogate-assisted optimization of problems with heterogeneous evaluation times: Expensive objectives and inexpensive constraints[END_REF][START_REF] Chugh | Surrogate-assisted evolutionary biobjective optimization for objectives with non-uniform latencies[END_REF].

In this paper, we investigate multi-objective problems with related, but heterogeneous landscapes. We are in particular interested in problems where the multi-modality and the ruggedness significantly differ among the objectives. As a representative class of problems, we select multi-objective NK landscapes [START_REF] Aguirre | Insights and properties of multiobjective MNKlandscapes[END_REF][START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF][START_REF] Knowles | Quantifying the effects of objective space dimension in evolutionary multiobjective optimization[END_REF][START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF], and we build on recent results analyzing and characterizing this problem class [START_REF] Cosson | Decomposition-based multi-objective landscape features and automated algorithm selection[END_REF][START_REF] Daolio | Problem features versus algorithm performance on rugged multiobjective combinatorial fitness landscapes[END_REF][START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multi-objective optimization[END_REF][START_REF] Liefooghe | On Pareto local optimal solutions networks[END_REF]. The parameter 𝐾 of the NK model critically influences the characteristics of the fitness landscape. As 𝐾 increases, the ruggedness and the multi-modality of the landscape also increase. Therefore, we investigate multi-objective NK landscapes where each objective has a different degree of variable interactions 𝐾. Such a model serves as a framework to investigate the impact of heterogeneous objectives, by allowing us to evaluate how the difference in the amount of ruggedness among the NK models translates into the difficulty of the MOP. For example, we can quantify how the increment in the number of single-objective local optima for one of the objectives affects the number of multiobjective local optima for MNK landscapes. Although our analysis could be extended to an arbitrary number of objectives, in this paper we focus on bi-objective NK landscapes. Furthermore, we present the analysis from the perspective of optimization strategies based on a Pareto dominant approach.

To evaluate the impact of heterogeneity, we propose to construct a rank-annotated neighborhood network with labeled local optimal solutions. This network is inspired by neighborhood graphs traditionally applied to fitness landscapes [START_REF] Sutton | An analysis of combinatorial search spaces for a class of NP-hard problems[END_REF], and by Pareto local optimal solutions networks [START_REF] Liefooghe | On Pareto local optimal solutions networks[END_REF] previously used to extract informative features from multi-objective NK landscapes. This annotated neighborhood network comprises both local optimal and non-local optimal solutions. Furthermore, local optimal solutions are labeled to indicate the different classes they belong to: Pareto local optimal solutions [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF], and (single-objective) local optima for each objective. In addition, the nodes from the network are layered according to the rank of their corresponding solution with respect to non-dominated sorting [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF], similarly to the strategy followed in Pareto plateau connection graphs [START_REF] Garrett | Plateau connection structure and multiobjective metaheuristic performance[END_REF].

Finally, we conduct exhaustive experiments on heterogeneous bi-objective NK landscapes with different degrees of variable interactions per objective, and different degrees of heterogeneity. Using traditional metrics from landscape analysis, and others extracted from the labeled rank-annotated neighborhood network, we analyze the heterogeneous landscapes in terms of multi-modality and of the Pareto set structure. We further analyze to what extent these metrics are impacted by the degree of variable interactions and by the degree of heterogeneity among the objectives.

The paper is organized as follows. Section 2 presents the background on multi-objective NK landscapes and reviews related work. Section 3 introduces the heterogeneous MNK model, proposes ways to quantify heterogeneity, and discusses its usability. Section 4 describes the rank-annotated neighborhood network with labeled local optimal solutions. Section 5 gives the experimental setup, and presents and discusses numerical results. Section 6 concludes the paper and mentions a number of open research lines.

BACKGROUND AND RELATED WORK 2.1 NK and MNK Landscapes

Let X = (𝑋 1 , . . . , 𝑋 𝑁 ) denote a vector of binary variables. We use x = (𝑥 1 , . . . , 𝑥 𝑁 ) to denote an assignment to the variables. 𝑆 denotes a set of indices in {1, . . . , 𝑁 }, and 𝑋 𝑆 (respectively 𝑥 𝑆 ) a subset of the variables of X (respectively x) determined by the indices in 𝑆.

The model for NK landscapes was originally introduced to study the effect of different patterns of variable interactions in the landscape [START_REF] Kauffman | Origins of Order[END_REF]. It allows to explore the way in which the structure and the strength of interactions among variables determine the ruggedness of the landscape. For given parameters, the problem consists of finding the global maximum of the function.

An NK landscape is defined with the following components:

• Number of variables, 𝑁 .

• Number of interactions per variable, 𝐾.

• A set of 𝐾 interactions Π(𝑋 𝑗 ) ∈ X, for 𝑋 𝑗 , 𝑗 ∈ {1, . . . , 𝑁 }.

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

The objective function 𝐹 , to be maximized, is defined as follows:

𝐹 (x) = 1 𝑁 𝑛 ∑ 𝑗=1 𝑓 𝑗 (𝑥 𝑗 , Π(𝑥 𝑗 )). (1)
Multi-objective NK (MNK) landscapes [START_REF] Aguirre | Insights and properties of multiobjective MNKlandscapes[END_REF][START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF] provide an extension of the NK model to the multi-objective case. An MNK landscape is defined as a vector function mapping binary vectors into 𝑀 real numbers F(.) = (𝐹 1 (.), 𝐹 2 (.), . . . , 𝐹 𝑀 (.)) : B 𝑁 ← R 𝑀 , where 𝑁 is the number of variables, 𝑀 is the number of objectives, 𝐹 𝑖 (.) is the 𝑖-th objective function, and B = {0, 1}. K = (𝐾 1 , . . . , 𝐾 𝑀 ) is a tuple of integers where 𝐾 𝑖 is the size of the neighborhood in the 𝑖-th landscape 1 . Each 𝐹 𝑖 (x) is defined similarly to Equation (1).

Definitions

Let 𝑍 denote the mapping of solutions to MNK landscapes in the objective space. Given two objective vectors 𝑧, 𝑧 ′ ∈ 𝑍 , 𝑧 is dominated by 𝑧 ′ iff for all 𝑖 ∈ {1, . . . , 𝑀 } 𝑧 ′ 𝑖 ≥ 𝑧 𝑖 , and there is a 𝑗 ∈ {1, . . . , 𝑀 } such that 𝑧 ′ 𝑗 > 𝑧 𝑗 . Similarly, given two solutions 𝑥, 𝑥 ′ ∈ 𝑋 , 𝑥 is dominated by 𝑥 ′ iff 𝐹 (𝑥) is dominated by 𝐹 (𝑥 ′ ). An objective vector 𝑧 ★ ∈ 𝑍 is non-dominated if there does not exist any 𝑧 ∈ 𝑍 such that 𝑧 ★ is dominated by 𝑧. A solution 𝑥 ★ ∈ 𝑋 is Pareto optimal, or non-dominated, if 𝐹 (𝑥) is non-dominated. The set of Pareto optimal solutions is the Pareto set (PS); its mapping in the objective space is the Pareto front (PF). One of the main challenges in multi-objective optimization is to identify the PS, or a good approximation of it for large-size and complex problems. A number of MOEAs have been designed to this end since the late eighties [START_REF] Coello | Evolutionary Algorithms for Solving Multi-objective Problems[END_REF][START_REF] Deb | Multi-objective Optimization Using Evolutionary Algorithms[END_REF].

Related Work

Among the research lines related to our work are approaches that study the influence of different factors in MNK landscapes, and works that propose metrics that serve to characterize MOPs, particularly metrics that result from a network-based landscape analysis. In this section, we cover some of the research in these directions.

Aguirre et al. [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF] present an exhaustive investigation of how the parameters of MNK landscapes influence several characteristics of the landscape, including the size of the Pareto front and the number of fronts. They extract a number of valuable conclusions related to the effects that the parameters of the MNK landscape have in the complexity of the problems. Since both the variable interactions' structure and the local functions of the NK model influence the ruggedness of the landscape, the impact of these factors have also been investigated from different perspectives. In [START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF][START_REF] Verel | Local optima networks of NK landscapes with neutrality[END_REF], the NK model has been used to create 𝜌MNK landscapes, where the same variable interactions' structure is used for each objective and the local sub-functions are modified to enforce desired correlations among the objectives. The model has been applied to investigate MOEAs in different scenarios [START_REF] Derbel | On the impact of multiobjective scalarizing functions[END_REF][START_REF] Marquet | Shake them all! rethinking selection and replacement in MOEA/D[END_REF]. Knowles and Corne [START_REF] Knowles | Quantifying the effects of objective space dimension in evolutionary multiobjective optimization[END_REF] use the MNK model to investigate the effect of augmenting the number of objectives considering different strength of correlation between them. Another way to extend the NK model to the multi-objective domain is by using a single and fixed set of local functions for each objective, and exploring the difficulty of the search by modifying the variable interactions, as proposed in [START_REF] Santana | Evolving MNK-landscapes with structural constraints[END_REF]. It is important to notice that all previous analyses of MNK landscapes have considered a fixed 𝐾 for all the involved objectives. In this paper, we are interested in MNK models with different values for 𝐾.

There has been an increasing number of works that propose features which inform about the characteristics or the difficulty of combinatorial MOPs [START_REF] Cosson | Decomposition-based multi-objective landscape features and automated algorithm selection[END_REF][START_REF] Daolio | Problem features versus algorithm performance on rugged multiobjective combinatorial fitness landscapes[END_REF][START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multi-objective optimization[END_REF][START_REF] Liefooghe | On Pareto local optimal solutions networks[END_REF]. In addition to providing tools for a better understanding of MOPs, these features have been applied to explain algorithm performance or implement criteria for algorithm selection [START_REF] Cosson | Decomposition-based multi-objective landscape features and automated algorithm selection[END_REF][START_REF] Liefooghe | Landscapeaware performance prediction for evolutionary multi-objective optimization[END_REF]. We show in this paper how several of these metrics can describe particular aspects of heterogeneous MOPs. We complement the information provided by traditional features with other metrics extracted from a network representation of the proposed heterogeneous problem.

In landscape analysis, the landscape is often defined by means of a neighborhood graph, where there is a vertex for each solution and edges connect neighbors. Local optima networks (LONs) [START_REF] Ochoa | A study of NK landscapes' basins and local optima networks[END_REF] serve as a network characterization of fitness landscapes. Each vertex corresponds to a local optimum, and edges represent basin adjacency between two local optima. They were originally applied to the NK model, and have been later extended to multi-objective NK landscapes in [START_REF] Liefooghe | On Pareto local optimal solutions networks[END_REF]. In the Pareto local optimal solutions network (PLOSnet), vertices correspond to Pareto local optima solutions [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF], and edges are constructed between mutually non-dominated neighbors. By definition, a PLOS-net is a sub-graph of the neighborhood graph. Network metrics have been proposed for both LONs and PLOS-nets [START_REF] Chicano | Local optima networks, landscape autocorrelation and heuristic search performance[END_REF][START_REF] Herrmann | Communities of local optima as funnels in fitness landscapes[END_REF].

There are other approaches that represent the multi-objective landscape as graphs. This is the case of Pareto plateau connection graphs [START_REF] Garrett | Plateau connection structure and multiobjective metaheuristic performance[END_REF], where the nodes from the network are layered according to the rank of their corresponding solution with respect to non-dominated sorting [START_REF] Goldberg | Genetic Algorithms in Search, Optimization, and Machine Learning[END_REF]. Also, methods that extract informative network features from the analysis of the interactions between the decision variables have been proposed for characterizing the behavior of multi-objective estimation of distribution algorithms [START_REF] Santana | Network measures for information extraction in evolutionary algorithms[END_REF] and other MOEAs [START_REF] Martins | Analysis of Bayesian network learning techniques for a hybrid multi-objective Bayesian estimation of distribution algorithm: a case study on MNK landscape[END_REF][START_REF] Santana | Mining probabilistic models learned by EDAs in the optimization of multi-objective problems[END_REF].

MULTI-OBJECTIVE NK LANDSCAPES WITH HETEROGENEOUS OBJECTIVES

The parameter 𝐾 of an NK model critically influences the characteristics of the landscape: the ruggedness and the number of local optima typically increase with 𝐾. Therefore, we propose to use MNK models, where each objective is defined using a different parameter 𝐾, as a framework to investigate heterogeneous objectives. Such a model allows us to evaluate how the difference in the amount of ruggedness among the single-objective NK models translates into the difficulty of the MOP. We remark that previous research on MNK landscapes has focused on problems where all the objectives share the same parameter 𝐾.

Although our analysis could be extended to an arbitrary number of objectives, in this paper we focus on bi-objective MNK landscapes. Furthermore, we present the analysis from the perspective of optimization strategies based on a Pareto dominant approach.

Constructing Heterogeneous MNK Landscapes

Using 𝐾 as a way to define the amount of heterogeneity is a simple way to control the characteristics of the problem. However, another important question is the way in which the objectives are related.

In this section, we present the way individual objectives have been generated for the heterogeneous MNK model. We focus on an exhaustive exploration of a small, enumerable, heterogeneous MNK landscape. To reduce the effect that the variability of the variable interactions' structure may have in the characteristics of the landscape, we start from an initial instance with 𝐾 = 𝐾 𝑖𝑛𝑖𝑡 and progressively construct instances with lower 𝐾 values. This process is called NK model reduction and can be implemented in different ways. For the experiments shown in this paper, we select, for the interactions of variable 𝑥 𝑖 , the interacting variable with the smallest mean squared distance between each possible configuration of the other variables, that is:

𝑥 𝑗 =arg min 𝑥 𝑗 ∈Π (𝑥 𝑖 ) x ∈𝑥 𝑖 ∪Π (𝑥 𝑖 ) /𝑥 𝑗 (𝑓𝑖 (𝑥 𝑖 , x∪𝑥 𝑗 =1)-𝑓 𝑖 (𝑥 𝑖 , x∪𝑥 𝑗 =0)) 2 (2)
where x ∈ 𝑋 𝑖 ∪ Π(𝑋 𝑖 )/𝑋 𝑗 represents a configuration of the interactions of 𝑖 excluding 𝑗, and 𝑓 𝑖 ( x ∈ 𝑋 𝑖 ∪ Π(𝑋 𝑖 )/𝑋 𝑗 ) is the value of function 𝑓 𝑖 for a particular configuration of the interacting variables.

Once the interacting variable to be removed is selected, a new function 𝑓 𝑖 with 2 𝐾-1 values is randomly generated from a uniform distribution. This process guarantees that, for any pair of NK models (𝑁 𝐾 1 , 𝑁 𝐾 2 ) generated from a common initial instance and such that 𝐾 1 < 𝐾 2 , the Π 𝑁 𝐾 1 (𝑋 𝑗 ) variable interactions of a variable 𝑋 𝑗 are a subset of Π 𝑁 𝐾 2 (𝑋 𝑗 ).

Characterizing Heterogeneity

Considering bi-objective MNK landscapes, benchmark parameters 𝐾 1 and 𝐾 2 configure the heterogeneity of the instance. Using 𝐾 1 = 𝐾 2 to construct a homogeneous instance. We define the total degree of variable interaction as 𝑇 = 𝐾 1 + 𝐾 2 , and the degree of heterogeneity as 𝐷 = |𝐾 2 -𝐾 1 |. It is important to notice that, given 𝑇 and 𝐷, the sum 𝑇 + 𝐷 = max(𝐾 1 , 𝐾 2 ) • 2 must be even in order to ensure feasibility, resulting in only two possible symmetrical pairs (𝐾 1 , 𝐾 2 ). Fixing 𝐾 1 < 𝐾 2 reduce the number of possible configurations to one, therefore the total degree of variable interactions 𝑇 and the degree of heterogeneity 𝐷 are two alternative parameters to characterize MNK landscapes. To observe the impact of heterogeneity on MNK landscapes, we consider the following scenarios:

• Fixed total degree of variable interactions (𝑇 ),

• Fixed degree of heterogeneity (𝐷),

• Increasing heterogeneity without fixing 𝑇 nor 𝐷.

RANK-ANNOTATED NEIGHBORHOOD NETWORK WITH LABELED LOCAL OPTIMAL SOLUTION

In order to evaluate how heterogeneous objectives influence the properties of the landscape, we focus on metrics that capture the composition of the search space according to the rank of solutions with respect to non-dominated sorting. These metrics attempt to capture the way in which solutions are ranked, and the likelihood for solutions of one rank to have neighbors with a better rank.

Constructing the Network

We define the neighborhood structure N as the usual 1-bit flip neighborhood structure; i.e. two solutions are neighbors if their Hamming distance is 1. A Pareto local optimal (PLO) solution is a solution 𝑥 ∈ 𝑋 for which there does not exist any neighboring solution 𝑥 ′ ∈ N (𝑥) such that 𝑥 is dominated by 𝑥 ′ [START_REF] Paquete | On local optima in multiobjective combinatorial optimization problems[END_REF]. Similarly, we use SLO 𝑖 to refer to a single-objective local optima with respect to objective 𝐹 𝑖 .

Solutions in a given rank can be classified into two groups: PLO solutions and solutions that are not locally Pareto optimal. Furthermore, we propose to further classify PLO solutions according to their local optimality for each separate objective. For bi-objective problems, we consider the following five types of solutions:

(𝑡 1 ) PLO solutions that are SLO 1 , (𝑡 2 ) PLO solutions that are SLO 2 , (𝑡 3 ) PLO solutions that are both SLO 1 and SLO 2 , (𝑡 4 ) PLO solutions that are not local optima of any objective, (𝑡 5 ) Solutions that are not PLO.

By construction of MNK landscapes, we assume there is no equivalent neighboring solutions, and therefore a SLO 𝑖 for any objective is always a PLO, whereas the opposite is not true. Furthermore, in the previous classification, groups 𝑡 1 and 𝑡 2 exclude elements from 𝑡 3 . The rationale behind this classification is to analyze the relationship between SLO and PLO, and to investigate what is the distribution of PLO and non-PLO solutions among solutions with the best ranks. Computing the SLO of a single-objective problem requires the evaluation of only one fitness function. Therefore, it is worth investigating whether the knowledge about the SLO of the objectives that comprise a MOP can be exploited to design more efficient strategies for MOPs. Similarly, an important question is how non-PLO solutions are distributed among the different ranks. While PLO have received most of the attention in the analysis of MOPs, non-PLO solutions can contribute to promising PF approximations.

Since our focus is on Pareto dominant approaches, solutions are further analyzed by organizing them according to non-dominated sorting, as in [START_REF] Garrett | Plateau connection structure and multiobjective metaheuristic performance[END_REF]. As such, non-dominated solutions are ranked 1. For this analysis, we assume that the complete landscape of the MOP is known. This is equivalent to an exhaustive enumeration of the solution space that allows us to associate each solution with its rank. The introduced classification of solutions is used to annotate all the solutions for each rank. As a convenient graphical representation of the search space, we use neighborhood networks augmented with labels that identify the type of solutions in terms of local optimality. The neighborhood network is stratified in layers where all solutions that belong to the same rank are located in the same layer.

Fig. 1 shows an example of a rank-annotated neighborhood network with labeled solutions. In the figure, PLO solutions are represented as circles and non-PLO solutions as squares. The different types of PLO solutions are represented with different colors. To ease the visualization, the example displays only some of the connections between the vertices. Notice that edges indicating the neighborhood relationships have the same interpretation than in the neighborhood network, but they can also be read in terms of rank improvement (or deterioration) when visiting the neighbors of a solution.

Landscape Characteristics

As previously discussed in Section 2.3, in multi-objective optimization, much research effort has been put to quantify the impact of the number of objectives, their degree of non-linearity, and the correlation among them on general-properties from the landscape. We summarize in Table 1 some existing metrics and introduce additional ones derived from the analysis of the rank-annotated neighborhood network. From the rank-annotated neighborhood The proportions are computed by first inspecting the ranks of all the neighbors for all solutions of type 𝑡, and then dividing by the total number of neighbors. For example, if there are 5 solutions of type 𝑡 1 and of rank 𝑟 , then the total number of neighbors is 30 and this value is used to normalize the number of connections to solutions with lower, equal, and higher ranks. The same value is used at the time of computing the average rank difference 𝐷 𝑡,𝑟 which is also divided by the highest rank in order to scale the average rank change to values in [-1, 1].

The metrics provide information about the composition of the ranks and allow us to identify which type of PLO solutions are more "promising" in terms of improving the quality of the current solutions. For instance, we can evaluate how likely for solutions with sub-optimal ranks is to move to solutions with better ranks, and whether is is better to try first PLO solutions that are SLO for a single objective or, by contrast, whether PLO solutions that are not SLO should be given priority. Our final goal is to determine how the characteristics of the landscape change with heterogeneity, [START_REF] Knowles | Instance generators and test suites for the multiobjective quadratic assignment problem[END_REF] hv hypervolume covered by solution from R 1

[2] cc proportion of connected components in R 1 [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF] sing proportion of isolated solutions (singletons) in R 1 [START_REF] Paquete | Clusters of non-dominated solutions in multiobjective combinatorial optimization: An experimental analysis[END_REF] lcc prop. size of the largest connected component (lcc) in R 1 [START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF] lcc_hv proportion of hypervolume covered by the lcc in R 1 [START_REF] Liefooghe | Dominance, indicator and decomposition based search for multi-objective QAP: landscape analysis and automated algorithm selection[END_REF] plo proportion of Pareto local optimal solutions (PLO) [ and to what extent these metrics can provide information about which search strategy is more promising for different degrees of heterogeneity.

EXPERIMENTAL ANALYSIS

In this section, we evaluate the effect of the objectives' heterogeneity in the characteristics of MNK landscapes. We first present the experimental framework, including a description of the instance generation process. Then, we present an initial assessment of traditional metrics used to investigate MNK landscapes. Finally, an analysis of the landscapes based on the augmented neighborhood networks is presented.

Benchmark Problems

The benchmark used for our experiments have been generated using the reduction method explained in Section 3.1 with 𝑁 = 14 and 𝐾 ∈ {1, . . . , 12}, such that 𝐾 𝑖𝑛𝑖𝑡 = 12. For each 𝐾 value, we generate 50 random instances, and apply the reduction for each of them obtaining 50 × 12 = 600 (single-objective) NK models. Finally, in order to create bi-objective MNK landscapes, we pair each possible combination of 𝐾 values for each of the original instances, ending up with 50 instances of the 12×11 2 = 66 bi-objective problems. For each bi-objective problem, we report average values among the 50 folds.

Intrinsic Characteristics of Heterogeneous MNK Landscapes

5.2.1 Visual Inspection of the Objective Space. We start by analyzing the impact of heterogeneity on the structure of the objective space. Given a fixed value of 𝑇 = 𝐾 1 + 𝐾 2 = 13, we show in Fig. 2 three exemplary heterogeneous MNK landscapes in which the level of heterogeneity varies, from slightly heterogeneous objectives (left) to highly heterogeneous objectives (right). First of all, we observe that the distribution of objective values is different for different levels of heterogeneity. For nearly homogeneous objectives (left), the ellipse surrounding the objective space looks like a circle, where the range of objective values is the same for both objectives, similar to what we observe for homogeneous MNK landscapes [START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF]. However, as the degree of heterogeneity increases, the position of solutions in the objective space seems to squeeze for the objective with a smaller 𝐾. We attribute this to the underlying distribution of objective values for NK models, for which the range is known to increase with K.

Fig. 2 additionally highlights Pareto optimal solutions, SLO for each objective and remaining PLO in different colors. We observe that the PF size and the number of PLO solutions do not seem to significantly fluctuate with the level of heterogeneity. However, the more homogeneous the objectives, the better distributed the SLO among both objectives (left). By contrast, for highly heterogeneous objectives (right), almost all local optima are with respect to the second objective (𝑓 2 ), resulting in a proportional increase of PLO solutions on top of the objective space. 5.2.2 Structure of the Pareto set. Fig. 3 (top right and center) summarizes some metrics related to the structure of the Pareto set for all considered heterogeneous bi-objective problems: the number of solutions with the best rank (solutions from R 1 ), and the number of connected components. The analysis of the figure reveals that the number of Pareto optimal solutions does not vary significantly among the considered problems. For a fixed 𝐾 1 , it seems to slightly decrease with 𝐾 2 > 𝐾 1 , which corroborates results from homogeneous MNK landscapes [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF][START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF]. However, this observation does not hold when fixing 𝐾 2 and varying 𝐾 1 < 𝐾 2 . This suggests that the "most-difficult" objective has more impact on the cardinality of the Pareto set than the "easier" objective for MNK landscapes. We also observe that there are slightly fewer Pareto optimal solutions as the degree of heterogeneity among the objectives increases. Similarly, the proportion of supported solutions (supp) and the global hypervolume (hv) do not reveal any significant variations (not reported here due to space restriction). They mostly relate to the shape of the PF, which do not reveal significant changes in the benchmark. We can however notice that hv seems to increase with the number of variable interactions, as for homogeneous MNK landscapes [START_REF] Aguirre | Working principles, behavior, and performance of MOEAs on MNK-landscapes[END_REF].

In terms of connectedness of the Pareto set, it is reported to decrease with the number of variable interactions for homogeneous MNK landscapes in [START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF]. This is what we observe in Fig. 3, where the number of connected components (cc) increases with 𝐾 1 and 𝐾 2 .

Although not reported due to space restriction, we observe that related metrics, such as the proportion of isolated Pareto-optimal solutions, and the size and hypervolume of the largest connected component, vary accordingly. Interestingly, we do not observe any significant change with respect to the difference between 𝐾 1 and 𝐾 2 , suggesting that the connectedness is not impacted by the level of heterogeneity among the objectives.

Finally, we remark that the number of ranks does not seem to be affected by 𝐾 1 and 𝐾 2 , since it exhibits only very minor variations (not reported here due to space restriction). These initial observations show that, as far as the structure of the Pareto set is concerned, the information that can be extracted from the way in which heterogeneity influences problem difficulty is limited.

Multi-modality.

Let us now analyze how the non-linearity of objectives, and their level of heterogeneity influences the multimodality of the landscape. In Fig. 3 (top left, and bottom), we report additional metric values related to the number of local optimal solutions for the 50 folds of each pair of (𝐾 1 , 𝐾 2 ) values. Due to space restriction, only a subset of metrics are reported.

For single-objective NK landscapes, it is expected that higher 𝐾 values result in more local optima [START_REF] Kauffman | Origins of Order[END_REF]. For homogeneous MNK landscapes, the number of PLO solutions is also known to increase with 𝐾 [START_REF] Verel | On the structure of multiobjective combinatorial search space:MNK-landscapes with correlated objectives[END_REF]. Overall, all metrics related to the multi-modality show a regular variation according to both 𝐾 values : the larger 𝐾, the rugger the landscape. The number of PLO solutions (plo, top right) varies significantly and, as expected, it clearly depends on the value of 𝐾 for both objectives. However, the number of PLO solutions remains constant whatever the level of heterogeneity among the objectives. Similarly, the mean proportion of SLO per objective (mean_slo, bottom left) significantly increases along 𝐾 1 and 𝐾 2 . Our results reveal that the variation can be significant. As already illustrated in Fig. 2, the number of SLO 𝑖 for objective 𝐹 𝑖 increases with 𝐾 𝑖 , and so does its proportion among PLO solutions. Interestingly, we also observe that the latter increases a bit more significantly as objectives get more heterogeneous, that is when there is a large difference between 𝐾 1 and 𝐾 2 (see, e.g., p_slo_2, bottom center).

The metrics discussed above, such as plo or mean_slo but also cc, mostly capture this general characteristic of the landscape. Similar to the number of PLO solutions, the number of PLO solutions that are not SLO for any objective plo_not_slo increases with the total degree of interactions 𝑇 . However, its normalized variant p_plo_not_slo decreases with 𝑇 . This can be interpreted as follows: for smoother landscapes with (𝐾 1 = 1, 𝐾 2 = 2), 5% of PLO solutions are also an SLO with respect to one objective, whereas for highly rugged landscapes with (𝐾 1 = 11, 𝐾 2 = 12), more than 90% of PLO solutions come from SLO solutions with respect to 𝐹 1 or 𝐹 2 .

The total degree of variables interactions 𝑇 captures the general trend of the ruggedness, but fails to grasp the disparity among the objectives. For a low degree of heterogeneity, we already observed in Fig. 2 that SLO are well-balanced between the objectives. To go further, we compute the deviation between the number of SLO for each objective slo_dev, together with its normalized variant, as a proportion among PLO solutions, p_slo_dev. The latter is reported in Fig. 3 (bottom right) for different 𝐾 1 and 𝐾 2 values. Increasing the heterogeneity 𝐷 such as 𝐾 1 < 𝐾 2 result in a strong concentration of SLO on the second objective. In other words, slo_dev increases with the objectives' heterogeneity. Contrary to the proportion of PLO which are not SLO, that decreases with the total degree of interactions 𝑇 = 𝐾 1 + 𝐾 2 , the proportional variant p_slo_dev is mostly impacted by the degree of heterogeneity 𝐷 = |𝐾 1 -𝐾 2 |.

Correlation Analysis.

To push our analysis further, we report in Fig. 4 the correlation between the 17 considered metrics and (1) the number of variable interactions for the first objective 𝐾 1 , (2) the number of variable interactions for the second objective 𝐾 2 , (3) the total number of variable interactions for both objectives 𝑇 = 𝐾 1 + 𝐾 2 , and (4) the degree of heterogeneity among the objectives 𝐷 = |𝐾 1 -𝐾 2 |. The correlations reveal five categories of observations. Firstly, as expected, the number of SLO (slo_i) and their proportion within PLO (p_slo_i) are highly correlated with 𝐾 𝑖 , forming two small clusters in the center. The latter is normalized by the number of PLO solutions, which results in an anti-correlated relationship between p_slo_i and the cluster with 𝐾 of the other objective. The number of PLO (plo), the average number of SLO per objective (mean_slo) and the number of PLO which are not SLO (plo_not_slo) are all strongly correlated with 𝑇 = 𝐾 1 + 𝐾 2 . So are the proportion of isolated solutions (sing) and the proportion of connected components (cc) among Pareto optimal solutions, forming all together a large cluster in the bottom right of the figure. This cluster is related to both previous clusters, where 𝐾 1 and 𝐾 2 appear. Two metrics are in conflict with the ones above : as expected, the proportion of PLO which are not SLO (p_plo_not_slo) decreases with 𝑇 = 𝐾 1 + 𝐾 2 . Similarly, the size of the largest connected component (lcc) follows an opposite trend to the number of singletons (sing) and the proportion of connected component (cc). In a fourth cluster, we observe that the deviation of SLO (slo_dev) and its proportion (p_slo_dev) are strongly impacted by the degree of heterogeneity 𝐷 = |𝐾 1 -𝐾 2 |, when slo_dev also reveals a medium correlation with other benchmark parameters (𝐾 1 , 𝐾 2 , and 𝑇 = 𝐾 1 + 𝐾 2 ). As such, for the same degree of heterogeneity, the range of the deviation of SLO for each objective is impacted by the total degree of variable interactions. Finally, as previously mentioned, the proportion of Pareto-optimal solutions and of supported solutions (supp), together with the hypervolume of the largest connected component (lcc_hv) do not change much with benchmark parameters.

Differential Analysis of PLO and non-PLO Solutions

The analysis presented in the previous sections shows that a repertoire of traditional metrics sheds light on how the choice of the parameters 𝐾 1 and 𝐾 2 shapes the characteristics of heterogeneous MNK landscapes. In this section, we approach the analysis of heterogeneous MNK landscapes using metrics derived from the rankannotated neighborhood network. These metrics simultaneously provide two perspectives of heterogeneity: the differential role played by the rank, and the distribution of solutions across ranks. We focus on the three scenarios described in Section 3.2. For each scenario, we select three heterogeneous MNK instances and compute the metrics. Fig. 5 shows, for the 9 instances, and for the set of the first 100 ranks, the proportion of neighbors from solutions in each rank that belong to better ranks. This metric provides a measure of how likely it is to find better solutions by locally exploring solutions of a given rank. In Fig. 5, the first three rows represent the scenario of increasing heterogeneity without fixed 𝑇 and 𝐷, the three rows in the middle correspond to a fixed total degree of variable interactions, and the bottom three rows correspond to a fixed degree of heterogeneity.

The analysis of the first scenario in Fig. 5 (left) reveals that, as the total amount of interactions increases, it becomes less likely to improve (jump to a better rank) from a PLO with the best ranks. For the first 50 ranks, the probability of finding better neighbors for (𝐾 1 = 1, 𝐾 2 = 12) is higher than for (𝐾 1 = 11, 𝐾 2 = 12). It is not only that there are more local optima for the second objective, but also that they are distributed in a different way along the ranks. The analysis of the proportion of better neighbors for SLO 1 and SLO 2 show that, in the first ranks, it seems easier to improve solutions from SLO 2 than those for SLO 1 . Notice, that while the second objective is fixed (𝐾 2 = 12), the probability of improving a SLO 2 of the different ranks change as both the heterogeneity and the total degree of variable interactions increase.

When the total amount of interactions 𝑇 = 𝐾 1 + 𝐾 2 is fixed, the proportion of PLO neighbors with a better rank show significant differences for the three MNK instances considered. So does the proportion of SLO 1 . For (𝐾 1 = 5, 𝐾 2 = 8), it is very difficult to find, in the first ranks, better neighbors for SLO 1 . However, as 𝐷 = |𝐾 1 -𝐾 2 | increases, the probability also increases.

When the degree of heterogeneity 𝐷 = |𝐾 1 -𝐾 2 | is fixed, the proportion of neighbors in the first ranks also shows differences among the instances. The increase is related to the increase of the total amount of variable interactions 𝑇 . A common characteristic for both the second and third scenarios is that the probability to improve SLO 1 solutions is limited to a relatively short range of rank values.

CONCLUSIONS

In multi-objective optimization, the heterogeneity among the objectives can be manifested in a variety of ways. One of these ways is the different amount of ruggedness of the landscape induced by each objective. Most previous research on combinatorial MOPs, and in particular on MNK landscapes, have focused on homogeneous objectives. Furthermore, the landscape analysis from the literature of combinatorial MOPs have neither been conceived taking into consideration heterogeneity nor tested on this type of problems. The work presented in this paper aims at providing a better understanding of the impact of heterogeneity on combinatorial MOPs. Our main contribution consists in introducing heterogeneous MNK landscapes with different values of 𝐾 as a benchmark to investigate this question. We further investigated the impact of heterogeneity by means of existing multi-objective landscape metrics and of new metrics that integrate information about the dominance relation among solutions and their type in terms of local optimality. We also introduced a rank-annotated neighborhood network with labeled local optimal solutions, which allows us to extract a finer description of how heterogeneity influences the dominance relations among solutions. Our detailed analysis about local optima for each objective reveals their impact on heterogeneity in the characteristics of the landscape.

An important direction for future research is how to incorporate information about the heterogeneity, as provided by the problem parameters and metrics, into more efficient optimization methods. In the context of local search algorithms, estimating the heterogeneity of the distribution of local optima could help to take proper decisions on the trade-off of exploration and intensification between all objectives. Different strategies could be followed, such as distributing the search effort with respect to the proportion of local optima per objective. Another way to take the impact of heterogeneity into account is by exploiting the fact that, in different ranks, some types of local optima are more likely to be improved than others.

For the analysis included in this paper, we have considered the case of relatively small 𝐾 values. It is possible to extend the analysis to large values of 𝐾 but in this case, it is not feasible to store the fitness values for all sub-functions. Fitness components could be computed as they are called and combined with a caching approach that store partial configurations already. Similar approaches have been previously applied in evolutionary computation [START_REF] Altenberg | Evolving better representations through selective genome growth[END_REF][START_REF] Altenberg | Handbook of Evolutionary Computation, chapter NK fitness landscapes[END_REF]. Similarly, an aspect that is kept for future work is the investigation of heterogeneous MNK landscapes with more than two objectives. Addressing many-objective problems constitutes a formidable challenge [START_REF] Allmendinger | What if we increase the number of objectives? theoretical and empirical implications for many-objective optimization[END_REF] and heterogeneous MOPs do not seem to be an exception in this respect.

Figure 1 :

 1 Figure 1: Rank-annotated neighborhood network with labeled local optimal solutions.

Figure 2 :

 2 Figure 2: Objective space for three heterogeneous MNK landscapes with fixed total degree of variable interaction 𝑇 = 𝐾 1 +𝐾 2 = 13. The degree of heterogeneity increases from left to right. Red points are Pareto optimal solutions, blue points are SLO w.r.t. the first objective (which are not Pareto optimal), green points are SLO w.r.t. the second objective, orange points are PLO solutions (that are not Pareto optimal, nor SLO for any objective), and grey points are non-local optimal solutions.

Figure 3 :

 3 Figure 3: Statistics computed from heterogeneous MNK landscapes.

Figure 4 :

 4 Figure 4: Spearman correlation between landscape metrics, 𝐾 1 , 𝐾 2 , the total degree of interactions 𝐾 1 + 𝐾 2 , and the degree of heterogeneity |𝐾 1 -𝐾 2 |.

Figure 5 :

 5 Figure 5: Proportion of neighbors with better rank for different PLO solutions in three different scenarios.

Table 1 :

 1 General landscape characteristics extracted from MNK landscapes.

	name description
	PF_size proportion of solutions in R 1 (best rank)	[2, 22]
	supp proportion of supported solutions in R 1

  30] mean_slo prop. single-objective local optima (SLO) per objective [26] slo_i proportion of SLO for objective 𝑓 𝑖 p_slo_i same as above, proportional to the number of PLO plo_not_slo proportion of PLO that are not SLO for any objective p_plo_not_slo same as above, proportional to the number of PLO slo_dev deviation of the number of SLO per objective: | slo_1 -slo_2 | p_slo_dev same as above, proportional to the number of PLO

Although in previous work on MNK landscapes, all 𝐾 𝑖 's were assumed to be the same, we use this notation to emphasize that this does not always have to be the case.

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


Hugo Monzón
email: hugo91@gmail.com

Hernán Aguirre

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

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

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

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

Dynamic Compartmental Models for Algorithm Analysis and Population Size Estimation

Keywords: Theory of computation → Evolutionary algorithms; Empirical study, Working principles of evolutionary computing, Genetic algorithms, Multi-objective optimization, Compartmental Models, Modeling

Dynamic Compartmental Models (DCM) can be used to study the population dynamics of Multi-and Many-objective Optimization Evolutionary Algorithms (MOEAs). These models track the composition of the instantaneous population by grouping them in compartments and capture their behavior in a set of values, creating a compact representation for analysis and comparison of algorithms. Furthermore, the use of DCMs is not limited to analysis, by creating models of the same algorithm with different configurations is possible to extract new models by interpolation, and use them to explore fine-grained configurations lying between the ones used as a base. We illustrate the use of the model on some Multi-and Many-objective algorithms, run on enumerable MNK-Landscapes instances with 6 objectives for the analysis, and 5 objectives when used as a tool to do configuration.

INTRODUCTION

In recent years Multi-and Many-objective Evolutionary Algorithms (MOEAs) have become established tools to solve optimization problems in several areas [START_REF] Coello | An updated survey of evolutionary multiobjective optimization techniques: state of the art and future trends[END_REF]. However, until now their use and development lacks a deeper understanding of their dynamics, what makes them successful sometimes, and what makes them fail. Modeling is a powerful tool that could improve this situation, capturing MOEAs properties in parameters that ease their analysis and could lead towards improvements or new algorithms. Predictions obtained could serve for algorithm selection or exploring configurations.

Dynamic Compartmental Models (DCMs) [START_REF] Monzón | Closed state model for understanding the dynamics of MOEAs[END_REF], are linear models inspired in epidemiological compartmental models [START_REF] Godfrey | Compartmental Models and Their Application[END_REF]. Similarly to the way these models track the health status of individuals in a population, DCMs do the same using the non-domination status of solutions to capture the dynamics of MOEAs. In these models, it is assumed that a MOEA population can be split into two or more non-overlapping groups by some set of rules based on Pareto dominance. The population size always remain constant, while group sizes variate in a linear manner and the model tracks these changes. By doing so, given the compartments composition at any time t, the model can estimate how they will change in time t + 1 and so on.

The key element of the model is that translates algorithm dynamics' in terms of compartments and parameters which can be used for analysis or configuration of algorithms as we will show in this work. For analysis we can look at the relationship between compartments and the parameters, using them to draw conclusions on the algorithm behavior during different stages of the search process. For configuration, we present a methodology that allows sampling the algorithm's configuration space through the models. Each model and its parameters are linked to an algorithm and its configuration. We propose to interpolate the model's parameters to obtain models that represent configurations that lie in between the ones used as samples. To illustrate this concept, here we focus only on population size.

DYNAMIC COMPARTMENTAL MODELS

In this paper we will focus on a three compartmental DCM. The population is split into three non-overlapping groups or compartments. At each generation t, the values x t , y t and z t represent the proportion of the population that belongs to each compartment, while fulfilling at all times x t + y t + z t = 1. If we were to represent the DCM in a time discrete manner, we would arrive at the following equation system.

             x t +1 = (1 -(α + β))x t + ᾱy t + βz t y t +1 = αx t + (1 -( ᾱ + γ ))y t + γ z t z t +1 = βx t + γy t + (1 -( β + γ ))z t x t + y t + z t = 1 (1)
where α and β are coefficients that represent the loss in x t which becomes a gain for y t and z t , respectively. ᾱ and γ represent the loss in y t which becomes a gain for x t and z t , respectively. Finally, β and γ represent the loss in z t .

In this work, we will use a DCM able to track performance by following the discovery of new Pareto Optimal solutions. The first rule separates solutions in two groups, the ones belonging to the Pareto Optimal (PO) set and the ones that do not. For the second rule, we take the first group and check when a solution appeared in the current generation. Solutions only found in this and only this generation belongs to the first compartment: PO Absolutely new (POA). Solutions that have appeared in any generation counting from the initial one to t -1 belong to the second compartment: PO Not new (PON). Finally, the third compartment contains all the solutions that are not part of the PO set: Not PO (NPO).

Having determined what solutions each compartment will track, we need to find a set of parameters that minimizes the difference between the values produced by the model D estimated and the actual measured values D measured that comes from the algorithms runs. One easy way is to turn it into an optimization problem where we want to minimize 1 n n i=1 (D estimated -D measured ) 2 in other words the mean square error mse. We solve this problem by using the Covariance Matrix Adaptation Evolutionary Strategy [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF] (CMA-ES), a single-objective numerical optimizer.

ALGORITHMS AND TEST PROBLEM

The algorithms selected are four representatives Multi-and Manyobjective optimization algorithms. The Non-dominated Sorting Genetic Algorithm-II (NSGA-II), Adaptive ϵ-Sampling ϵ-Hood (AϵSϵH), Indicator Based EA (IBEA) with the hypervolume (HV) indicator, and MOEA based on Decomposition (MOEA/D).

MNK-landscapes [START_REF] Aguirre | Insights on properties of multiobjective MNK-landscapes[END_REF], a multi-objective test problem generator, is used in all experiments. The parameter K controls the ruggedness of the landscape, M is the number of objectives, and N the number of variables. The instances used here have K = 1 bit, N = 20 bits and M = 6, for the analysis example and M = 5 objectives for the population estimation example.

All algorithms were run 30 times, with population size 200 and given 100 generations for the 6 objective instance. As for the 5 objective instance, population size ranges from 150 to 750 in 50 increments, and the number of generations depends on the assigned budget of Function Evaluations (FE) and can be calculated as Gen Max = F E/Pop. The maximum number of FE was set to 10000 simulating a limited budget equivalent to 1% of the whole search space (2 20 total possible values).

ALGORITHM ANALYSIS WITH DCMS 4.1 Quality of the found models

Following the procedure described in Section 2, create a model for each algorithm on the problem instance with 6 objectives. We conducted a visual inspection of the model estimation against the measured data. For this, we take the compartments counts on the initial generation for each of the runs and use the model to estimate how they change from generation 1 to 100. Figure 1 to 3 shows the estimation for the POA, PON and NPO compartments in red, while in black is the measured values from the algorithms runs. From the plots we can see that the models are capturing the changes in each compartment for all algorithms, following the mean of all the runs.

Interpretation of the DCM parameters

Parameters in DCM relate two compartments, and in order to properly interpret them, we need to know that they encapsulate according to the equations. Let us focus α and consider the DCM model with the following compartments: POA, PON, NPO. Here this parameter relates the number of POA and PON solutions. If we were to expand the first two equations of the system (1).

POA t +1 = POA t -αPOA t + • • • PON t +1 = PON t + αPOA t + • • • .
Note that the values of POA and PON at time t +1 change in opposite directions, both proportionally to POA at time t by αPOA t . Hence, a loss in one compartment becomes gain in the other. However, α can take positive or negative values. Thus, when α > 0, the model tells us that the number of POA solutions reduces with time proportionally to its previous value and the same amount increases in PON. On the other hand, when α < 0, the model says that the number of POA solutions increases proportionally to its previous value and the same amount reduces in PON. In the above example, α > 0 can be interpreted as negative feedback of POA on itself, whereas α < 0 a positive feedback of POA on itself.

As an example, we are going to concentrate on two algorithms, AϵSϵH and MOEA/D. We write the equations with the proper set of parameters for each algorithm and simplify them to show how the compartment depends on itself and the two other compartments. This allows us to understand better the overall dynamics of the DCMs. In the following, consider the equations marked with A as the ones for AϵSϵH and with M the ones for MOEA/D.

POA A t +1 = 0.5967POA A t + 0.0489PON A t + 0.0768NPO A t POA M t +1 = 0.9621POA M t -0.0288PON M t + 0.0319NPO M t PON A t +1 = 0.9287PON A t + 0.5305POA A t -0.0061NPO A t PON M t +1 = 0.9684PON M t + 0.0001POA M t + 0.0292NPO M t NPO A t +1 = 0.9293NPO A t -0.1272POA A t + 0.0224PON A t NPO M t +1 = 0.9389NPO M t + 0.0378POA M t + 0.0604PON M t
In 6 objectives, even random sampling is able to find some PO solutions, so it is expected to find some of them in the initial population. According to the data, there are some PO solutions, so for our model they will be absolutely new since it is the initial population, being them less than 1%, PON is 0 since there are no previous generations, and the remaining ones are NPO.

We notice that AϵSϵH seems to retain more than 53% of POA solutions as PON solutions, as indicated by the coefficient value of the POA t component in the PON equation, while for MOEA/D this value is very small, less than 0.01%. Thus MOEA/D is dropping substantially more solutions than AϵSϵH. However, note also in POA equations that POA t coefficient is 0.59 for AϵSϵH and 0.96 for MOEA/D, which indicates that MOEA/D finds many more Pareto absolutely new solutions than AϵSϵH. It is worth noting that POA solutions in this instance are found mostly from previous POA solutions, rather than from PON or NPO solutions. Note that the analysis done here through the parameters can also be done by inspecting the plots in Figures 123.

ESTIMATION OF POPULATION SIZES 5.1 Using splines to interpolate models

Our models are composed of six parameters (α, β, γ , ᾱ, β, γ ), which are linked to a particular combination of population size, algorithm and problem, so for different configurations we would have other set of parameters. We could then interpolate the values in the space of parameters, and discover this way other configurations, however for this work we only focus on the population size. In order to interpolate the models, we make a spline for each parameter, taking as data points the pair population size and parameter value.

The first step is selecting a range to explore and some sample population sizes on the range. We run the algorithms to generate the data to construct the models. Next, we fit the models to the data and get the parameters. With this data, we now create splines for each parameter. Finally, we use the splines to explore the range, generating models in between population sizes used as samples.

Here for our example, the problem was a 5 objective MNK-Landscape instance, and we decided to explore the range {150, 750} with {150, 350, 550, 750} as sample points. We created then, new points every 50, exploring 3 new population sizes between the ones we actually had to run the algorithm. We pay the full cost of creating four models, and obtain nine more very cheaply.

Selecting Population Size under a Budget

To do selection, we consider the number of Accumulated PO (Acc. PO) solutions obtained by the algorithm with a budget of 10000 FE. We obtain this value by accumulating the estimated POA count from the model associated with each population size. For verification, we analyze the estimated Acc. PO at the end of the budget from fitted models and interpolated models, and compare then with the actual measured values that we obtain running the algorithm. We want to see if the models will allow us to discover which interval between the sampled population shows promise.

Figure 4 shows on the left for sample population sizes, the average number of Acc. PO solutions at the end of 10000 FE for 30 runs of the algorithm. This sample population sizes were used as knots for the spline. On the center we include a similar plot with values produced by the fitted models in red, while the ones from interpolated models are in black. Finally on the right side we have now the measured values for all population sizes in the range. Red ones represent the population sizes in the range used as knots for the spline, and the ones in blue the interpolation ones. We use this plot to check against the trends found by the models. In all plots the error bars show the 95% confidence interval for each mean.

Looking at Figure 4 left and center plots, we can visually compare the fitted model estimation to the measured values for the population sizes used as knots. We see that for all sizes except 150 the fitted model gets very close to the measured average, in the particular case of 750 is even a little higher. Nevertheless, as we will see, even with an underestimating model as a knot, the trend is still correctly replicated. Now we shift the focus to the center and right figures, mainly the estimations of our models obtained through interpolation. We clearly see the influence of our first knot, 150, pulling down the estimations of the interpolated models between [150, 350]. Nevertheless in this section the growing trend is still correctly maintained when compared to the plot on the right with all the measured results.

Moving to the next interval of [350, 550], we find that our models estimations get closer to the measured results, and also follows the trend detecting that 400 and 500 obtain better values than 450. Finally in the last interval of [550, 750], the trend is not correctly replicated, as 700 places a bit higher than 650, which according to the right side of the figure should be below. Since the fitted model of 750 actually overestimates a little in this part, it seems to pull the interval making 700 seem a better option. However, the relative position of the knots points is still correct, as two of them place higher than 750 as it should be, and neither of the group is higher than the knot on left, 550.

Summarizing, if we were to trust only the model results, the information obtained from the interpolated models will guide us towards choosing a population size for our budget better than the one used in our sampling. If we partially trust the model, it still will guide us to a region worth exploring in detail.

CONCLUSION

In this work we presented some possible uses of Dynamic Compartmental Models, with an example as an analysis tool and also a small step forward towards algorithm configuration with an example that estimates population size given a budget. The key in both examples is the ability of parameters to capture the complex dynamics in a numerical way which opens many possibilities, specially for algorithm configuration.

In future works, we would like to propose new sets of rules that does not involve PO solutions, to be able to test the DCMs on real world problems and test the limits of this new tool.

Figure 1 :

 1 Figure 1: Model estimation vs measured PO Absolutely new data.

Figure 2 :

 2 Figure 2: Model estimation vs measured PO Not new data.

Figure 3 :

 3 Figure 3: Model estimation vs measured Non PO data.

Figure 4 :

 4 Figure 4: [Budget: 10000 FE] Average Accumulated PO over population size. Left: Sampled population sizes used as knots for the spline. Center: Model estimation for the knots and interpolated population sizes. Right: Measured values for all population sizes.


Antoine Castillon
email: antoine.castillon.etu@univ-lille.fr

Julien Baste
email: julien.baste@univ-lille.fr

Hamida Seba
email: hamida.seba@univ-lyon1.fr

Mohammed Haddad
email: mohammed.haddad@univ-lyon1.fr

Quasi-Clique Mining for Graph Summarization

come   L'archive ouverte pluridisciplinaire

Introduction

Graphs are a powerful modelling tool allowing to focus on the interactions between various entities. Today, several large-scale systems such as social networks or the link structure of the World Wide Web involve graphs with millions and even billions of nodes and edges [START_REF] Boldi | The webgraph framework i: Compression techniques[END_REF]. The analysis of such graphs can prove to be very difficult given the huge amount of information contained in them. We note that a graph G = (V, E) is represented by a set of vertices V and a set of edges E that link the vertices.

The goal of graph summarizing (or graph compression) is to provide, given an input graph, a smaller summary [START_REF] Liu | Graph summarization methods and applications: A survey[END_REF]. Depending of the situation, this summary can assume different roles. Either it exists uniquely to save storage space and can be decompressed to obtain the original graph with or without loss of information, or it can be used directly to obtain information on the input graph such as paths, neighbourhoods and clusters. When summarizing a graph, one can focus on the optimisation of the size of the summary introducing a size threshold [START_REF] Ahmad | Ssag: Summarization and sparsification of attributed graphs[END_REF][START_REF] Lee | Ssumm: Sparse summarization of massive graphs[END_REF] or an information loss threshold [START_REF] Navlakha | Graph summarization with bounded error[END_REF]. Others adapt the summary so as to optimise some kind of queries [START_REF] Ahmad | Ssag: Summarization and sparsification of attributed graphs[END_REF][START_REF] Daniel | Compact representations of separable graphs[END_REF].

Even if these compression methods perform well in term of execution time and compression rate, they do not retain all the information available in the input graph, especially the dense components. Indeed, with these methods, such information is usually not accessible in the summary and requires a decompression in order to be retrieved. However, it appears that these dense components are essential in the analysis of many real life networks, such as social networks or protein connection networks [START_REF] Bader | An automated method for finding molecular complexes in large protein interaction networks[END_REF][START_REF] Truong | Graph methods for social network analysis[END_REF][START_REF] Ucar | Improving functional modularity in protein-protein interactions graphs using hub-induced subgraphs[END_REF]. Hence, in this paper, we focus on how well a summary retain information about the dense subgraphs of its original graph. This property can be measured by the concept of visibility introduced in [START_REF] Wang | Redundancy-aware maximal cliques[END_REF] and defined in Section 3. In [START_REF] Wang | Dense subgraphs summarization: An efficient way to summarize large scale graphs by super nodes[END_REF], the authors propose the Dense Subgraph Summarization method (DSS), a lossless compression scheme which addresses this issue, allowing a direct access to the dense subgraphs in the summary. However, this method requires the computation of some dense components of the Graph. In [START_REF] Wang | Dense subgraphs summarization: An efficient way to summarize large scale graphs by super nodes[END_REF], the dense components are formally defined as quasi-cliques (cf. Definition 1) and the computation is done by solving the Quasi-Clique Enumeration Problem (cf. Definition 1).

[γ quasi-clique] Given G = (V, E), and γ ∈ [0, 1], a subset of vertices

Q ⊆ V is called a γ quasi-clique if for all v ∈ V , |N (v) ∩ Q| ≥ γ(|Q| -1)
, where N (v) represents the neighbours of v.

[Quasi-Clique Enumeration Problem] Given G = (V, E), a density threshold γ ∈ [0, 1] and a size threshold m, the Quasi-Clique Enumeration Problem (QCE) consists of finding all maximal γ quasicliques in G of size greater than m.

DSS [START_REF] Wang | Dense subgraphs summarization: An efficient way to summarize large scale graphs by super nodes[END_REF] uses the Quick algorithm [START_REF] Liu | Effective pruning techniques for mining quasi-cliques[END_REF] to solve the QCE problem. This algorithm performs a Depth-First Search on the solution space tree along with several pruning techniques.

In this paper we introduce new pruning techniques which, when added to the Quick algorithm, allow an important speed up. Also we describe how to adapt the search for quasi-cliques to this summarization approach by solving best suited relaxations of the Quasi-Clique Enumeration Problem. When compared to the DSS method alone the adding of our new techniques allows an important improvement of the performances both in time and compression rates and more importantly it retains most of the information regarding dense subgraphs within the summary.

The remainder of this document is organised as follows. Section 2 describes the DSS compression method and the Quick algorithm. Section 3 introduces the new pruning techniques and presents a relaxation of the quasi-clique enumeration problem. Section 4 compares the performances of the different quasi-clique mining schemes and their application to dense subgraph based compression.

Background and Motivation

Dense Subgraph Summarization

Unlike other compression methods [START_REF] Navlakha | Graph summarization with bounded error[END_REF], the DSS method [START_REF] Wang | Dense subgraphs summarization: An efficient way to summarize large scale graphs by super nodes[END_REF] is not designed to optimise the size of the summary. Thus it is, generally speaking, not the most efficient in terms of compression rates. The purpose of this method is to have in the summary a direct access to most of the original graph information, especially its dense subgraphs. We will discuss in Section 3.2 how to measure the retention of the dense subgraphs in the summary using the notion of visibility.

The DSS compression scheme works as follows: a graph G is compressed into a supergraph SG in several steps. First the dense subgraphs of G are computed using a quasi-clique enumeration algorithm (detailed in Section 2.2). Second the dense subgraphs are compressed into supernodes: each supernode Sn i records all the nodes of G it contains in the set node set i . Also to ensure a lossless compression, Sn i keeps track of all non-edges inside node set i in the set loss edges i (an example of the supernodes transformation is given in Figure 1). Third, a superedge is added between two supernodes with overlapping node sets (as shown in Figure 2). A superedge records its two endpoints along with the overlapping part of the two node sets in a set connect set. Also, to ensure a lossless compression scheme, a set remaining edges is added to store the edges that are not already covered by the supernodes. However, this process can create redundant superedges for instance the supernodes {1, 2, 3}, {1, 2, 4} and {1, 2, 5} are connected with three superedges but, since they share the same overlapping part, only two superedges are required to represent the connection between them. In order to process these redundant edges a tracker is added. The tracker prevents the creation of a superedge between two supernodes already connected to superedges with the same connect set.

Quasi-Clique Enumeration

As explained in the previous section, the DSS method requires a dense sub-graph mining usually done by solving QCE problem (cf. Definition 1). This problem is well known to be difficult (the associated decision problem being NP-hard [START_REF] Baril | Hardness and tractability of the γ-complete subgraph problem[END_REF]), and have already been studied a lot trough several aspects: exhaustive enumerations [START_REF] Liu | Effective pruning techniques for mining quasi-cliques[END_REF][START_REF] Zeng | Out-of-core coherent closed quasiclique mining from large dense graph databases[END_REF], or only top-k enumeration [START_REF] Sanei-Mehri | Enumerating top-k quasi-cliques[END_REF][START_REF] Yuan | Diversified top-k clique search[END_REF]. We focus on the exhaustive enumeration and its most classic method: the Quick algorithm [START_REF] Liu | Effective pruning techniques for mining quasi-cliques[END_REF]. Since there is no known polynomial time algorithm, the method to tackle this problem is to review every subsets of V . Given a graph G, a density threshold γ and a size threshold m, the Quick algorithm performs a depth-first exploration of the solution space tree (as in Figure 3). In order to avoid the exponential time complexity and to speed up the computation, the algorithm prunes the tree: branches which cannot lead to a new solution are skipped. During the search, the algorithm keeps in memory X the current subset and a set cand formed by vertices outside X which are the next candidates to add to X. The pruning occurs according to these sets and several criteria, and actually corresponds to the removal of vertices from cand. Some of the most efficient criteria used in Quick are listed here.

Diameter:

If γ ≥ 0.5, then the diameter1 of a quasi-clique is at most 2. So, during the execution one can restrict cand to the nodes that are at most two hops apart of all nodes in X.

Degree threshold:

All vertices in a quasi-clique Y have at least ⌈γ(|Y | -1)⌉ neighbours in Y . So, during the execution one can abort if ∃u ∈ X such that d X∪cand (u) < γ(|X| -1), and restrict cand to the u such that d X∪cand (u) ≥ γ(|X| -1).

Upper/Lower bounds:

Given X and cand one can compute L X (resp. U X ) a lower (resp. upper) bound on the number of vertices of cand that must be added to X concurrently to form a quasi-clique. Given these numbers one can prune more efficiently2 .

Critical vertices:

If there is v ∈ X such that d X∪cand (v) = ⌈γ(|X| + L X -1)⌉ then we can add all neighbours of v in X.

Our Contributions

New Pruning Techniques

While the Quick algorithm is efficient in the general case, it is not hard to point out cases where it is not. One can, for instance, take the example depicted in Figure 4. The reason that explains why Quick in not efficient on this example is because all pruning criteria of the algorithm are based only on the vertices and their degree. In the case where each vertex is contained in a quasi-clique, the pruning techniques cannot be properly applied and the algorithm proceeds to the exploration of almost the whole solution space tree. Given a graph G, γ ∈ [0, 1], Y a γ quasi-clique it holds that:

∀X ⊆ Y : u∈X N Y (u) ≥ |X|⌈γ(|Y | -1)⌉ -(|X| -1)|Y |.
We provide a proof of the lemma by induction on |X|.

• If |X| = 1: it is the definition of a γ quasi-clique. • For X ⊆ Y such that |X| < |Y |. Let v ∈ Y \ X, N Y (v) ∩ u∈X N Y (u) = N Y (v) \ V \ u∈X N Y (u) , so by induction on X, u∈X∪{v} N Y (u) ≥ ⌈γ(|Y | -1)⌉ -|Y | + |X|⌈γ(|Y | -1)⌉ -(|X| -1)|Y | ≥ (|X| + 1)⌈γ(|Y | -1)⌉ -|X||Y |.
• The results holds for all X ⊆ Y and the lemma follows.

Using Lemma 3.1, we propose the following rule: given the density and size threshold γ and m, if {u, v} ∈ E is such that u and v share strictly less than 2⌈γ(m -1)⌉ -m common neighbours then we can remove {u, v} from E. Applying this rule does not change the quasi-cliques of G because, thanks to Lemma 3.1 such edges are not contained in any quasi-cliques. This pre-processing is depicted in Algorithm 1.

Algorithm 1 Pre processing(G, γ, m) repeat delete all u ∈ V s.t. d(u) < ⌈γ(m -1)⌉ delete all {u, v} ∈ E s.t. |N (u) ∩ N (v)| < 2⌈γ(m -1)
⌉ -m until G has not been modified during the last loop Actually, this rule does not use the lemma at its full potential. Indeed the lemma provides a much more general result which could be used to add a new pruning technique based on forbidden sets of vertices. For instance, when parsing the solution space tree, if the label set X of the current node does not verify the lemma then the branch starting from this node can be pruned right away. However, such pruning technique is very hard to adapt in practice since it would require to test at each node X all the subsets of X. Also, as discussed in Section 4, the performances of the pre-processing alone are already very promising. Hence, we limit ourselves to the pre-processing algorithm, and a pruning technique making a full use of the lemma could be the object of further research.

Redundancy Aware Maximal Quasi-Cliques

While the Quick algorithm is an effective method to tackle the QCE problem, it is important to note that this problem is not perfectly suited to the DSS method. Indeed, having access to all quasi-cliques can create redundancy in the summary and impact the performances. Figure 5 gives such an example: the graph G has four very similar 0.85 quasi-cliques, so, in the summary the fact that [ [START_REF] Ahmad | Ssag: Summarization and sparsification of attributed graphs[END_REF][START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF]] forms a clique is stored four times (one per supernodes). This redundancy does not appear with only one quasi-clique which leads to a better compression.

In order to detect and avoid redundant quasi-cliques, [START_REF] Wang | Redundancy-aware maximal cliques[END_REF] introduce the notion of visibility 3 . The visibility of a maximal quasi-clique Q with respect to a set of maximal quasi-cliques S is a value between Sn 1 = ([ [START_REF] Ahmad | Ssag: Summarization and sparsification of attributed graphs[END_REF][START_REF] Liu | Graph summarization methods and applications: A survey[END_REF]], ∅)

Sn 2 = ([[1, 8]] ∪ {11, 12}, ∅) Sn 3 = ([[1, 8]] ∪ {13, 14}, ∅) Sn 4 = ([[1, 8]] ∪ {15, 16}, ∅) graph G
dense components of G Summary of size 64 0 and 1 which describes how well Q is represented by the quasi-cliques of S. If the visibility is 1 then Q is in S, on the contrary if the visibility is 0 then Q does not overlap with any quasi-clique of S. [START_REF] Wang | Redundancy-aware maximal cliques[END_REF] also formally introduce a variant of the QCE problem (Definition 3.2) where instead of finding all the quasi-cliques one only have to ensure that all maximal quasi-cliques have a sufficient visibility.

[Visibility] Given a graph G, γ

∈ [0, 1], Q ∈ M γ (G) and S ⊆ M γ (G) (the set of all maximal γ quasi-cliques of G), the S-visibility of Q: V S (Q), is defined as: V S (Q) = max Q ′ ∈S |Q ∩ Q ′ | |Q| . [τ -visible QCE] Given G, γ ∈ [0, 1] and τ ∈ [0, 1], the τ -visible Quasi-Clique Enumeration Problem is to compute S ⊆ M γ (G) such that ∀Q ∈ M γ (G), V S (Q) ≥ τ .
Note that this problem is well suited to the DSS method for a few reasons. First, the purpose of the DSS method is to have a direct access to the dense subgraphs in the summary. A summary obtained by the DSS method with S a solution of the τ -visible QCE problem as input, provides a direct access to all the elements of S which represent well all maximal quasi-cliques. Second, avoiding redundant quasi-cliques improves not only compression rates but, if done during the search, also allows to speed up the computation, pruning branches corresponding to redundant quasi-cliques.

Due to the similarity between this problem and the QCE problem, our approach to solve it is to adapt the Quick algorithm, adding techniques allowing to skip redundant quasi-cliques during the search. The first algorithm we use is Algorithm 2, introduced in [START_REF] Wang | Redundancy-aware maximal cliques[END_REF]. This algorithm adds to an enumeration algorithm, a pruning technique based on visibility. This pruning occurs randomly during the exploration of the tree depending on the visibility of the set of the current vertices, X ∪ cand, with respect to the last quasi-clique found. Due to the structure of the solution space tree, similar sets are usually close to each other in the tree. Thus, redundant quasi-cliques are often skipped. This approach has several advantages. First, as expected the skipping allows this algorithm to run significantly faster than Quick and gives better compression rates (performances are detailed in Section 4). Also, as proven in [START_REF] Wang | Redundancy-aware maximal cliques[END_REF], the set S obtained by this algorithm verifies that for all Q ∈ M γ (G), the expected value of V S (Q) is greater than τ which ensures a good representation of all quasi-cliques of G. Finally, the parameter τ represents a trade-off between the runtime and the visibility, allowing the user to adapt the algorithm according to its preferences. With lower values of τ the algorithm runs faster but ensures less visibility. On the contrary higher values of τ ensure better visibility but slow down the runtime.

Regarding our second approach, we note that given a graph G and two quasi-cliques

Q, Q ′ . If Q and Q ′ are very similar (V {Q} (Q ′ ) is close to 1) then most of the edges in G[Q ′ ] are already covered in G[Q]. On the contrary if Q and Q ′ are very different then G[Q] and G[Q ′ ]
do not share many edges. Hence, to avoid redundancy, after finding the quasi-clique Q one can remove from G all the edges inside G[Q], thus the quasi-cliques of G similar to Q become sparse and are skipped by the pruning techniques of the Quick algorithm. Removing such edges usually does not impact other different quasi-cliques since these Algorithm 2 Quick redundancy aware(G, γ, m, τ ):

Let T be the solution space tree of V (nodes are labelled with X and cand)

Let S = ∅ Let Q ′ = ∅ Depth first search on T X, cand = current node l = |X| + |cand| Prune T with probability 1 -l s |X∩Q ′ | l 4
if T is not pruned pursue the depth first search . . . (same as Quick)

if X is a quasi-clique and with probability

|X| s |X∩Q ′ | |X| Add X to S Q ′ = X return S
quasi-cliques do not share many edges with Q. This method depicted in Algorithm 3 outputs a set of quasi-cliques S which covers all edges contained in quasi-cliques. Similarly both the runtime and the compression rate are better with this method than with the Quick algorithm (Section 4). Even if this algorithm do not ensure any lower bound on the visibility (contrary to Algorithm 2), the experimental results presented in Section 4 show that the visibility obtained is mostly acceptable.

Algorithm 3 Quick delete covered edges(G, γ, m) run Quick(G, γ, m) each time a quasi-clique X is found: remove from G all edges of G[X]
4 Performances

We compare the performances of four different quasi-clique mining algorithms: the Quick algorithm presented in Section 2.2, Algorithms 2 and 3 presented in Section 3.2 and finally a greedy quasi-clique mining algorithm, i.e., Algorithm 4, which serves as a basis of comparison. These four algorithms have been implemented in Python 3 and executed on a Windows 10 computer with a 2.5 GHz Intel(R) Core (TM) i5 processor and 8 GB RAM. Each algorithm is tested with and without the pre-processing algorithm, namely Algorithm 1, on several types of graphs. Since the DSS method is especially interesting for social networks, we use two types of synthetic graphs which model well social networks: the community graph model, developed in Section 4.1, and the LFR graph model, developed in Section 4.2. In Section 4.4, we also use graphs from social networks such as Facebook, Twitter and Google+. As for the performances we want to measure not only the runtime of the algorithms but also how well they are suited to the DSS method. Hence, three aspects have been retained: the runtime, the size of the summary obtained with the DSS method and the visibility of all quasi-cliques in the input graph.

Algorithm 4 Greedy quasi cliques(G, γ)

Qclqs = ∅ non covered edges = E while non covered edges ̸ = ∅ do Let {u, v} in non covered edges X = {u, v} cand = {w ∈ V \ X|X ∪ {w} is a quasi-clique } while cand ̸ = ∅ do Let z in cand X = X ∪ {z} cand = {w ∈ V \ X|X ∪ {w} is a quasi-clique } end while Qclqs = Qclqs ∪ {X}
remove from non covered edges all edges inside X end while The first model we consider is the community graph model described in Algorithm 5. As input, the user specifies n the number of vertices, k the number of communities, s the size of the communities and two parameters α, β ∈ [0, 1]. Given these inputs the algorithm constructs first the communities and then the graph. The edges inside the communities are added to the graph with probability α, usually close to 1, and the ones outside the communities are added with probability β, usually close to 0. In order to model further social networks we add the possibility for communities to intersect one another and also to have different sizes replacing the parameter s with two parameters s 1 and s 2 respectively the lower and upper bound of the community sizes.

Community Graph Model

On this type of graphs, since the communities are mostly disjoints, the pre-processing algorithm is expected to be very effective. An instance of a community graph is given in Figure 6. Note that the preprocessing algorithm manages to delete all edges outside the communities while edges inside communities are mostly preserved. Also, the communities are separated in different connected components which corresponds to the best case scenario, the quasi-cliques of G = (V, E) are computed in time O(|V | × |E|) (complexity of the pre-processing algorithm).

The performances of the algorithms on community graphs are given in Figures 78. In Figure 8 there is the comparison in time, compression rates and visibility (respectively, first, second and third column) of all algorithms. Note that the Quick algorithm achieves the slowest runtime, the worst compression rates and the best visibility, that is always 1. This is due to the fact that this algorithm actually enumerates all quasi-cliques (which explains the visibility). Also, as discussed in Section 3.2, this causes the algorithm to be slower than the others (which prune more the search tree) and creates redundancy in the summary and thus, leads to worst compression rates. Regarding the other algorithms, it seems there is a trade-off between runtime, compression rate and visibility: the greedy algorithm which achieves the best runtimes and the second best compression rates have the worst visibility (around 0.4 against 0.9 for the others). This trade-off explains why, the DSS method requires complex quasi-clique mining algorithms like Quick (or Quick redundancy aware or Quick delete covered edges) and why it cannot properly work with simple greedy algorithms.

Figure 7 presents the results of the comparison between the performances with and without the use of our pre-processing algorithm, Algorithm 1. It appears that the use of this algorithm lowers significantly the runtimes, up to a factor 7 for values of β greater than 0.4, avoiding a lot of cases like Figure 4 and thus is very useful for speeding up quasi-clique mining algorithms. Also, since this algorithm does not delete edges inside quasi-cliques, the performances in compression rates and visibility are not really impacted by its use, making this algorithm also very interesting for the DSS method.

LFR Graph Model

Even if the results with the community graph model are promising, the model itself is very simplistic. Indeed, for instance two nodes in a community C have mostly the same number of neighbours inside and outside C (≈ α|C| and β(n -|C|)), also the communities can be very similar (as in Figure 6). On the contrary, real networks often have a much more heterogeneous distributions of nodes degree and community sizes [START_REF] Lancichinetti | Benchmark graphs for testing community detection algorithms[END_REF].

The second model we consider, the LFR graph model [START_REF] Lancichinetti | Benchmark graphs for testing community detection algorithms[END_REF], addresses these issues using power laws for the degree and the community size distributions (which models well social networks). Given as input n the number of vertices, τ 1 , τ 2 > 1 the power law exponents of respectively the degree distribution and the community size distribution, µ ∈ [0, 1] and average degree the average degree of the vertices, the LFR benchmark algorithm proceeds as follows. First, the algorithm finds a degree sequence with a power law distribution which respects the average degree value. Then, the algorithm generates the sizes of each community according to another power law and each node is randomly assigned to a community according to the community sizes. Finally, the edges are computed according to the degree of each node and to the community, each node u having The performances of the algorithms on LFR graphs are given in Figures 910. The results obtained are very similar to the results with the community graph model. They show as previously that the use of the pre-processing algorithm improves significantly the runtimes without affecting the compression rates nor the visibility. Additional remarks can be pointed out, for instance with values of µ greater than 0.4 the visibility of all algorithms is constant to 1. This is due to the fact that for these values of µ, the communities are sparse (the density of edges inside them is around 1 -µ) and thus, no quasicliques exist, and with no quasi-cliques, the visibility is always 1. Also, with the same values of µ, the compression rates of all algorithms are the same and mostly stable because with no quasi-cliques the summary graph is empty and the original graph is entirely stored in the set remaining edges. Thus, all summarizations have the same size, the size of the input graph (since, the parameter µ does not affect the number of edges, the size of the input graph is mostly constant). Finally, note that with the Quick and the Quick redundancy aware algorithms, the compression rates first drop for values of µ from 0.2 to 0.3, increases for values of µ in [0.3, 0.4] and then stabilises. This is caused again by redundancies: for small values of µ (around 0.2) the communities are dense and contain a lot of quasi-cliques which creates redundancies (note that with the Quick algorithm redundancies can cause the summary to be larger than the input graph). With larger values of µ (around 0.3), fewer quasi-cliques remain which prevents redundancy. And as explained before with µ larger than 0.4 the input graph does not contain any quasi-clique and thus no compression is achieved.

Social Networks

Since our methods and algorithms have proven to be effective on synthetic graphs and in order to test further their reliability we also evaluate their performances on real social networks. The networks used are available on the Stanford large network dataset collection [START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF]. The performances on social networks are given in Figures 1112. When comparing the performances with and without the use of the pre-processing algorithm, these results prove its utility showing again a significant decrease of the execution while not impacting the compression rates nor the visibility. However, when comparing the different algorithms, they give a more mixed result. Indeed, first the redundancy aware algorithm does not seems to be much more efficient than Quick both in time and compression rates, where the gain should be. Also, even if Algorithm 3 seems efficient regarding execution time and compression rates, the visibility obtained with this algorithm is worst than the others dropping suddenly below 0.8 while the others stay close to 1. Finally, the greedy algorithm seems to perform surprisingly well on these graphs, being the most efficient both in time and compression rate while keeping a very high visibility. However, we think that this high visibility rate is due to the small number of different quasi-cliques existing in the graphs. Indeed, even with 1800 nodes the graphs only contain highly similar quasi-cliques, for instance five quasi-cliques are enough to achieve a visibility of 0.95, and with three of them the visibility stays higher than 0.93. Hence, given the other results, we think that the visibility rate of the greedy algorithm will drop significantly when adding more nodes (and thus more quasi-cliques) to the graph.

Discussion

While the previous results have shown that our pre-processing algorithm and our Quick variants are very effective for the DSS method some social networks present patterns unsolvable with the Quick algorithm (and its variants). For instance, in Figure 13 there is a large sparse quasi-clique which is very problematic for the Quick algorithm and its variants. Even if the density of the quasi-clique is quite low, the large number of vertices prevents the pruning, making all techniques (including the pre-processing) useless. Hence, the Quick algorithm proceeds to the review of every subsets, thus running indefinitely without finding any quasi-clique. Since, the redundancy aware and delete covered edges algorithms run similarly as Quick as long as no quasi-clique is found, the variants have exactly the same issue. Finally, this case cannot be solve by a better management of the arbitrary parameters γ and m. Indeed, as shown in [START_REF] Liu | Effective pruning techniques for mining quasi-cliques[END_REF] the Quick algorithm is always faster with larger values of γ, also only values of γ > 0.5 achieve interesting compression rates with the DSS method (with γ < 0.5 there are more lost edges than compressed edges). Also, in this example, with a size threshold m greater than 121, there is no quasi-clique and thus no compression and with a size threshold m lower than 120 the issue occurs.

Conclusion

In this paper, we focused on graph compression using dense subgraphs and more precisely quasi-cliques, trough the example of the DSS method. This approach requiring a quasi-clique mining algorithm, we show how to improve the already existing quasi-clique enumeration algorithms introducing a very effective pre-processing technique. Also, we adapt and introduce other quasi-clique mining algorithms designed to avoid redundancy (without loosing the purpose of the DSS method) improving performances both in runtime and compression rate. The experimental results show that if our algorithms are very promising, some specific patterns still cannot be handled and require new techniques. This can be the subject of future work.

1 = (node set 1 ,Figure 1 :

 111 Figure 1: A dense node set and the corresponding supernode

2 Se--Figure 2 :

 22 Figure 2: A graph G and the associated supergraph with supernode set {Sn 1 , Sn 2 } and superedge set {Se}.

Figure 3 :

 3 Figure 3: The solution space tree for V = {a, b, c}

Figure 4 :

 4 Figure 4: Example of the inefficiency of the Quick algorithm.

Sn 1 Figure 5 :

 15 Figure 5: Comparison of two summarizations of G, one with all the dense components and the other with only one of them.

Figure 6 :

 6 Figure 6: Left: a community graph; Right: the same graph after the execution of the pre-processing algorithm.

Algorithm 5

 5 Community graph(n, k, s, α, β) G = ([[1, n]], ∅) partition V G into k + 1 sets R, C 1 , ..., C k (|C 1 | = ... = |C k | = s) for i ∈ [[1, k]] dofor {u, v} ⊆ C i do add uv with probability α end for end for for {u, v} ⊆ V G s.t. u and v are not both in any C i do add u, v with probability β end for return G

  Quick algorithm: Quick redundancy aware algorithm: Quick delete covered edges algorithm: Greedy quasi cliques algorithm:

Figure 7 :

 7 Figure 7: Comparative performances of the algorithms with and without the pre-processing on community graphs with n = 250 nodes, k = 8 communities of size 30, α = 0.9 and β ∈ [0, 0.05] corresponds to the values on the x-axis.

Figure 8 :

 8 Figure 8: Performances of all algorithms on community graphs with n = 250 nodes, k = 8 communities of size 30, α = 0.9 and β ∈ [0, 0.05] corresponds to the values on the x-axis.

  (1 -µ) • d(u) neighbours in the same community and µ • d(u) neighbours in other communities. Quick algorithm: Quick redundancy aware algorithm: Quick delete covered edges algorithm: Greedy quasi cliques algorithm:

Figure 9 :

 9 Figure 9: Comparative performances of the algorithms with and without the pre-processing on LFR graphs with n = 300, average degree = 20, 16 communities of size 19 (on average) and µ ∈ [0.2, 0.8] corresponds to the value on the x-axis.

Figure 10 :

 10 Figure 10: Performances of all algorithms on LFR graphs with n = 300, average degree = 20, 16 communities of size 19 (on average) and µ ∈ [0.2, 0.8] corresponds to the value on the x-axis.

Figure 11 :

 11 Figure 11: Performances of all algorithms on social networks with n ∈ [0, 2000] nodes, the value of n is represented by the x-axis.

Figure 12 :

 12 Figure 12: Comparative performances of the algorithms with and without the pre-processing on social networks with n ∈ [0, 2000] nodes, the value of n is represented by the x-axis.

Figure 13 :

 13 Figure 13: A social network from Facebook with 4039 nodes and 88234 edges.

The diameter of a graph is defined as the length of the longest induced path.

For instance in the second pruning technique γ(|X| -1) can be replaced by γ(|X| + L X -1).

The results of[START_REF] Wang | Redundancy-aware maximal cliques[END_REF] hold for cliques but can be easily adapt to quasi-cliques

s : r ∈ [0, 1] → (1 -r)(2 -τ ) (2 -r -τ ), note that s is decreasing, s(0) = 1 and s(1) = 0.

Acknowledgement:

This work is funded by Agence Nationale de la Recherche (ANR) under grant ANR-20-CE23-0002.