Sonia Cafieri
email: sonia.cafieri@enac.fr

Pierre Hansen
email: pierre.hansen@gerad.ca

Leo Liberti
email: liberti@lix.polytechnique.fr

A Locally Optimal Heuristic for Modularity Maximization of Networks

)] and with the hierarchical agglomerative heuristic of Clauset, Newman and Moore [Phys. Rev. E 70, 066111 (2004)]. Computational results are given for a series of problems of the literature with up to 4941 vertices and 6594 edges. They show that the proposed divisive heuristic gives better results than the divisive heuristic of Newman and than the agglomerative heuristic of Clauset et al.

Résumé

L'identification de communautés dans les réseaux se fait actuellement par des heuristiques hiérarchiques agglomératives ou divisives, ansi que par des heuristiques de partitionnement, des hybrides et, dans quelques articles, par des algorithmes exacts. Nous considérons ici le cas de réseaux hiérarchiques dans lequel les communautés doivent être détectées et nous proposons une heuristique divisive qui est localement optimale, dans le sens de ce que chacune des bipartitions successives se fait de manière optimale. Cette heuristique est comparée à l'heuristique hiérarchique divisive spectrale de Newman [Proceedings of the National Academy of Sciences, USA 103, 8577 (2006)] et à l'heuristique hiérarchique agglomerative de Clauset et al. [Phys. Rev. E 70, 066111 (2004)]. On donne des résultats de calcul sur une série de problèmes de la litérature. Ils montrent que l'heuristique divisive proposée donne de meilleurs résultats que la précédente heuristique divisive et que l'heuristique agglomérative cités.

Introduction

Networks, or graphs, are a powerful and versatile tool for the study of complex systems, with many applications in computer science, engineering, transportation, sociology, political science, biology, chemistry and other fields. A network consists of a set of vertices (or nodes) and a set of edges (or lines). Vertices are represented by points and associated with entities, such as customers, users of the World Wide Web, employees in an organization, transmitters, road-crossings, railway stations and atoms. Edges are pairs of vertices and represented by a line joining them. The shape of this line is irrelevant; only its presence or absence matters. Edges represent relationships between the entities associated with the vertices: communication, collaboration, existence of a connection such as a road or a railway line and chemical bonds. A detailed introduction to networks has recently been given by Newman [START_REF] Newman | Networks: an introduction[END_REF].

A very important and much studied problem in network science and its applications is the detection of communities (also called modules or clusters). These are sets of entities, or vertices, which are likely to have some common function. Usually, the number of inner edges, i.e., edges joining two vertices of the same community, is larger than the number of outer edges, i.e., edges joining two vertices of different communities. An in-depth survey of the problem of community detection in graphs has recently been given by Fortunato [START_REF] Fortunato | Community detection in graphs[END_REF]. There are several precise definitions of communities, and corresponding criteria, and many more heuristics as well as a few exact algorithms to find partitions or sets of nested partitions into communities. A heuristic finds a near optimal partition (or sometimes an optimal partition but without proof of its optimality) in moderate time. An exact algorithm finds an optimal partition, with proof of its optimality, hopefully in reasonable time. The most used definition for the quality of a community or of a partition into communities is that of modularity, proposed by Newman and Girvan [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. Modularity of a community is defined as the difference between the number of edges it contains and the expected number of edges that it would contain if all edges were drawn at random, keeping the same distribution of degrees. The modularity of a partition is the sum of the modularities of its communities. See e.g. [START_REF] Fortunato | Community detection in graphs[END_REF][START_REF] Fortunato | Resolution limit in community detection[END_REF][START_REF] Cafieri | Loops and multiple edges in modularity maximization of networks[END_REF] for a discussion of the strengths and weaknesses of the modularity function. Given a network and a partition, modularity can be viewed as a measure of the extent to which the classes of the partition can be considered to be communities. Alternatively, given a network, modularity can be maximized to find an optimal partition, together with its number of clusters and their modularities.

As traditional in Data Analysis, given a set of n entities, clustering heuristics are either hierarchical, i.e., they aim at finding a set of nested partitions, or partitioning schemes, i.e., they aim at finding a single partition (or possibly several partitions into given numbers of clusters). In turn, hierarchical heuristics are divided into agglomerative and divisive ones. Hierarchical agglomerative heuristics [START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF][START_REF] Clauset | Finding community structure in very large networks[END_REF][START_REF] Danon | The effect of size heterogeneity on community identification in complex networks[END_REF][START_REF] Wakita | Finding community structure in mega-scale social networks[END_REF][START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] proceed from an initial partition with n communities each containing a single entity and iteratively merge the pair of entities for which this operation increases most the objective function (e.g., modularity), until all entities belong to the same community. Thus, they find 2n-1 communities which are pairwise disjoint or included one into the other. Hierarchical divisive heuristics [START_REF] Newman | Modularity and community structure in networks[END_REF] proceed from an initial partition containing all entities and iteratively divide a community into two in such a way that the increase in the objective function value (e.g. modularity) is the largest possible, or the decrease in the objective value is the smallest possible. Bipartitions are iterated until a partition into n communities having each a single entity is obtained. Thus, once more, 2n -1 communities are obtained, which are pairwise disjoint or included one into the other. Note that for some objectives, including modularity, mergings or bipartitions can be ended once they do not improve the objective function value any more.

The partitioning and hybrid heuristics rely upon simulated annealing [START_REF] Guimerà | Functional cartography of complex metabolic networks[END_REF][START_REF] Massen | Identifying communities within energy landscapes[END_REF][START_REF] Medus | Detection of community structures in networks via global optimization[END_REF], mean field annealing [START_REF] Lehmann | Deterministic modularity optimization[END_REF], genetic search [START_REF] Tasgin | Community detection in complex networks using genetic algorithms[END_REF], extremal optimization [START_REF] Duch | Community identification using extremal optimization[END_REF], linear programming followed by randomized rounding [START_REF] Agarwal | Modularity-maximizing graph communities via mathematical programming[END_REF], dynamical clustering [START_REF] Boccaletti | Detecting complex network modularity by dynamical clustering[END_REF], multilevel partitioning [START_REF] Djidjev | A scalable multilevel algorithm for graph clustering and community structure detection[END_REF], contraction-dilation [START_REF] Mei | Revealing network communities through modularity maximization by a contraction-dilation method[END_REF], multistep greedy search [START_REF] Schuetz | Efficient modularity optimization by multistep algorithm and vertex mover refinement[END_REF], quantum mechanics [START_REF] Niu | Detecting the community structure in complex networks based on quantum mechanics[END_REF] and many more sources of inspiration [START_REF] Chen | A fast and efficient heuristic algorithm for detecting community structures in complex networks[END_REF][START_REF] Sun | Improved community structure detection using a modified fine-tuning strategy[END_REF][START_REF] Ruan | Identifying network communities with a high resolution[END_REF]17,[START_REF] Kumpula | Limited resolution and multiresolution methods in complex network community detection[END_REF].

Hierarchical heuristics are in principle devised for finding a hierarchy of partitions implicit in the given network when it corresponds to some situation where hierarchy is observed or postulated. Such situations include the description of hierarchies in social organization and networks describing evolutionary processes. Results are presented on a dendrogram which displays visually mergings or divisions of communities together with the values of a characteristic of each community (a variant of the dendrogram, called espalier, allows displaying simultaneously two characteristics of the communities [START_REF] Hansen | Espaliers: A generalization of dendrograms[END_REF]).

The subproblem of choosing at each iteration which pair of communities should be merged is easy. It suffices to consider all O(n 2 ) merging of pairs of entities and compute each time the objective function value for the new community. Moreover, a careful use of data structures often reduces complexity. Low order polynomial hierarchical agglomerative heuristics have been obtained in several classical papers. These include O(n 2 ) algorithms for single-linkage [START_REF] Gower | Minimum spanning trees and single linkage cluster analysis[END_REF], complete linkage and minimum variance [START_REF] Benzécri | Construction d'une classification ascendante hierarchique par la recherche en chaîne des voisins reciproques[END_REF] and several others have a complexity of O(n 2 log n) [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF]. However, in a divisive hierarchical heuristic, the subproblem of finding a bipartition locally optimizing the adopted criterion is more difficult. For some criteria there exists a polynomial algorithm for bipartitioning. For instance, this is the case for the minimum diameter criterion for which there is a O(n 2 ) algorithm. With careful use of data structures, this gives a O(n 2 log n) locally optimal algorithm for hierarchical divisive clustering with the minimum diameter criterion [START_REF] Guénoche | Efficient algorithms for divisive hierarchical clustering with the diameter criterion[END_REF]. The situation is less favorable for the maximum modularity criterion. Indeed, this problem is NP-hard even in the case of two clusters [START_REF] Brandes | On modularity clustering[END_REF]. Nevertheless, as shown below, a non-polynomial algorithm can solve instances with up to 4941 vertices.

We consider only networks with unweighted and undirected edges in the present paper. Its purpose is to propose a locally optimal divisive heuristic for the maximum modularity criterion. To that effect, in the next section the bipartition subproblem is expressed as a quadratic mixed-integer program with a convex relaxation. This problem can then be solved by the CPLEX solver [START_REF] Ilog | ILOG CPLEX 11.0 User's Manual[END_REF]. In Sect. 3 the full hierarchical divisive algorithm is described as well as the previous spectral hierarchical divisive heuristic of Newman [START_REF] Newman | Modularity and community structure in networks[END_REF] and the fast hierarchical agglomerative heuristic of Clauset, Newman and Moore [START_REF] Clauset | Finding community structure in very large networks[END_REF]. A computational comparison of the three heuristics, detailing also the respective contributions of Newman's spectral results and the Kernighan-Lin heuristic, is given in Sect. 4. Conclusions are drawn in Sect. 5.

An exact algorithm for bipartition

We present in this section an exact algorithm for bipartition with maximization of modularity. We model this bipartitioning problem using binary variables to identify to which community each vertex and each edge belongs. In this respect, our model is similar to that of Xu et al. [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF]. These authors proposed in 2007 a modularity maximization model to obtain a partition (generally with more than two communities) of a network. This model is expressed as a mixed integer convex quadratic program. Xu et al. were able to solve exactly instances with up to 104 vertices.

Let G = (V, E) be a graph, or network, with set of vertices V of order n = |V | and set of edges E of size m = |E|. We next recall two equivalent definitions of modularity Q. In the first one, it is expressed as a sum over communities of their modularities [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]:

Q = s [a s -e s ],
where a s is the fraction of all edges that lie within community s and e s is the expected value of the same quantity in a graph in which the vertices have the same degrees but edges are placed at random. In the second one, modularity Q is expressed as a function, for each community, of its number of inner edges and of the sum of degrees of its vertices:

Q = s m s m - d s 2m 2 , (1)
where m s denotes the number of edges in community s, i.e., the subgraph G s = (V s , E s ) with set of vertices V s ⊂ V and set of edges E s having both vertices in V s , and d s denotes the sum of degrees k i of the vertices of community s. Since we aim to find a bipartition, only two sub-modules of the original community have to be considered, i.e., s ∈ {1, 2}. We can express the sum of degrees d 2 of vertices belonging to the second community as a function of the sum of degrees d 1 of vertices belonging to the first one:

d 2 = d t -d 1 , (2)
where d t is the sum of degrees in the community to be bipartitioned and it is equal to 2m at the outset. We rewrite (1) for s ∈ {1, 2}, using (2):

Q = m 1 + m 2 m - d 2 1 4m 2 - d 2 2 4m 2 = = m 1 + m 2 m - d 2 1 4m 2 - d 2 t + d 2 1 -2d t d 1 4m 2 = = m 1 + m 2 m - d 2 1 2m 2 - d 2 t 4m 2 + d t d 1 2m 2 .
(3)

We then introduce binary variables X r1 , X r2 and Y i1 to model the assignment of vertices and edges to the two communities of the bipartition. These variables are defined as follows:

X rs = 1 if edge r belongs to community s 0 otherwise (4)
for r = 1, 2, . . . m and s = 1, 2 and

Y i1 = 1 if vertex i belongs to community 1 0 otherwise, i.e. if vertex i belongs to community 2 (5)
for i = 1, 2, . . . n. Two sets of variables X r1 and X r2 are needed as an edge may belong to the first community, or to the second one, or be a bridge between both of them. One set of variables Y i1 suffices as any vertex which does not belong to the first community must belong to the second.

Moreover, we impose for consistency that any edge r = {v i , v j } with end vertices indiced by i and j can only belong to community s if both of its end vertices belong also to that community:

X r1 ≤ Y i1 ∀r = {v i , v j } ∈ E X r1 ≤ Y j1 ∀r = {v i , v j } ∈ E (6)
and

X r2 ≤ 1 -Y i1 ∀r = {v i , v j } ∈ E X r2 ≤ 1 -Y j1 ∀r = {v i , v j } ∈ E. (7)
Furthermore, we exploit the following expressions in terms of variables X and Y for the number of edges of each of the two communities and the sum of vertex degrees of the first one:

m s = r X rs ∀s ∈ {1, 2}, (8)
d 1 = i∈V1 k i Y i1 . (9)
The sum of vertex degrees of the first community only is needed, because of expression (2).

Maximizing modularity (3) subject to constraints ( 6)-( 7) and ( 8)-( 9) gives a quadratic convex mixedinteger program that can be solved by CPLEX [START_REF] Ilog | ILOG CPLEX 11.0 User's Manual[END_REF]. Indeed, this model contains a single nonlinear but concave term, i.e., -d 2 1 /2m 2 , in the objective function, which is to be maximized. Hence, its continuous relaxation, obtained by removing the integrality constraints on the variables, is a convex quadratic program and easy to solve.

Note that in Xu et al.'s [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] model a number of other constraints are imposed. For instance, constraints are used to express that community s can be nonempty only if community s -1 is so, and lower and upper bounds on the cardinality of the modules are added as an option. Furthermore, symmetry-breaking constraints avoid the computation of equivalent alternative optimal solutions.

Heuristics

In this section, we first recall the hierarchical agglomerative heuristic of Clauset, Newman and Moore [START_REF] Clauset | Finding community structure in very large networks[END_REF] (CNM), and the hierarchical divisive spectral heuristic of Newman [START_REF] Newman | Modularity and community structure in networks[END_REF] (divisive spectral), to which we compare the heuristic of the present paper. We then describe the new heuristic itself.

The CNM heuristic fits into the general scheme for hierarchical agglomerative heuristics of cluster analysis [START_REF] Murtagh | A survey of recent advances in hierarchical clustering algorithms[END_REF]. It can therefore be implemented with a complexity of O(n 2 log n) in worst case. However, a careful exploitation of sparsity of the graph under study reduces its worst-case complexity to O(m log n) and its complexity in practice to O(n log 2 n), which is close to linear time. Results in terms of partitions and the dendrograms obtained are the same for both implementations. The CNM heuristic proceeds from an initial partition with n clusters each containing a single entity. Then, it merges iteratively the two clusters in the current partition for which the modularity increases the most. The formula

∆Q ij = 1/2m -k i k j /(2m) 2 if v i , v j are connected 0 otherwise
is used initially, then the ∆Q ij for all edges are updated each in constant time. Mergings take place as long as the best of them increases modularity (or, in other words, there is a positive ∆Q ij ).

Two remarks are in order. First, in hierarchical agglomerative heuristics, errors, i.e., assignment of two entities to the same cluster while they should be in different clusters in the (or all) optimal partition(s), are never corrected. Second, as n is usually much larger than the number of clusters in the optimal partition, there are close to n mergings before reaching the best partition obtained by the heuristic and hence many occasions of error.

The hierarchical divisive heuristic of Newman proceeds from an initial partition containing all entities by iteratively splitting one of its clusters, as long as this operation increments modulatity. Two questions have to be answered to specify this heuristic: which cluster should be selected for splitting at each iteration and how should the splitting be made. The answer to the former question is unimportant as the best partition obtained does not depend on the order of the splittings but only on the way they are performed. The second question is difficult: indeed, finding the optimal, i.e., modularity maximizing, splitting of a cluster contains the problem of maximum modularity bipartition, which, as shown by Brandes et al. [START_REF] Brandes | On modularity clustering[END_REF], is NP-hard.

For large instances, splitting will have to be done in a heuristic way. Newman [START_REF] Newman | Modularity and community structure in networks[END_REF] proposes two ways to do so. The first is based on spectral graph theory. The first eigenvector of the modularity matrix

B = (b ij ) with b ij = a ij -k i k j /2m
is computed. The entities corresponding to positive components of this eigenvector form one community and the remaining ones the other.

Results of splitting according to the first eigenvector can be improved by a variety of heuristics. Newman suggests the use of the Kernighan-Lin heuristic [START_REF] Kernighan | An efficient heuristic procedure for partitioning graphs[END_REF]. This heuristic proceeds from an initial bipartition to a sequence of re-assignments of one entity from a community to the other. At each step the re-assignment which improves most, or deteriorates least, the objective function value, is selected, performed and further re-assignments of the moved entity forbidden. Once no more re-assignments are allowed, the best partition found among the n partitions considered in the sequence is selected as new initial partition. The whole procedure stops when a full sequence of n re-assignments does not lead to any improvement.

Again two remarks are in order. First, errors, i.e., assigning two entities to different communities when both belong to the same community in the optimal partition are never corrected. Second, as the number of communities in the optimal partition tends to be small, few splittings will take place and thus there are few occasions for errors to be made.

The heuristic of the present paper proceeds along the same lines as the hierarchical divisive heuristic of Newman. The difference is that the exact bipartitioning method of Sect. 2 is used for the splitting step. Due to the difficulty of this subproblem and present limitations of nonlinear integer programming, the exact bipartitioning could be applied only to small or medium-size networks with up to 4941 vertices. The proposed heuristic is locally optimal in that the splitting step is done optimally, but not globally optimal in that the new heuristic is a greedy one, i.e., each splitting step is done without considering the consequences on further steps, but better results might be obtained if this was done. Note that there are some cases where a greedy heuristic is optimal, the best known one in the field of cluster analysis being the single-linkage algorithm which maximizes the split (or minimum dissimilarity between pairs of entities in different communities) of the partitions obtained at all levels [START_REF] Gower | Minimum spanning trees and single linkage cluster analysis[END_REF][START_REF] Delattre | Bicriterion cluster analysis[END_REF]. This is not the case for the proposed hierarchical divisive heuristic, as shown by an example in Sect. 4.

Computational results

In all experiments, we considered a set of 12 well known problems from the literature. They are listed together with their order n and size m in the first three columns of Table 1. They are all undirected and unweighted networks without loops. All data are available on sites listed in [START_REF] Noack | Multi-level algorithms for modularity clustering[END_REF]. Optimal solutions of these test problems are given in [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF]. The bipartition subproblem is solved using CPLEX [START_REF] Ilog | ILOG CPLEX 11.0 User's Manual[END_REF]. Other solvers for convex mixed-integer quadratic programs did not perform as well in preliminary tests.

In a first series of experiments, we compared the maximum modularity values and the corresponding number of clusters for the hierarchical agglomerative heuristic of Clauset et al. [START_REF] Clauset | Finding community structure in very large networks[END_REF], and for the locally optimal divisive heuristic of the present paper. We also list maximum modularity and number of clusters obtained by an exact algorithm of [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF], when possible, i.e., for datasets from 1 to 11. The modularity obtained by the proposed divisive heuristic for the 12th dataset is 0.9394. This appears to be the best value currently known for this dataset. Indeed, using the best heuristic from the extensive comparison by Noack and Rotta [START_REF] Noack | Multi-level algorithms for modularity clustering[END_REF], gives a value of 0.93854.

In the following, we refer to the locally optimal divisive heuristic of the present paper as divisive CHL.

It appears that:

• the locally optimal hierarchical divisive heuristic CHL always gives partitions with smaller modularity than the exact algorithm. The difference, however, is moderate and goes from 0.18018% to 2.36603% of the optimal value. The average error is 0.82540%. • The agglomerative heuristic always gives partitions with smaller modularity than the divisive one and, by transitivity, than the exact one. The difference this time is much more substantial and goes from 1.21494% to 12.9848%, relative to the exact solution. The average error is 5.52342%. The error between the divisive and the agglomerative heuristics goes from 0.780214% to 10.9013% of the value obtained by the divisive heuristic. The average error is 4.75179%, which is again substantial. • The divisive heuristic CHL gives a partition with fewer communities than the exact one in 3 cases out of 11, the same number in 3 cases and a larger number in the remaining 5 cases. The average number of communities obtained with the divisive heuristic is 9.27273 vs. 8.81818 for the exact method, so it is slightly in excess. • The agglomerative heuristic gives a partition with fewer communities than the exact one in 7 cases out of 11, a partition with the same number of communities than the exact one in 2 cases out of 11, and a partition with more communities in the remaining 2 cases. The average number of communities obtained with the agglomerative heuristic is 8 vs. 8.81818 for the exact method, so somewhat smaller.

In a second series of experiments we compared the divisive heuristic CHL of the present paper with the previous divisive heuristic of Newman [38]. To better understand the performance of Newman's heuristic, we consider three versions of it. In the first one, the bipartition at each iteration is done solely on the basis of the leading eigenvector of the modularity matrix. In the second one, the Kernighan-Lin heuristic is applied at each iteration to a randomly generated solution. In the third one, the Kernighan-Lin heuristic is applied at each iteration to the solution given by the leading eigenvector. Results are presented in Tables 2, 3, 4. It appears that Table 1: Comparison of results of Clauset et al.'s heuristic (CNM), the proposed locally optimal divisive heuristic (divisive CHL), and an exact algorithm for modularity maximization [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF] (exact) on real world datasets. M denotes the number of communities and Q the modularity value of the best solution found. error(%) denotes the percentage error for the two heuristics with respect to the exact algorithm. Average values are given for the first 11 datasets as the 12-th one cannot be solved exactly in reasonable time. n and m are the number of vertices and the number of edges of the networks. We consider Zachary's karate club dataset [START_REF] Zachary | An information flow model for conflict and fission in small group[END_REF] describing friendship relationships between members of a club, Lusseau's dolphins dataset [START_REF] Lusseau | The bottlenose dolphin community of doubtful sound features a large proportion of long-lasting associations. can geographic isolation explain this unique trait?[END_REF] describing communications between dolphins in Doubtful Sound New Zealand, Hugo's Les Misérables dataset describing characters in Victor Hugo's masterpiece and their interactions, compiled by Knuth [START_REF] Knuth | The Stanford GraphBase: A Platform for Combinatorial Computing[END_REF], a dataset main) on classes and relationships from a software project related to Graph Drawing [22], a network dealing with protein interactions [START_REF] Dartnell | Self-similar community structure in a network of human interactions[END_REF], Krebs' political books dataset [START_REF] Krebs | [END_REF], a dataset representing the schedule of football games between American college teams [START_REF] Girvan | Community structure in social and biological networks[END_REF], another dataset on classes and relationships from a software project [22], a network dealing with connections between US airports [46], a dataset on a coauthorship network of scientists working on network theory and experiment, compiled by M. Newman [START_REF] Newman | Finding community structure in networks using the eigenvectors of matrices[END_REF], a network describing electronic circuits [52] and a network representing the topology of the Western States Power Grid of the United States [START_REF] Watts | Collective dynamics of 'small-world' networks[END_REF]. • the first version always gives partitions with substantially smaller modularities than the divisive heuristic of the present paper or the exact algorithm. The error relative to the partition given by the exact algorithm goes from 5.81428% to 18.5189%. The average error is 10.7419%. The error relative to the partition given by the divisive heuristic of the present paper goes from 5.63859% to 18.0227%. The average error is 9.99712%. • The second version gives results even worse than those of the first version. The error relative to the partition given by the exact algorithm goes from 5.95934% to 28.6442%. The average error is 14.5856%.

dataset n m agglomerative CNM divisive CHL exact M Q error(%) M Q error(%) M Q

• The third version gives better results

The error relative to the partition given by the exact algorithm goes from 0.09863% to 6.07995%. The average error is 3.02627%. The error relative to the partition given by the divisive heuristic of the present paper goes from -0.081704% to 5.46305%. The average error is 2.21592%. Observe that in the case of the political books instance, the modularity obtained with the divisive spectral + KL heuristic, i.e., 0.52672, is slightly better than that one obtained with the heuristic of the present paper, i.e, 0.52629. This is not a numerical error but illustrates, as mentioned above, that the proposed heuristic is locally but not globally optimal. The observation that version 3 is better than version 1 and that version 1 is better than version 2 was already made by Newman [START_REF] Newman | Modularity and community structure in networks[END_REF]. • Computing times of the heuristic proposed in this paper and an exact column generation algorithm for modularity maximization of [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF] are given in Table 5. As the computers used are not the same for the two cases, these results should only be considered as indicative. Ratios between these computing times are given in the penultimate column and similar ratios, corrected by dividing those times by the clock frequency of the computer used, in the last column. The time of the heuristic CHL is less than the time of the exact algorithm in 9 cases over 12 and in several cases very substancially so. Moreover, the heuristic could solve the last problem, which is about 10 times larger than the penultimate, in reasonable time, while the exact algorithm could not.

• Currently, it is not possible to solve very large instances with the heuristic of the present paper or large ones with the exact column generation algorithm for modularity maximization of [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF]. One cannot therefore be sure that the differences in performances between the heuristics described above and/or the exact algorithm extend to the resolution of large instances, except in a few cases. A comparison between the agglomerative heuristic CNM and the complete version of the divisive heuristic of Newman is part of Table 1 of [START_REF] Newman | Modularity and community structure in networks[END_REF]. The corresponding columns are reproduced in Table 6, together with the percentage errors in the modularities obtained. This confirms that the agglomerative heuristic of Clauset et al.

gives poor results for small as well as for large instances.

Table 2: Comparison of results of Newman's spectral divisive heuristic (divisive spectral), the proposed locally optimal divisive heuristic (divisive CHL), and an exact algorithm for modularity maximization [START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF] (exact) on real world datasets. M denotes the number of communities and Q the modularity value of the best found solution. error(%) denotes the percentage error for the two heuristic with respect to the exact algorithm. Average values are given for the first 11 datasets as the 12-th one cannot be solved exactly in reasonable time. n and m are the number of vertices and the number of edges of the networks.

dataset n m divisive spectral divisive CHL exact M Q error(%) M Q error(%) M Q

Conclusions

In this paper we presented a hierarchical divisive heuristic for modularity maximization which is locally optimal, i.e., such that the bipartition obtained at each iteration is guaranteed to be optimal. This heuristic can be used for two purposes: approximate modularity maximization of general networks or modularity maximization of networks which are known to correspond to some hierarchy either natural or man made. In the former case, experimental results show that the partitions obtained with the proposed divisive heuristic tend to have a modularity value close to that of optimal partitions (recall that the average error observed in our experiments is 0.82540%). So, the partition found can be considered as a fairly good approximation Table 6: Comparison of results [START_REF] Newman | Modularity and community structure in networks[END_REF] of the Clauset et al.'s agglomerative heuristic (agglomerative CNM) and Newman's spectral divisive heuristic with the Kernighan-Lin refinement (divisive spectral + KL). Q denotes the modularity value of the best found solution. error(%) denotes the percentage error for the first heuristic with respect to the second one. n is the number of vertices. Jazz dataset describes musicians which worked together [47], metabolic dataset describes chemical reactions as well as the regulatory interactions that guide these reactions in C. elegans [3], e-mail dataset describes e-mail interchanges between members of a university [START_REF] Guimerà | Self-similar community structure in a network of human interactions[END_REF], key signing (PGP) datadset describes the giant component of the network of users of the Pretty-Good-Privacy algorithm for secure information interchange [7], physicists (cond-mat) dataset describes a collaboration network of scientists posting preprints on the condensed matter archive [43]. of the optimal one, or at least as a tentative solution to be improved upon by various local improvement heuristics.

dataset
Our experiments also show that modularity maximization with the hierarchical agglomerative heuristic of Clauset, Newman and Moore tends to have a much larger error (with an average error of 5.52342% for our experiments).

It is for the latter case that our heuristic is tailored. We therefore compared it with three versions of the divisive heuristic of Newman. As observed by that author, a two-phase method in which a first bipartition is made in a splitting step on the basis of the first eigenvector of the modularity matrix followed by application of the Kernighan-Lin heuristic gives the best results. Still, the proposed heuristic appears to be even better, indeed the average error relative to the exact solution is reduced more than threefold, i.e., from 3.02627% to 0.82540%.

Table 3 :

 3 Comparison of results of the divisive Kernighan-Lin based heuristic (KL), the proposed locally optimal divisive heuristic (divisive CHL), and an exact algorithm for modularity maximization[START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF] (exact) on real world datasets. M denotes the number of communities and Q the modularity value of the best found solution. error(%) denotes the percentage error for the two heuristics with respect to the exact algorithm. Average values are given for the first 11 datasets as the 12-th one cannot be solved exactly in reasonable time. n and m are the number of vertices and the number of edges of the networks. Note that values of Q and error percentage are given with less digits that elsewhere in this paper due to the fact that the Kernighan-Lin heuristic depends upon the random partition to which it is applied.

	dataset	n	m		KL			divisive CHL		exact
				M	Q	error(%)	M	Q	error(%)	M	Q
	karate	34	78	2	0.372	11.43	4	0.41880	0.23583	4	0.41979
	dolphin	62	159	4	0.477	9.745	4	0.52646	0.38977	5	0.52852
	les miserables	77	254	3	0.489	12.73	8	0.54676	2.36603	6	0.56001
	A00 main	83	135	3	0.450	15.31	7	0.52806	0.53494	9	0.53090
	p53 protein	104	226	11	0.402	24.83	7	0.52843	1.25203	7	0.53513
	political books	105	441	4	0.496	5.959	4	0.52629	0.18018	5	0.52724
	football	115	613	5	0.538	11.048	10	0.60091	0.60539	10	0.60457
	A01 main	249	635	6	0.512	19.023	15	0.62877	0.65255	14	0.63290
	usair97	332	2126	6	0.339	7.992	8	0.35959	2.33840	6	0.36820
	netscience main	379	914	8	0.606	28.644	20	0.84702	0.18619	19	0.84860
	s838	512	819	5	0.707	13.726	15	0.81663	0.33805	12	0.81940
	power	4941	6594	8	0.649	-	40	0.93937	-	-	-
	average			5.2	0.48970	14.5856	9.3	0.57525	0.82540	8.8	0.57957

Table 4 :

 4 Comparison of results of Newman's spectral divisive heuristic with the Kernighan-Lin refinement (divisive spectral + KL), the proposed locally optimal divisive heuristic (divisive CHL), and an exact algorithm for modularity maximization[START_REF] Aloise | Column generation algorithms for exact modularity maximization in networks[END_REF] (exact) on real world datasets. M denotes the number of communities and Q the modularity value of the best found solution. error(%) denotes the percentage error for the two heuristics divisive spectral + KL and divisive with respect to the exact algorithm. error div(%) denotes the percentage error for the divisive spectral + KL with respect to the proposed locally optimal heuristic. Average values are given for the first 11 datasets as the 12-th one cannot be solved exactly in reasonable time. n and m are the number of vertices and the number of edges of the networks. Note that values of Q and error percentage are given with less digits that elsewhere in this paper due to the fact that the Kernighan-Lin heuristic depends upon the random partition to which it is applied.

	dataset	n	m		divisive spectral + KL			divisive CHL		exact
				M	Q	error div(%) error(%) M	Q	error(%) M	Q
	karate	34	78	4	0.419	0	0.236	4	0.41880	0.23583	4	0.41979
	dolphin	62	159	5	0.508	3.415	3.792	4	0.52646	0.38977	5	0.52852
	les miserables	77	254	7	0.538	1.533	3.862	8	0.54676	2.36603	6	0.56001
	A00 main	83	135	7	0.527	0.199	0.733	7	0.52806	0.53494	9	0.53090
	p53 protein	104	226	6	0.518	1.930	3.158	7	0.52843	1.25203	7	0.53513
	political books	105	441	4	0.527	-0.081	0.099	4	0.52629	0.18018	5	0.52724
	football	115	613	8	0.579	3.638	4.221	10 0.60091	0.60539	10 0.60457
	A01 main	249	635	16	0.594	5.463	6.080	15 0.62877	0.65255	14 0.63290
	usair97	332	2126	7	0.358	0.501	2.827	8	0.35959	2.33840	6	0.36820
	netscience main	379	914	23	0.820	3.191	3.371	20 0.84702	0.18619	19 0.84860
	s838	512	819	13	0.779	4.587	4.910	15 0.81663	0.33805	12 0.81940
	power	4941 6594	8	0.791	-		40 0.93937	-	-	-
	average			9.09 0.560731	2.21592	3.02627	9.3 0.57525	0.82540	8.8 0.57957

Table 5 :

 5 Comparison of times for the proposed locally optimal heuristic and an exact algorithm for modularity maximization. All results are in seconds of CPU. Solutions given by the proposed heuristic were obtained on a 2.4 GHz Intel Xeon CPU of a computer with 8GB RAM shared by three other similar CPU running Linux. Solutions given by the exact algorithm were obtained on a dual processor computer Intel Pentium with 3.20 GHz, 3GB RAM running Linux. In the penultimate column, ratios between computing times of the proposed heuristic and the exact algorithm are reported. In the last column, ratios between these computing times are corrected by dividing CPU times by the clock frequency of the computer used.

	dataset	n	m	time divisive CHL	time exact	CHL/exact	(CHL/exact)'
	karate	34	78	0.42	0.34	1.2353	1.6471
	dolphin	62	159	1.40	7.75	0.1806	0.2409
	les miserables	77	254	4.52	7.26	0.6226	0.8301
	A00 main	83	135	0.89	3.66	0.2432	0.3242
	p53 protein	104	226	7.81	11.60	0.6733	0.8977
	political books	105	441	12.09	45.65	0.2648	0.3531
	football	115	613	338.19	249.41	1.3559	1.8979
	A01 main	249	635	656.54	1014.48	0.6472	0.8629
	usair97	332	2126	33157.85	16216.77	2.0447	2.7262
	netscience main	379	914	22.84	1615.14	0.0141	0.0188
	s838	512	819	38.37	7655.56	0.0050	0.0067
	power	4941	6594	4498.21	-	-	-


Daniel Aloise
email: daniel.aloise@gerad.ca

Sonia Cafieri
email: sonia.cafieri@enac.fr

Gilles Caporossi
email: cgilles.caporossi@gerad.ca

Pierre Hansen
email: pierre.hansen@gerad.ca

Leo Liberti
email: dliberti@lix.polytechnique.fr

Sylvain Perron
email: sylvain.perron@gerad.ca

Column Generation Algorithms for Exact Modularity Maximization in Networks

Finding modules, or clusters, in networks currently attracts much attention in several domains. The most studied criterion for doing so, due to Newman and Girvan [Physical Review E 69, 026113 (2004)], is modularity maximization. Many heuristics have been proposed for maximizing modularity and yield rapidly near optimal solution or sometimes optimal ones but without a guarantee of optimality. There are few exact algorithms, prominent among which is a paper by Xu et al. [Eur. Physical Journal B 60, 231-239 (2007)]. Modularity maximization can also be expressed as a clique partitioning problem and the row generation algorithm of Grötschel and Wakabayashi [Mathematical Programming 45, 59-96 (1989)] applied. We propose to extend both of these algorithms using the powerful column generation methods for linear and non linear integer programming. Performance of the four resulting algorithms are compared on problems from the literature. Instances with up to 512 entities are solved exactly for the first time. Moreover, the computing time of previously solved problems are reduced substantially.

Résumé

L'identification de modules, ou groupes, dans les réseaux attire actuellement beaucoup d'attention dans plusieurs domaines. Pour arriver à identifier les modules, le critère le plus utilisé jusqu'à présent est celui de la maximisation de la modularité, critère introduit par Newman et Girvan [Physical Review E 69, 026113 (2004)]. Plusieurs méthodes heuristiques ont été proposées pour ce critère et procurent rapidement une solution quasi optimale ou parfois une solution optimale sans toutefois une preuve de son optimalité. Parmi le très faible nombre de méthodes exactes proposées, la plus importante est très certainement celle proposée par Xu et al. [Eur. Physical Journal B 60, 231-239 (2007)]. La maximisation de la modularité peut également être exprimée comme un problème de partionnement en cliques pouvant être résolu par l'algorithme de génération de contraintes de Grötschel et Wakabayashi [Mathematical Programming 45, 59-96 (1989)]. Nous proposons d'étendre ces deux algorithmes en utilisant les techniques avancées de la génération de colonnes pour la programmation linéaire et non linéaire entière. La performance des quatre algorithmes est comparée sur des problèmes provenant de la littérature. Des instances contenant jusqu'à 512 entités sont résolues exactement pour la première fois. De plus, les temps de résolution pour des problèmes déjà résolus sont réduits substantiellement.

Introduction

Clustering is an important chapter of data analysis and data mining with numerous applications in a variety of fields. It aims at solving the following general problem: given a set of entities, find subsets, or clusters, which are homogeneous and/or well-separated. As the concepts of homogeneity and of separation can be made precise in many ways, there are a large variety of clustering problems [START_REF] Hansen | Cluster analysis and mathematical programming[END_REF][START_REF] Jain | Data clustering: A review[END_REF][START_REF] Kaufman | Finding Groups in Data. An Introduction to Cluster Analysis[END_REF][START_REF] Mirkin | Clustering for Data Mining: A Data Recovery Approach[END_REF]. These problems in turn are solved by exact algorithms or, more often and particularly for large data sets, by heuristics. An exact algorithm provides, hopefully in reasonable computing time, an optimal solution together with a proof of its optimality. A heuristic provides, usually in moderate computing time, a near optimal solution or sometimes an optimal solution but without proof of its optimality.

In the last decade, clustering on networks has been extensively studied, mostly in the physics and computer science research communities. Rather than using the term cluster, the words module or community are often adopted in the physics literature. Recall that a network, or graph, G = (V, E) is composed of a set V of n vertices and a set E of m edges, which join pairs of vertices. A vertex v j is represented by a point and an edge e ij = {v i , v j } by a line joining its two end vertices v i and v j . The shape of this line does not matter, only the presence or absence of an edge is important. In a simple graph, there is at most one edge between any pair of vertices, otherwise one has a multigraph. With a slight abuse of set notation, a loop e ii = {v i , v i } is an edge for which both end vertices coincide. The degree k i of a vertex v i ∈ V is the number of edges incident with v i . A subgraph G S = (S, E S ) of a graph G = (V, E) induced by a set of vertices S ⊆ V is a graph with vertex set S and edge set E S equal to all edges with both vertices in S. Such a subgraph corresponds to a cluster (or module or community) and many heuristics aim at finding a partition of V into pairwise disjoint nonempty subsets V 1 , V 2 , . . . , V N inducing subgraphs of G. Various objective functions have been proposed for evaluating such a partition. Roughly speaking, one seeks modules which contain more inner edges (with both vertices in the same module) than cut edges (with vertices in different modules). The degree k i of the vertex v i can be split in two: the indegree k in i or number of neighbors within its community and the outdegree k out i or number of neighbors outside its community. Several concepts of community follow.

In 2004, Radicchi et al. [START_REF] Radicchi | Defining and identifying communities in networks[END_REF] defined a community in the strong sense as a subgraph all vertices of which have larger indegree than outdegree and a community in the weak sense as a subgraph for which the sum of vertex indegrees is larger than the sum of vertex outdegrees. As an inner edge contributes by two to the sum of the indegrees and a cut edge contributes by one to the sum of outdegrees, the number of inner edges in a community in the weak sense must be at least as large as half the number of cut edges. As cut edges contribute to the sum of degrees of two communities, this definition entails that for the network as a whole the number of inner edges is larger than the number of cut edges. Recently, several extensions of the definition of the community in the weak sense have been proposed. One may consider the difference for each community of the sum of indegrees and the sum of outdegrees. Then summing these contributions for all communities gives a multiway cut problem. Another approach is to normalize the contribution of each community by dividing it by its number of vertices [START_REF] Li | Quantitative function for community detection[END_REF]. The resulting function, to be maximized, is called modularity density. Alternatively, contributions of communities may be divided by their number of edges [START_REF] Medus | Alternative approach to community detection in networks[END_REF]. Finally, one may consider maximizing, in a divisive hierarchical method, the minimum ratio of the number of edges in a community divided by the number of cut edges [START_REF] Cafieri | Edge ratio and community structure in networks[END_REF].

A different, and currently mainstream, approach was inaugurated by Newman and Girvan also in 2004 [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. They propose to find a partition of V which maximizes the sum, over all modules, of the number of inner edges minus the expected number of such edges assuming that they are drawn at random with the same distribution of degrees as in G. In [START_REF] Newman | Finding and evaluating community structure in networks[END_REF] the following precise definition of modularity is given:

Q = s [a s -e s ] , (1)
where a s is the fraction of all edges that lie within module s and e s is the expected value of the same quantity in a graph in which the vertices have the same expected degrees but edges are placed at random. A maximum value of Q near to 0 indicates that the network considered is close to a random one (barring fluctuations), while a maximum value of Q near to 1 indicates strong community structure. Observe that maximizing modularity gives an optimal partition together with the optimal number of modules.

The modularity maximization problem has been extensively studied both from the algorithmic and from the applications viewpoints. Some papers discuss a few of its mathematical properties, among which are the following.

(i) Even if the networks under study most of the time have no loops nor multiple edges, the expected number of loops may be positive. Moreover, for some pairs of end vertices, their expected number of edges may be greater than 1 [START_REF] Massen | Identifying communities within energy landscapes[END_REF]. (ii) Modularity maximization suffers from a resolution limit [START_REF] Fortunato | Resolution limit in community detection[END_REF], i.e., when the network is large small modules can be absorbed by larger ones even if they are very dense. (iii) The partition with optimal or near optimal modularity may contain modules which are not communities in the weak sense, or even in the most weak sense [START_REF] Zhang | Modularity optimization in community detection of complex networks[END_REF] that twice the number of inner edges is never less than the number of edges joining this module to another one.

Modification to the model and/or to heuristics have been proposed to address these problems. Changing slightly the null model (i.e. the value of e s in Eq. ( 1)), by simulation [START_REF] Massen | Identifying communities within energy landscapes[END_REF] or by using an analytical formula [START_REF] Cafieri | Loops and multiple edges in modularity maximization of networks[END_REF] allows the removal of loops and multiple edges. A parameter can also be introduced together with multiresolution heuristics [START_REF] Kumpula | Limited resolution and multiresolution methods in complex network community detection[END_REF][START_REF] Zhang | A Constrained Optimization Method for Community Detection[END_REF]. Nevertheless, the original modularity maximization appears presently to be much used. In this paper, we limit ourselves to this classical problem.

Brandes et al. [START_REF] Brandes | On modularity clustering[END_REF] have shown that modularity maximization is NP-hard.

Numerous heuristics have been proposed to maximize modularity. They are based on divisive hierarchical clustering, agglomerative hierarchical clustering, partitioning, and hybrids. They rely upon various criteria for agglomeration or division [START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF][START_REF] Clauset | Finding community structure in very large networks[END_REF][START_REF] Danon | The effect of size heterogeneity on community identification in complex networks[END_REF][START_REF] Wakita | Finding community structure in mega-scale social networks[END_REF][START_REF] Blondel | Fast unfolding of communities in large networks[END_REF], simulated annealing [START_REF] Massen | Identifying communities within energy landscapes[END_REF][START_REF] Guimerà | Functional cartography of complex metabolic networks[END_REF][START_REF] Medus | Detection of community structures in networks via global optimization[END_REF], mean field annealing [START_REF] Lehmann | Deterministic modularity optimization[END_REF], genetic search [START_REF] Tasgin | Community detection in complex networks using genetic algorithms[END_REF], extremal optimization [START_REF] Duch | Community identification using extremal optimization[END_REF], spectral clustering [START_REF] Newman | Modularity and community structure in networks[END_REF][START_REF] Richardson | Spectral tripartitioning of networks[END_REF][START_REF] Sun | Improved community structure detection using a modified fine-tuning strategy[END_REF], linear programming followed by randomized rounding [START_REF] Agarwal | Modularity-maximizing graph communities via mathematical programming[END_REF], dynamical clustering [START_REF] Boccaletti | Detecting complex network modularity by dynamical clustering[END_REF], multilevel partitioning [START_REF] Djidjev | A scalable multilevel algorithm for graph clustering and community structure detection[END_REF], contraction-dilation [START_REF] Mei | Revealing network communities through modularity maximization by a contraction-dilation method[END_REF], multistep greedy search [START_REF] Schuetz | Efficient modularity optimization by multistep greedy algorithm and vertex mover refinement[END_REF], quantum mechanics [START_REF] Niu | Detecting the community structure in complex networks based on quantum mechanics[END_REF] and other approaches [START_REF] Kumpula | Limited resolution and multiresolution methods in complex network community detection[END_REF][START_REF] Blondel | Fast unfolding of communities in large networks[END_REF][START_REF] Sun | Improved community structure detection using a modified fine-tuning strategy[END_REF][START_REF] Chen | A fast and efficient heuristic algorithm for detecting community structures in complex networks[END_REF][START_REF] Ruan | Identifying network communities with a high resolution[END_REF][START_REF] Fan | Accuracy and precision of methods for community identification in weighted networks[END_REF].

In contrast, papers proposing exact algorithms or using mathematical programming are rare for modularity maximization. There are two approaches. In the first one, the original graph G = (V, E) is replaced by a complete weighted graph K n = (V, E ′ ) of order n = |V | as G and such that for any pair of vertices v i , v j ∈ G, K n has an edge e ′ ij with a weight equal to the modularity of that edge. In addition to the complexity result mentioned above, Brandes et al. [START_REF] Brandes | On modularity clustering[END_REF] give an integer programming formulation for modularity maximization and mention that the optimal solution of two test problems with 34 and with 105 entities were determined. Their approach is in fact close to the work of Grötschel and Wakabayashi [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF][START_REF] Grötschel | Facets of the clique partitioning polytope[END_REF] on clique partitioning. It is discussed in the next section. The second approach works directly on the original graph G = (V, E). Xu, Tsoka and Papageorgiou [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] propose a mixed integer convex quadratic programming model, discussed below. They solve exactly four test problems with up to 104 entities.

The purpose of the present paper is to assess and advance the state of the art of algorithms for exact modularity maximization. To this effect, we discuss and compare four exact algorithms, two of which are new. Two of these algorithms work on a reduction of modularity maximization to clique partitioning; the other two work on the direct formulation. They are: (i) the row generation algorithm of [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF], which subsumes the algorithm of [START_REF] Brandes | On modularity clustering[END_REF], (ii) a new stabilized column generation algorithm for clique partitioning which enhances the efficiency of that approach, (iii) the mixed integer convex quadratic programming approach of [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF], (iv) another new stabilized column generation algorithm which enhances the efficiency of the second approach.

Column generation algorithms implicitly take into account all possible communities (or in other words all subsets of the set of entities under study). They replace the problem of finding simultaneously all communities in an optimal partition by a sequence of optimization problems for finding one community at a time, or more precisely a community which improves the modularity of the current solution. So, problems are solved much faster than with previous algorithms and larger instances can be tackled, the largest to date having 512 entities.

Clearly, in many applications of modularity maximization, instances are still larger and sometimes very much larger, than those which can be solved exactly with the proposed algorithms. Nevertheless, it is our belief that exact algorithms are worthy of study for several reasons.

(i) Instances which are more than toy problems can presently be solved exactly. As shown below, these include many problems used to illustrate the performance of various heuristics. (ii) As mentioned in several papers, finding the significance of the detected communities is difficult. Indeed, various heuristics often lead to partitions that disagree on two or more of their communities. Having an exact solution solves the problem of separating possible inadequacies of the model from eventual errors resulting from the use of heuristics. Unsuspected communities may be interpreted with more confidence and proposed to the user for a substantive analysis. This is one of the main aims of clustering. (iii) If solving a given instance proves to be too time consuming, the exact algorithm may often be stopped and the best solution found considered as a heuristic one. It is not uncommon that the optimal solution is found at an early stage of the resolution. Then, the problem of maximizing modularity will, in fact, be solved but without a proof of optimality. (iv) An exact algorithm can provide a benchmark of exactly solved instances which can be used to compare heuristics and fine tune them. More precisely, the comparison of the symmetric differences between the optimal solution and the heuristically obtained ones may suggest additional moves which improve the heuristic under study. Iterating this approach with several heuristics may lead to performant hybrids. As a rule of thumb, a sophisticated heuristic should be able to find quickly an optimal solution for most or possibly all practical instances which can be solved exactly with a proof of optimality. (v) Using projection, i.e., fixing some of the communities found by a heuristic, an exact algorithm can be applied to this reduced network in order to improve the heuristic solution [START_REF] Cafieri | Improving heuristics for network modularity maximization using an exact algorithm[END_REF]. (vi) Getting improved heuristics will in turn lead to more efficient exact algorithms. Indeed, there are usually some steps in such algorithms which can be solved heuristically, always or most of the time, without forfeiting the guarantee of optimality. In column generation algorithms, as discussed below, knowledge of a good or possibly optimal initial solution enhances stabilization substantially. Heuristics are also useful in solving the auxiliary problem. (vii) An exact algorithm can also be used as a tool for the theoretical study of maximum modularity partitions, e.g., the fact that they satisfy or not conditions on the number of inner and of cut edges of communities, or conditions of robustness [START_REF] Zhang | Modularity optimization in community detection of complex networks[END_REF]. (viii) With active research, regularly improved mathematical programming packages such as CPLEX [START_REF] Ilog | ILOG CPLEX 11.0 User's Manual[END_REF] and increasing computer power, the size of problems solved exactly is likely to increase substantially with time. A comparison with a central problem of Operations Research, i.e., the traveling salesman problem, illustrates this point. Given a set of cities and pairwise distances between them, the traveling salesman problem is to find a minimum length tour visiting once and only once each city. In 1954, Dantzig, Fulkerson, and Johnson [START_REF] Dantzig | Solution of a large-scale traveling salesman problem[END_REF] solved optimally an instance with 49 cities. In 2009, Applegate et al. [START_REF] Applegate | Certification of an optimal tsp tour through 85,900 cities[END_REF] were able to solve to optimality an instance with 85 900 cities.

2 Modularity maximization as clique partitioning

Row generation

The modularity function was expressed above as a sum of values over all communities. As shown in [START_REF] Newman | Modularity and community structure in networks[END_REF], modularity can also be written as a sum of values over all edges of the complete graph K n .

Q = 1 2m i,j∈V a ij - k i k j 2m δ(c i , c j ) (2)
where m corresponds to the cardinality of E, k i , k j are the degrees of vertices i and j, respectively, a ij is the ij component of the adjacency matrix of G, equal to 1 if vertices i and j are adjacent, and to 0 otherwise. Finally, δ(c i , c j ) is the Kronecker symbol equal to 1 if the communities c i and c j , to which i and j belong, are the same, and to 0 otherwise. The quantity kikj 2m is the expected number of edges between vertices i and j in a null model where edges are placed at random, while the distribution of degrees remains the same.

Introducing binary variables x ij equal to 1 if vertices i and j belong to the same module and 0 otherwise, and setting

w ij = 1 m a ij - k i k j 2m , (3)
modularity maximization can be reformulated as a clique partitioning problem. Since K n is complete, it is a clique and any of its induced subgraphs are cliques also. Partitioning G is thus equivalent to partitioning K n into cliques. The resulting partition is an equivalence relation, i.e., reflexive, symmetric and transitive.

From reflexivity (or the fact that each entity is in the same module as itself) all x ii = 1. Thus the sum of elements on the main diagonal is a constant for all partitions, equal to

-C = - i∈V k i k i 2m . (4)
From symmetry, using x ij = x ji , one can eliminate variables corresponding to values of indices i > j. The model can then be written, as in [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF][START_REF] Grötschel | Facets of the clique partitioning polytope[END_REF],

max i<j∈V w ij x ij -C s.t. x ij + x jk -x ik ≤ 1 for all 1 ≤ i < j < k ≤ n, x ij -x jk + x ik ≤ 1 for all 1 ≤ i < j < k ≤ n, (5)
-

x ij + x jk + x ik ≤ 1 for all 1 ≤ i < j < k ≤ n, x ij ∈ {0, 1} for all 1 ≤ i < j ≤ n.
This model has n(n-1)

2 variables and 3 n 3 = n(n-1)(n-2) 2 = O(n 3 ) constraints.
The first three sets of constraints express transitivity, i.e., if entities i and j are in the same module and entities j and k are in the same module, then entities i and k must be in the same module. The fourth (and last) set of constraints expresses integrality, i.e., edges are present in the solution entirely or not at all. Problem (5) is a linear program in 0-1 variables and thus small instances may be solved by integer linear programming packages such as CPLEX. In the solution process, the linear programming relaxation obtained by replacing the constraints x ij ∈ {0, 1} by x ij ∈ [0, 1] is first solved. If the optimal solution of this relaxation is in integers, which is often the case, it corresponds to a partition of maximum modularity. Should the solution of the continuous relaxation be fractional, one can branch, or add one or several cutting planes, i.e., additional linear constraints which cut off the current fractional solution but do not eliminate any feasible integer solution. Branching is done by setting a fractional 0-1 variable x ij to 1 or 0, i.e., imposing on the one hand that entities i and j belong to the same community and on the other hand that they belong to different communities. Two linear programming subproblems are thus obtained and their solution gives bounds valid for the former and the latter case respectively. These bounds are not larger than the bound given by the solution of the subproblem on which branching took place. Branching on a subproblem stops if the value of the corresponding bound is smaller than that of the best solution known, or incumbent. As the possible number of branching choices is finite, the algorithm converges. Choosing a variable x ij with value closest to 1 2 balances improvement of both bounds. These rules and others are discussed in [START_REF] Achterberg | Branching rules revisited[END_REF]. In practice, branching is done by CPLEX. Several families of cutting planes were obtained in [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF][START_REF] Grötschel | Facets of the clique partitioning polytope[END_REF] and could be used as an alternative or complement to branching.

Note that a similar formulation is presented in [START_REF] Brandes | On modularity clustering[END_REF] but with constraints of the type x ij + x jk -2x ik ≤ 1 instead of x ij + x jkx ik ≤ 1, and so forth. While this does not change the set of feasible solutions in 0 -1 variables, the continuous relaxation in the algorithm of [START_REF] Brandes | On modularity clustering[END_REF] will be less tight than that of [START_REF] Radicchi | Defining and identifying communities in networks[END_REF]. Consequently, the algorithm is more time consuming than that of [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF] as it is shown in our computational experiments (see Section 4 below).

Unfortunately, the number of constraints of problem (5) grows rapidly with n. Not all of these constraints will be tight at the optimum. Those which are not could be deleted without changing the optimal solution but it is not known a priori which ones they are. Grötschel and Wakabayashi [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF] therefore proposed to add the constraints progressively. More precisely, they add by batches of 300 those constraints which are most violated by the current solution. In their experiments on a variety of clustering problems, they found that only a small number of such constraints were tight at the optimum and that the solution of the continuous relaxation was very often integer. Our computation experiments (see Section 4 below) show this appears also to be the case for modularity maximization.

Column generation

Column generation is a powerful technique of linear programming which allows the exact solution of linear programs with a number of columns exponential in the size of the input (there may be billions of them and, in some cases, much more). To this effect, it follows the usual steps of the simplex algorithm, apart from finding an entering column with a positive reduced cost in case of maximization which is done by solving an auxiliary problem. The precise form of this last problem depends on the type of problem under study. It is often a combinatorial optimization or a global optimization problem. It can be solved heuristically as long as a column with a reduced cost of the required sign can be found. When this is no longer the case, an exact algorithm for the auxiliary problem must be applied either to find a column with the adequate reduced cost sign, undetected by the heuristic, or to prove that there is no such column and hence the linear programming relaxation is solved. Column generation has proven to be very useful in the solution of large clustering problems, e.g., minimum sum-of-squares clustering [START_REF] Hansen | Cluster analysis and mathematical programming[END_REF][START_REF] Du Merle | An interior point algorithm for minimum sum-ofsquares clustering[END_REF][START_REF] Aloise | An Improved Column Generation Algorithm for Minimum Sum-of-squares Clustering[END_REF].

For modularity maximization clustering, as for other clustering problems with an objective function additive over the clusters, the columns correspond to all subsets of V , i.e., to all nonempty modules.

To express this problem, define a it = 1 if vertex i belongs to module t and to a it = 0 otherwise. One can then write the model as

max t∈T c t z t -C (6) s.t. t∈T a it z t = 1 ∀i = 1, . . . , n (7)
z t ∈ {0, 1} ∀t ∈ T, (8)
where c t = i j>i w ij a it a jt , i.e., the value of the module indexed by t with t = 1 . . . 2 n -1.

The objective function [START_REF] Li | Quantitative function for community detection[END_REF] of the primal problem expresses that modularity is equal to the sum of modularities of all selected modules minus a constant corresponding to the diagonal terms. The first set of constraints [START_REF] Medus | Alternative approach to community detection in networks[END_REF] expresses that each entity must belong to one and only one module and the second set of constraints that modules must be selected entirely or not at all. If the integrality constraints (8) are replaced by

z t ≥ 0,∀t ∈ T, (9)
the upper bound z t ≤ 1 being implied by constraint [START_REF] Medus | Alternative approach to community detection in networks[END_REF], one obtains a relaxation of ( 6)-( 8) which is a linear program.

Recall that to any primal linear program is associated another linear program called its dual. This dual program has as many variables as the primal has constraints and as many constraints as the primal has variables.

The dual of the relaxation of ( 6)-( 7),( 9) can be written

min n i=1 λ i -C (10) s.t. n i=1 a it λ i ≥ c t ∀t ∈ T, (11)
λ i ∈ R ∀i = 1 . . . n. (12)
The objective function [START_REF] Massen | Identifying communities within energy landscapes[END_REF] of the dual problem ( 10)-( 12) is equal to the sum of all dual variables minus a constant C. The constraints [START_REF] Fortunato | Resolution limit in community detection[END_REF] express that the sum of dual variables associated with the entities of any community must be at least as large as its modularity. Finally, the constraints (12) express the fact that the dual variables are unrestricted in sign.

From the duality theorem of linear programming, the optimal solutions (z * 1 , z * 2 , . . . , z * T ) of the primal and (λ * 1 , λ * 2 , . . . , λ * n ,) of the dual have the same value:

t∈T c t z * t = n i=1 λ * i . (13)
Moreover, any feasible solution of the primal has smaller or equal value than any solution of the dual. Geometrically, the sets of feasible solutions of the primal and dual are polyhedra. Finding the optimal solution of the primal or of the dual amounts to finding an optimal vertex and the corresponding cone in either of these polyhedra. Problem ( 6)-( 7),( 9) is called the master problem. It is the relaxation of a partitioning problem which can in principle be solved by a package such as CPLEX. However, as the number of columns is exponential, this is possible only for very small n. To overcome these difficulties, one resorts to column generation. Following that approach a reduced master problem with considerably fewer columns is solved instead of the full master problem ( 6)-( 7), [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. As usual in the branch and bound approach to mixed integer programming, one first solves the continuous relaxation of this restricted master problem. One begins with a relaxed problem containing some feasible columns and possibly artificial variables. Then improving columns will be added progressively. Finding such column(s) is the auxiliary problem whose role is to find a column with positive (negative) reduced cost in case of maximization (minimization). For problem ( 6)-( 7), [START_REF] Newman | Finding and evaluating community structure in networks[END_REF], the reduced cost associated with column t will be equal to c ti λ i a it where the λ i are the current values of the dual variables of the continuous relaxation of problem ( 6)-( 7), [START_REF] Newman | Finding and evaluating community structure in networks[END_REF]. Replacing the coefficients a it by binary variables y i leads to the following expression of the auxiliary problem: max

y∈B n i j>i w ij y i y j - i λ i y i
This is a quadratic program in 0-1 variables with a 100% dense matrix of coefficients. Many algorithms and numerous heuristics have been proposed to solve it. In our experiments, we use a Variable Neighborhood Search (vns) heuristic [START_REF] Mladenovic | Variable neighborhood search[END_REF][START_REF] Hansen | Variable neighborhood search: Principles and applications[END_REF] as long as it can find a column with positive reduced cost. vns is a metaheuristic, i.e., a framework for building heuristics, based on the idea of systematic change of neighborhood during the search. It explores progressively larger neighborhoods of the incumbent (or best known) solution in a probabilistic way. Therefore, often favorable characteristics of the incumbent will be kept and used to obtain promising neighboring solutions. vns applies a local search routine repeatedly to get from these neighboring solutions to local optima.

When vns fails to find an improving column, we use as exact method a simple branch and bound algorithm [START_REF] Hansen | A simple enumerative algorithm for unconstrained 0-1 quadratic programming[END_REF] or a recent algorithm using bounds based on semidefinite programming [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF].

It is well known that column generation algorithms suffer from slow convergence particularly when the optimal solution is degenerate, i.e., when such a solution has many variables equal to 0, which is the case for clustering problems. Column generation algorithms also suffer from the plateau effect, i.e., the optimal solution keeps the same value for several or many iterations [START_REF] Vanderbeck | A generic view of dantzig-wolfe decomposition in mixed integer programming[END_REF].

To alleviate these defects, one can use a variant of the stabilization methods for column generation due to du Merle et al. [START_REF] Du Merle | Stabilized column generation[END_REF], which we call focussed column generation.

The principle is to identify, from a heuristic solution, a small region in the space of dual variables hopefully containing the optimal solution (λ * 1 , λ * 2 , . . . , λ * n ,). Then departures from this zone are penalized.

To that effect, we first seek a good heuristic solution of the modularity maximization problem. Recently, Noack and Rotta [START_REF] Noack | Multi-level algorithms for modularity clustering[END_REF] compared experimentally codes for eight heuristics. We used the SS+ML heuristic which is the best according to their experiments, i.e., Single-Step Greedy coarsening by Significance with Multi-Level Fast Greedy refinement (SS+ML). This heuristic solution can be further improved by a local application of modularity maximization to each pair of modules at a time. Then, a dual solution is derived from this last heuristic primal solution, and intervals hopefully containing the optimal values are determined. This is done by computing the increase (resp. the decrease) of the objective function value if a vertex is duplicated (resp. removed). Penalties for getting out of these intervals are imposed and progressively diminished until they get down to 0. Details on this method are given in [START_REF] Du Merle | Stabilized column generation[END_REF].

3 Modularity maximization by mixed 0-1 quadratic programming

Direct formulation

Maximizing modularity by the clique partitioning approach discussed in Section 2 has a drawback: it replaces a usually sparse matrix of coefficients by a 100% dense one. An alternative approach is to work directly with a graph G = (V, E) instead of the complete graph K n . This was done by Xu, Tsoka and Papageorgiou [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] and leads to a 0-1 mixed integer quadratic problem whose continuous relaxation is convex, and which can therefore be solved by CPLEX. We next recall the main elements of Xu et al.'s model as they provide the necessary background for a new column generation algorithm described in the second part of this section. Considering again the definition of Q as a sum over modules of their modularities rewrite Q as

Q = s [a s -e s ] = s m s m - D s 2m 2 (14)
where m s denotes the number of edges in module s, i.e., the subgraph induced by V s and D s denotes the sum of degrees k i of the vertices of module s. Binary variables are then used to identify the modules to which each vertex and each edge belongs. To this effect, list all edges with a single index r = 1, 2, . . . , m. Then, introduce the following variables, X rs = 1 if edge r belongs to module s 0 otherwise for r = 1, 2, . . . m and s = 1, 2, . . . M and

Y is = 1 if vertex i belongs to module s 0 otherwise.
The number of edges and sum of vertex degrees can then be expressed as:

m s = r X rs and D s = i k i Y is .
A second series of constraints express that each vertex belongs to exactly one module:

s Y is = 1 ∀i = 1, 2, . . . n.
A third series of constraints express that any edge r = {v i , v j } with end vertices indiced by i and j can only belong to module s if both of those end vertices belong to that module.

G-2010-36

Les Cahiers du GERAD

X rs ≤ Y is ∀r = {v i , v j } ∈ E X rs ≤ Y js ∀r = {v i , v j } ∈ E
Note that these constraints are part of Fortet's [START_REF] Fortet | L'algèbre de Boole et ses applications en Recherche Opérationnelle[END_REF] linearization of quadratic 0-1 programs. They impose that X rs = 0 if either Y is or Y js or both are equal to 0. They do not impose that X is is equal to 1 is Y is = Y js = 1. Although the constraints X rs ≥ Y is + Y js -1 could be added to impose that, this is not necessary as adding an edge between vertices i and j when they are in the same module increases m s and hence Q.

The number of modules is a priori unknown; indicator variables u s = 1 if module s is nonempty and u s = 0 otherwise are used. Then constraints

u s ≤ u s-1 ∀s ∈ 2, 3, . . . S,
where S is an upper bound on the number of modules, are added and express that module s can be nonempty only if module s -1 is so.

Consequently, r X rs ≥ u s and r X rs ≤ (n -s + 1)u s .
The value ns + 1 in the latter constraints is due to the fact that each of the modules 1, 2, . . . s -1 must be nonempty. In fact, Xu et al. use more general parametric formulae which allow imposing lower and upper bounds on the cardinality of the modules.

It is well-known that for any given solution to a clustering problem, alternative equivalent solutions can be obtained by simply re-indexing clusters. For an optimal solution with M modules, M ! equivalent solutions exist. Symmetry-breaking constraints for clustering problems can be found in the literature [START_REF] Klein | Optimal clustering: A model and method[END_REF][START_REF] Sherali | A global optimization rlt-based approach for solving the hard clustering problem[END_REF][START_REF] Plastria | Formulating logical implications in combinatorial optimisation[END_REF], and were used in [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF]. Such symmetry breaking constraints notwithstanding, this problem has M (n + m) binary variables and M continuous variables, subject to M (1 + 2m + 2n) linear constraints.

Column generation reformulation

Once again, a linear programming master problem will be solved by the simplex algorithm where the entering column will be determined by solving an auxiliary problem. The master problem will be the same as in the previous column generation algorithm, i.e., its equations are given in ( 6)-( 7),(9) of Section 2.2. The auxiliary problem will be different from the quadratic 0-1 program used in the clique partitioning formulation. It will be close to the formulation of Xu et al. [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] described in Section 3.1, but much simpler. As a single community is to be determined at a time, it can be written as follows:

max x∈B n ,D∈R r x r m - D 2m 2 - i λ i y i s.t. D = i k i y i x r ≤ y i ∀r = {i, j} ∈ E x r ≤ y j ∀r = {i, j} ∈ E.
As before, edges are indexed by r and vertices by i (or j). Variable x r is equal to 1 if edge r belongs to the community which maximizes the objective function and to 0 otherwise. Similarly, y i is equal to 1 if the i th vertex belongs to the community and 0 otherwise. The objective function is equal to the modularity of the community to be determined minus the scalar product of the current value λ i of the dual variables times the indicator variable y i . Observe that this last term is the same as in the objective function of the auxiliary problem for clique partitioning. This is a mixed integer quadratic problem with n + m binary variables and 1 continuous variable, in the objective function, subject to 2m + 1 linear constraints. In the objective function there is a single concave nonlinear term. Clearly, the size of this auxiliary problem is much smaller than that of the direct formulation described in Section 3.1, particularly for large number of communities M . This auxiliary problem is first solved with a vns heuristic as long as a column with a positive reduced cost can be found. When this is no more the case, CPLEX is called to find such a column or prove that there are no more.

Computational comparison

To compare the four algorithms described in the previous sections, we selected 11 test problems from the modularity maximization literature. Their names, orders and sizes are given in Table 1. The data of these test problems can be found in various databases mentioned in a recent paper from Noack and Rotta. Note that we always assume unit edge weights.

Results are given in Table 2. The first column gives the problem ID as in Table 1, the next six columns give the modularity value (Q) and the number of modules (M ) for: (a) the solution obtained by Noack-Rotta's heuristic (NR Sol); (b) the improved solution obtained from Noack-Rotta's solution (Imp Sol); (c) the optimal solution (Opt Sol).

The remaining columns summarize the performance of the different exact algorithms:

• CPRG: clique partitioning row generation of [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF];

• CPCG-HJM: clique partitioning column generation with the exact method of [START_REF] Hansen | A simple enumerative algorithm for unconstrained 0-1 quadratic programming[END_REF] for the auxiliary problem; • CPCG-BE: clique partitioning column generation with the exact method of [START_REF] Billionnet | Using a mixed integer quadratic programming solver for the unconstrained quadratic 0-1 problem[END_REF] for the auxiliary problem; 0-1 mixed integer programming formulation of [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF]; • 0-1 MICG: 0-1 mixed integer column generation.

For CPRG column, the CPU time is reported when possible or an "OM" is used to indicate that the program reaches the memory limit. The formulation of [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF] is implemented in AMPL [START_REF] Fourer | The AMPL Book[END_REF] and solved using the "lazy constraints" feature of CPLEX. We do not present the results for the formulation of [START_REF] Brandes | On modularity clustering[END_REF] as it is much more time consuming. For example, the solution of Problem 1 takes 674 seconds instead of 0.23 second.

For the 0-1 MIQP column, we give the CPU time reported in [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] or indicate by "NA" the problem not considered in [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF]. It was not possible for us to obtain comparable computing time as the parameters setting used to obtain these results were not detailed in the paper.

For the column generation algorithms, the table contains CPU time as well as, in parentheses, the number of calls for the exact algorithm for the auxiliary problem followed by the number of calls to the heuristic. The term "OT" means that the optimal solution could not be found after more than 100 000 seconds.

G-2010-36

Les Cahiers du GERAD All results are in seconds of CPU. Except for 0-1 MIQP, all results were obtained on a dual processor computer Intel Pentium computer with 3.20 GHz, 2 Mb cache memory, 3 GB RAM running under Linux. ILOG CPLEX 10.110 was used for the linear programming part of all algorithms.

From Table 2, it appears that:

• surprisingly, in view of the many heuristics proposed for modularity maximization, those applied here seldom reach the optimal value: in 2 cases out of 11 of Noack and Rotta and 5 cases out of 11 for the improved heuristic; • both row-generation and column-generation algorithm based on reformulation of modularity maximization as a clique partitioning problem are competitive for small instances, but become too time or memory consuming for larger ones due to the rapid increase in the number of variables and constraints; • for the largest instances, the 0-1 MICG algorithm appears to be the only one able to find an optimal solution; • the reduction in resolution times in comparison of the results of [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] are very substantial: the use of 0-1 MICG algorithm divides these times by a factor of 3 to 159; • the size of the largest problem solved is raised from 105 entities to 512 entities.

In summary, the column generation algorithm reformulated from [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF] direct formulation (0-1 MICG) appears to be the best choice since its computing time is comparable for small instances, lower for medium instances and is the only algorithm able to solve large instances to optimality.

Conclusions

In this paper, we have studied exact maximization of the modularity of a network according to the definition of Newman and Girvan. Two approaches were proposed previous to this work: on the one hand the row generation algorithm of [START_REF] Grötschel | A cutting plane algorithm for a clustering problem[END_REF] to which modularity maximization can be reduced and on the other hand the direct formulation of Xu et al. [START_REF] Xu | Finding community structures in complex networks using mixed integer optimization[END_REF]. We have proposed column generation algorithms based on these two approaches and performed a computation comparison on a series of well-known problems from the literature. These results show that: (i) while row generation is fast for small problems, it is outperformed by column generation for larger ones; (ii) the column generation algorithm based on direct formulation performs best for the larger problems; (iii) size of the larger problems has substantially augmented: the largest problem solved to date has 512 entities versus 105 entities before; (iv) computation times are substantially reduced. While it would of course be desirable to solve exactly larger problems, the results obtained by the proposed algorithm already appear to be useful for the several reasons detailed in the introduction.

It is easy to see that the first two approaches can also be applied to several variants of the standard modularity definition, i.e., weighted networks, directed networks, and networks which are both directed and weighted [START_REF] Fortunato | Community detection in graphs[END_REF]. Moreover, the column generation framework is easily adaptable to yield exact algorithms for other definitions of communities than those of Newman and Girvan, allowing comparison which do not depend on the specific heuristic used.

Among topics for future research are (i) the design of better heuristics, (ii) the use of cutting planes in the solution of the master problem and/or the auxiliary problem, (iii) a further study of the properties of the auxiliary problem which might lead to improvements in its resolution. The importance of this last point stems from the fact that the proportion of the resolution time devoted to the resolution of the auxiliary problem tends to increase with problem size.

Table 1 :

 1 Order and size of 11 test problems

	Problem ID	Name	n	m
	1	Zachary's karate club	34	78
	2	Dolphins social network	62	159
	3	Les Misérables	77	254
	4	A00 main	83	135
	5	protein p53	104	226
	6	Books about US politics	105	441
	7	American College Football	115	613
	8	A01 main	249	635
	9	USAir97	332	2126
	10	netscience main	379	914
	11	Electronic Circuit (s838)	512	819

Table 2 :

 2 Results of the comparison between algorithms for modularity maximization. Best results are reported in boldface.

	Pb.	NR Sol	Imp Sol	Opt Sol	CPRG	CPCG-	0-1 MIQP	0-1 MICG
	ID	Q	M	Q	M	Q	M		HJM	BE
	1	0.4198 4 0.4198 4 0.4198 4	0.23	0.18 (3/7)	0.73 (3/7)	1.03	0.34	(3/8)
	2	0.5238 4 0.5285 5 0.5285 5	9.70	7.49 (1/13)	18.69 (1/13)	197.89	7.75	(1/13)
	3	0.5600 6 0.5600 6 0.5600 6	5.07	6.89 (1/13)	12.92 (1/13)	55.58	7.26	(1/13)
	4	0.56749 9 0.5283 9 0.5309 9	17.09	3.57 (1/26)	262.13 (1/26)	NA	3.66	(1/26)
	5	0.5322 7 0.5350 6 0.5351 7 4164.15 11.33 (1/44)	351.87 (1/44)	1844.31	11.60	(1/44)
	6	0.5269 4 0.5272 5 0.5272 5	663.05 318.99 (2/13) 32287.23 (2/13)	NA	45.65	(2/13)
	7	0.6002 10 0.6046 10 0.6046 10 282.58	OT	593.01 (1/8)	NA	249.41	(1/8)
	8	0.6203 12 0.6203 12 0.6329 14	OM	OT	OT	NA	1014.48 (12/145)
	9	0.3658 6 0.3660 6 0.3682 6	OM	OT	OT	NA 16216.77	(3/104)
	10	0.8474 19 0.8485 19 0.8486 19	OM	OT	OT	NA	1615.14	(12/66)
	11	0.8162 16 0.8166 16 0.8194 12	OM	OT	OT	NA	7655.56 (140/225)

Financial support by Grants ANR 07-JCJC-0151 "ARS", Digiteo 2009-14D "RM-NCCO", Digiteo 2009-55D "ARM", are greatefully acknowledged. P.H. acknowledges support from Digiteo Foundation. S.P. has been supported by the nserc (Natural Sciences and Engineering Research Council of Canada) grant 327435-06. G.C., P.H. and S.P. would like to thank the research office of HEC Montréal for the participation to the 2009-2010 research workshop of HEC Montréal.