Etienne Côme

Pierre Latouche

Model selection and clustering in stochastic block models with the exact integrated complete data likelihood

Keywords: Random graphs, stochastic block models, integrated classification likelihood

The stochastic block model (SBM) is a mixture model used for the clustering of nodes in networks. It has now been employed for more than a decade to analyze very different types of networks in many scientific fields such as Biology and social sciences. Because of conditional dependency, there is no analytical expression for the posterior distribution over the latent variables, given the data and model parameters.

Therefore, approximation strategies, based on variational techniques or sampling, have been proposed for clustering. Moreover, two SBM model selection criteria exist for the estimation of the number K of clusters in networks but, again, both of them rely on some approximations. In this paper, we show how an analytical expression can be derived for the integrated complete data log likelihood. We then propose an inference algorithm to maximize this exact quantity. This strategy enables the clustering of nodes as well as the estimation of the number clusters to be performed at the same time and no model selection criterion has to be computed for various values of K. The algorithm we propose has a better computational cost than existing inference techniques for SBM and can be employed to analyze large networks with ten thousand nodes. Using toy and true data sets, we compare our work with other approaches.

Introduction

There is a long history of research on networks which goes back to the earlier work of Moreno [START_REF] Moreno | Who shall survive?: a new approach to the problem of Human interrelations[END_REF]. Because they are simple data structures yet capable of representing complex systems, they are used in many scientific fields [START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF][START_REF] Palla | Quantifying social group evolution[END_REF]. Originally considered in social sciences [START_REF] Fienberg | Categorical data analysis of single sociometric relations[END_REF] to characterize relationships between actors [START_REF] Holland | Stochastic block models : First steps[END_REF][START_REF] Boulet | Batch kernel som and related laplacian methods for social network analysis[END_REF], networks are also used to describe neural networks [START_REF] White | The structure of the nervous system of the nematode caenorhabditis elegans, Philosophical Transactions[END_REF], powergrids [START_REF] Watts | Collective dynamics of small-world networks[END_REF], and the Internet [START_REF] Adamic | The political blogosphere and the 2004 us election[END_REF][START_REF] Zanghi | Fast online graph clustering via erdos-renyi mixture[END_REF]. Other examples of real networks can be found in Biology with the use of regulatory networks to describe the regulation of genes by transcriptional factors [START_REF] Milo | Network motifs: simple building blocks of complex networks[END_REF] or metabolic networks to represent pathways of biochemical reactions [START_REF] Lacroix | Motif search in graphs: pplication to metabolic networks[END_REF]. As the number of networks used in practice has been increasing, a lot of attention has been paid on developing graph clustering algorithms to extract knowledge from their topology. Existing methods usually aim at uncovering very specific patterns in the data, namely communities or disassortative mixing. For an exhaustive review, we refer to [START_REF] Goldenberg | A survey of statistical network models[END_REF].

Most graph clustering algorithms look for communities, where two nodes of the same community are more likely to be connected than nodes of different communities. These techniques [START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF][START_REF] Newman | Modularity and community structure in networks[END_REF] often maximize the modularity score proposed by Girvan and Newman [START_REF] Girvan | Community structure in social and biological networks[END_REF] for clustering. However, recent work of Bickel and Chen [START_REF] Bickel | A non parametric view of network models and newmangirvan and other modularities[END_REF] showed that they were asymptotically biased and tended to lead to the discovery of an incorrect community structure, even for large graphs. Alternative strategies, see for instance [START_REF] Krivitsky | Representing degree distributions, clustering, and homophily in social networks with latent cluster random effects models[END_REF], are generally related to the probabilistic model of Handcock, Raftery and Tantrum [START_REF] Handcock | Model-based clustering for social networks[END_REF] which generalizes the work of Hoff, Raftery and Handcock [START_REF] Hoff | Latent space approaches to social network analysis[END_REF]. Nodes are first mapped into a a latent space and then clustered depending on their latent positions. Community structure algorithms are commonly used for affiliation network analysis. As mentioned in [START_REF] Newman | Mixture models and exploratory analysis in networks[END_REF], other graph clustering algorithms aim at uncovering dissasortative mixing in networks where, contrary to community structure, nodes mostly connect to nodes of different clusters. They are particularly suitable for the analysis of bipartite or quasi bipartite networks [START_REF] Estrada | Spectral measures of bipartivity in complex networks[END_REF].

In this paper, we consider the stochastic block model (SBM) proposed by Nowicki and Snijders [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] which is a probabilistic generalization [START_REF] Fienberg | Categorical data analysis of single sociometric relations[END_REF][START_REF] Holland | Stochastic block models : First steps[END_REF] of the work of White, Boorman and Breiger [START_REF] White | Social structure from multiple networks. i. blockmodels of roles and positions[END_REF]. As pointed out by Daudin, Picard and Robin [START_REF] Daudin | A mixture model for random graph[END_REF], SBM can be seen as a mixture model for graphs. It assumes that nodes are spread into K clusters and uses a K × K matrix Π to describe the connection probabilities between pairs of nodes. No assumption is made on Π such that very different structures can be taken into account. In particular, as shown in [START_REF] Latouche | Advances in Data Analysis Data Handling and Business Intelligence[END_REF], contrary to the methods mentioned previously, SBM can be used to retrieve both communities and disassortative mixing in networks.

Many extensions have been developed to overcome some limits of the standard SBM. For example, Mariadassou, Robin and Vacher [START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF] introduced recently a probabilistic framework to deal with valued edges, allowing covariates to be taken into account. While the first model they proposed explains the value of an edge, between a pair of nodes, through their clusters only, the second and third approaches do account for covariates through Poisson regression models. This framework is relevant in practice because extra information on the edges is sometimes available, such as phylogenetic distances in host-parasite networks or amounts of energy transported between nodes in powergrids.

Another drawback of SBM is that it assumes that each node belongs to a single cluster while many objects in real world applications belong to several groups or communities [START_REF] Latouche | Overlapping stochastic block models with application to the french political blogosphere[END_REF]. To tackle this issue Airoldi et al. [START_REF] Airoldi | Mixed-membership stochastic blockmodels[END_REF] proposed the mixed membership stochastic block model (MMSBM) [START_REF] Airoldi | Mixed membership stochastic block models for relational data with application to protein-protein interactions[END_REF][START_REF] Airoldi | Mixed membership analysis of high-throughput interaction studies: relational data[END_REF]. A latent variable π i , drawn from a Dirichlet distribution, is associated to each node i of a network. Given a pair (i, j) of nodes, two binary latent vectors Z i→j and Z i←j are then considered.

The vector Z i→j is assumed to be sampled from a multinomial distribution with parameters (1, π i ) and describes the cluster membership of i in its relation towards j.

By symmetry, Z i←j is drawn from multinomial distribution with parameters (1, π j ) and characterizes the cluster membership of j in its relation towards i. Thus, in MMSBM, since each node can have different latent vectors through its relations towards other nodes, it can belong to several clusters. The connection probability between i and j is finally given by p ij = Z i→j BZ i←j . The overlapping stochastic block model (OSBM) was proposed by Latouche, Birmelé and Ambroise [START_REF] Latouche | Overlapping stochastic block models with application to the french political blogosphere[END_REF] as an alternative probabilistic model for networks allowing overlapping clusters. Contrary to MMSBM, edges are influenced by the fact that some nodes belong to multiple clusters. Thus, each node i is characterized by a binary latent vector Z i sampled from a product of Bernoulli distributions. An edge between nodes i and j is then drawn from a Bernoulli distribution with parameter p ij=g(a Z i ,Z j ) . The function g(•)

is the logistic function while a Z i ,Z j is a real variable describing the interactions between the nodes, depending on the different clusters they are associated with. It is given by a Z i ,Z j = Z i WZ j + Z i U + V Z j + W * . Finally, we mention the work of Karrer and Newman [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF] who proposed an interesting extension of SBM to deal with node degree heterogeneity inside clusters. The model deals with valued edges and includes another set of parameters describing vertices attractiveness. Using the right constraints the model is identifiable (up to permutations of clusters) and the atractivity parameters can be directly related to vertices degree. This work was extended to oriented networks in [START_REF] Zhu | Oriented and degree-generated block models: Generating and inferring communities with inhomogeneous degree distributions[END_REF] and finally tools for model selection between different models are derived in [START_REF] Yan | Model selection for degree-corrected block models[END_REF].

In this paper, we wont consider the extensions for SBM we mentionned, we will rather focus on the standard SBM. Our goal here is not to propose new extensions, allowing a SBM like model to be applicable on specific types of networks, or to introduce new latent structures. Conversely, considering the standard SBM, which has been widely used in practice for network analysis, for more than a decade, we aim at developping a new optimization procedure, improving over existing inference strategies. In SBM, the posterior distribution over the latent variables, given the parameters and the observed data, cannot be factorized due to conditional dependency. Therefore, optimization techniques such as the expectation maximization (EM) algorithm cannot be used directly for clustering. To tackle this issue, Daudin, Picard and Robin [START_REF] Daudin | A mixture model for random graph[END_REF] proposed an approximation method based on a variational EM algorithm. Note that an online version of this algorithm exists [START_REF] Zanghi | Strategies for online inference of network mixture[END_REF]. A Bayesian framework was also considered by Nowicki and Snijders [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] where conjugate priors for the model parameters were introduced. Again, because the posterior distribution over the model parameters, given the data, is not tractable, approximation techniques were employed for inference. Thus, Nowicki and Snijders [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] used a Gibbs sampling procedure while Latouche, Birmelé and Ambroise [START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF] relied on a variational Bayes EM algorithm. To our knowledge, only two model selection criteria, the integrated classification likelihood (ICL) and the integrated likelihood variational Bayes (ILvb) have been developed for SBM in order to estimate the number K of clusters in networks. Standard criteria such as the Akaike information criterion (AIC) or bayesian information criterion (BIC) cannot be used because they rely on the SBM observed data log likelihood which is not tractable in practice (see for instance [START_REF] Latouche | Advances in Data Analysis Data Handling and Business Intelligence[END_REF]). ICL was originally developed by Biernacki, Celeux and Govaert [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] for Gaussian mixture models and then adapted by Daudin, Picard and Robin [START_REF] Daudin | A mixture model for random graph[END_REF] to SBM. It is based on Laplace and Stirling approximations of the integrated complete data log likelihood. As shown in [START_REF] Biernacki | Exact and monte carlo calculations of integrated likelihoods for the latent class model[END_REF], it tends to miss some important structures present in the data for small data sample, because of the asymptotic approximations. To tackle this drawback, Latouche, Birmelé and Ambroise proposed in [START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF] the ILvb criterion which relies on a variational Bayes approximation of the integrated observed data log likelihood.

In this paper, we show how an analytical expression of the integrated complete data log likelihood can be obtained in a Bayesian framework and that no asymptotic approximation is required. We call the corresponding criterion ICL ex where ex stands for exact. We then propose a greedy inference algorithm which maximizes this exact quantity. The strategy has three advantages compared to existing approaches. First, it maximizes an analytical criterion directly derived from SBM, while variational techniques for instance rely on lower bounds for approximation. Thus, the lower bound of the variational EM algorithm proposed by Daudin, Picard and Robin [START_REF] Daudin | A mixture model for random graph[END_REF] approximates the observed data log likelihood, while Latouche, Birmelé and Ambroise [START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF] introduced a lower bound to estimate the integrated observed data log likelihood. Second, ICL ex only depends on all the latent binary vectors Z i , stored in the matrix Z, and the number K of clusters, not on the model parameters which are marginalized out. Therefore, the optimization task focus on (K, Z) and is purely combinatorial. When using the Gibbs algorithm [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF], the successive samples for Z and model parameters are highly correlated. As a consequence, nodes tend to be stuck in clusters, after a few iterations. Similar remarks could be made for the variational EM and variational Bayes EM algorithms. In our case, because the parameters are marginalized out (collapsed) the method is expected to explore more easily the latent space of Z. This property is at the core of collapsing methods (for more details, we refer to [START_REF] Liu | The collapsed gibbs sampler in bayesian computations with applications to the gene regulation problem[END_REF]). Finally, our strategy enables the clustering of nodes as well as the estimation of the number of clusters to be performed at the same time and no model selection criterion has to be computed for various values of K. Starting from a complex model with K = K up clusters, (K up being an upper bound for K), the proposed algorithm swaps labels until ICL ex reaches a local maximum. During the process, clusters may disappear, i.e. their cardinality reaches zero. Such an approach leads to a simple and time attractive algorithm with complexity of O(L + N K 2 up ), with L the total number of edges in the network and N the number of vertices.

As we shall see through a series of experiments, the greedy algorithm takes benefit of computing the exact ICL and improves over existing methods, both in terms of clustering and model selection. It can also deal with large networks with tens of thousands of vertices.

The stochastic block model

We consider a binary network with N nodes represented by its adjacency matrix X such that X ij = 1 if there is an edge from node i to node j, 0 otherwise. In this paper, we focus on directed networks, i.e. relations are oriented. Therefore X is not symmetric. Moreover, we do not consider any self loop, that is an edge from a node to itself. We emphasize that all the optimization equations derived in this work can easily be adapted to deal with undirected networks or to take into account self loops.

Model and notations

The stochastic block model (SBM) introduced by Nowicki and Snijders [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] assumes that the nodes are spread into K clusters with connectivity patterns described by a K × K matrix Π. The cluster of each node is given by its binary membership vector Z i sampled from a multinomial distribution :

Z i ∼ M(1, α = (α 1 , . . . , α K )), K k=1 α k = 1,
such that Z ik = 1 if i belongs to cluster k and zero otherwise. Contrary to the work of Latouche, Birmelé and Ambroise [START_REF] Latouche | Overlapping stochastic block models with application to the french political blogosphere[END_REF], each node belongs to a single cluster, that is K k=1 Z ik = 1, ∀i. Given the vectors Z i and Z j , an edge between node i and j is then drawn from a Bernoulli distribution with probability Π kl :

X ij |Z ik Z jl = 1 ∼ B(Π kl ).
This leads to a simple yet flexible generative model for networks. First, all the vectors Z i are sampled independently. We denote Z the binary N × K matrix storring the Z i s as raw vectors :

p(Z|α) = N i=1 M(Z i ; 1, α) = N i=1 K k=1 α Z ik k . (1)
Then, given the latent structure Z, all the edges in X are drawn independently :

p(X|Z, Π) = N i =j p(X ij |Z i , Z j , Π) = N i =j K k,l B(X ij ; Π kl ) Z ik Z jl = N i =j K k,l Π X ij kl (1 -Π kl ) 1-X ij Z ik Z jl . (2)

Integrated classification likelihood criteria

In this paper, we consider the integrated complete data log likelihood log p(X, Z|K)

in order to focus on the inference of Z and K from the observed data X, all the SBM parameters (α, Π) being integrated out. We first recall existing approximations and then show in Section 2.2.2 how an analytical expression of this quantity can be derived.

Asymptotic ICL criterion

When considering a factorized prior distribution p(α, Π|K) = p(α|K)p(Π|K)

over the model parameters, as in [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF], the integrated complete data log likelihood easily decomposes into two terms:

log p(X, Z|K) = log α,Π p(X, Z, Π, α|K)dαdΠ = log Π p(X|Z, Π, K)p(Π|K)dΠ α p(Z|α, K)p(α|K)dα = log p(X|Z, K) + log p(Z|K). (3)
However, for an arbitrary choice of the priors p(α|K) and p(Π|K), the marginal distributions p(X|Z, K) as well as p(Z|K) are usually not tractable and (3) does not have any analytical form. To tackle this issue, Daudin, Picard and Robin [START_REF] Daudin | A mixture model for random graph[END_REF] proposed an asymptotic approximation of log p(X, Z|K), so called integrated classification likelihood criterion (ICL). Note that ICL was originally proposed by Biernacki, Celeux and Govaert [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF] for Gaussian mixture models. It was then adapted by Biernacki, Celeux and Govaert [START_REF] Biernacki | Exact and monte carlo calculations of integrated likelihoods for the latent class model[END_REF] to mixtures of multivariate multinomial distributions and to the SBM model by Daudin, Picard and Robin [START_REF] Daudin | A mixture model for random graph[END_REF]. In the case we consider of a directed graph without self-loop, ICL is given by:

ICL(Z, K) ≈ log p(X, Z|K) = max α,Π log p(X, Z|α, Π, K) - 1 2 K 2 log (N (N -1)) - K -1 2 log(N ).
For an extensive description of the use of Laplace and Stirling approximations to derive the ICL criterion, we refer to [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. Since it approximates the integrated complete data log likelihood, ICL is known to be particularly suitable when the focus is on the clustering task and not on the estimation of the data density. However, as shown in [START_REF] Biernacki | Exact and monte carlo calculations of integrated likelihoods for the latent class model[END_REF][START_REF] Mariadassou | Uncovering latent structure in valued graphs: a variational approach[END_REF], it tends to miss some important structures present in the data because of the (asymptotic) approximations.

We emphasize that ICL is only used in the literature as a model selection criterion. In practice, a clustering method such as an EM like algorithm for instance is employed to obtained some estimates Z of Z, for various values of the number K of classes. ICL is then computed for every pair ( Z, K) and the pair ( Z * , K * ) is chosen such that the criterion is maximized. Thus, ICL is optimized only through the results ( Z, K) produced by the clustering algorithm. Conversely, after having given an analytical expression ICL ex of the integrated complete data log likelihood in the next section, we will show in Section 3 how to optimize directly ICL ex with respect to Z and K.

Exact ICL criterion

We rely on the same Bayesian framework as in [START_REF] Nowicki | Estimation and prediction for stochastic blockstructures[END_REF] and [START_REF] Latouche | Advances in Data Analysis Data Handling and Business Intelligence[END_REF]. Thus, we consider non informative conjugate priors for the model parameters α and Π. Since α, describing the cluster proportions, parametrizes a multinomial distribution (1), we rely on a Dirichlet prior distribution:

p(α) = Dir α; n 0 = (n 0 1 , . . . , n 0 K ) .
A common choice consists in fixing the hyperparameters to 1/2, i.e. n 0 k = 1/2, ∀k. Such a distribution corresponds to a non informative Jeffreys prior which is known to be proper [START_REF] Jeffreys | An invariant form for the prior probability in estimations problems[END_REF]. A uniform distribution can also be obtained by setting the hyperparameters to 1.

Moreover, since the presence or absence of an edge between nodes is sampled from a Bernoulli distribution, we consider independent Beta prior distributions to model the connectivity matrix Π:

p(Π) = K k,l Beta(Π kl ; η 0 kl , ζ 0 kl ).
Again, if no prior information is available, all hyperparameters η 0 kl and ζ 0 kl can be set to 1/2 or 1 to obtain a Jeffreys or uniform distribution.

With these choices of conjugate prior distributions over the model parameters, the marginal distributions p(X|Z, K) as well as p(Z|K) in (3) have analytical forms, and so has the integrated complete data log likelihood, as proved in AppendixA. We call ICL ex the corresponding criterion, where ex stands for exact. It is given by:

ICL ex (Z, K) = log p(X, Z|K) = K k,l log Γ(η 0 kl + ζ 0 kl )Γ(η kl )Γ(ζ kl ) Γ(η kl + ζ kl )Γ(η 0 kl )Γ(ζ 0 kl ) + log Γ( K k=1 n 0 k ) K k=1 Γ(n k ) Γ( K k=1 n k ) K k=1 Γ(n 0 k )
, where the components n k are:

n k = n 0 k + N i=1 Z ik , ∀k ∈ {1, . . . , K},
and can be seen as pseudo counters of the number of nodes in each class. Moreover, the parameters (η kl , ζ kl ) are given by:

η kl = η 0 kl + N i =j Z ik Z jl X ij , ∀(k, l) ∈ {1, . . . , K} 2 ,
and

ζ kl = ζ 0 kl + N i =j Z ik Z jl (1 -X ij ), ∀(k, l) ∈ {1, . . . , K} 2 .
They represent pseudo counters of the number of edges and non-edges connecting nodes of class k to nodes of class l, respectively.

Note that the ICL ex criterion is related to the variational Bayes approximation of the integrated observed data log likelihood log p(X|K) proposed by Latouche, Birmelé and Ambroise [START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF]. The key difference is that the parameters (n k , η kl , ζ kl ) in ICL ex depend on the hard assignment Z of nodes to classes and not on approximated posterior probabilities τ . Moreover, ICL ex does not involve any entropy term as in [START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF].

Greedy optimization

Since the model parameters have been marginalized out, the ICL ex criterion only involves the cluster indicator matrix Z whose dimensionality depends on the number K of clusters. Thus, this integrated likelihood is only a function of a partition P,

i.e. an assignment of the vertices to clusters. Looking directly for a global maximum of ICL ex is not feasible because it involves testing every possible partition of the vertices with various values of K. However, this is a combinatorial problem for which heuristics exist to obtain local maxima. In this paper, we rely on greedy heuristics which have been shown to scale well with sample sizes [START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF]. These approaches have already been used for graph clustering using ad-hoc criteria such as modularity [START_REF] Newman | Fast algorithm for detecting community structure in networks[END_REF][START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] and are reminiscent of the well known iterated conditional modes algorithm of Besag [START_REF] Besag | On the statistical analysis of dirty pictures (with discussions)[END_REF] used for maximum a posteriori estimation in Markov random fields.

The algorithm (see Algorithm 1) starts with a SBM model with K = K up clusters, K up being an upper bound for the number of clusters. K up is assumed to be given as an input along with a N × K up matrix Z. In practice, K up is set to a large value using user knowledge on the problem at hand, while Z can be initialized with the methods described in the next section. The algorithm then cycles randomly through all the vertices of the network. At each step, a single node i is considered while all the membership vectors Z j for j = i are hold fixed. If i is currently in cluster g, the method looks for every possible label swapping, i.e. removing i from cluster g and assigning it to a cluster h = g, and computes the corresponding change ∆ g→h in the ICL ex criterion. Note that ∆ g→h takes two forms (see AppendixB) whether cluster g is empty after removing i or not. If no label swapping induces an increase of the criterion, the vector Z i is not modified. Otherwise, the label swapping with the maximal increase is applied and Z i is changed accordingly. During the process, clusters may disappear, i.e. their cardinality reaches zero. Each time one of these moves is accepted, the model is updated and the corresponding column is removed from the cluster indicator matrix Z. Finally, the algorithm stops if a complete pass over the vertices did not lead to any increase of the ICL ex criterion. Thus, the algorithm, automatically infers the number of clusters while clustering the vertices of the network. Starting with an over-segmented initial solution our approach simplifies the model until a local maximum is reached.

Complexity

In order to set up such an algorithm, it is sufficient to know how to compute the changes in the ICL ex criterion induced by the possible swap movements (from cluster g to cluster h) for a given node i, the others being kept fixed. Such changes can be computed efficiently (see AppendixB for details) and the complexity of finding the best swap movement for a node is in average O(l+K 2 ), where l is the average number of edges per node. Such complexity can be achieved, since good approximations of the logarithm of the gamma function are available with constant running time. The greedy algorithm has therefore a total complexity of O(N (l + K 2 up ) + L), since a swap movement cost is O(l + K 2 ); the initialization of the edges counters (η kl , ζ kl ) cost is L (the total number of edges in the graph) and several complete passes over the set of nodes will be performed (typically less than 10). Eventually, this can be simplified in O(N K 2 up + L) and compared to the complexity of O(LK 3 up ) achieved using a variational algorithm and a model selection criterion as in [START_REF] Daudin | A mixture model for random graph[END_REF][START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF]. Indeed, contrary to our approach which estimates the number of clusters in a single run, while clustering the nodes, these approaches are run multiple times for various values of K and K * is chosen such that the corresponding model selection criterion is maximized. Since each run costs O(LK 2 ), the overall complexity is O(LK 3 up ).

Initialization and restarts

Several solutions are possible for initializing the algorithm, a simple choice consisting in sampling random partitions while a more relevant though expensive starting point can be obtained with the k-means algorithm. One possible trade-off in terms of computational burden is to use only few iterations of k-means. We used the latter method in all the experiments that we carried out. Moreover, since our method is only guaranteed to reach a local optima, a common strategy is to run the optimization algorithm with multiple initializations and to keep the best one according to the ICL ex criterion.

Hierarchical clustering

Eventually, in a final step, it is possible to check that merge movements between clusters do not induce any increase of the objective function. This can be done with a greedy hierarchical algorithm which costs O(K 3 ) (see details in AppendixC). Since the labels swap algorithm usually greatly reduces the number of clusters (K << K up ), the computational cost of this last step is low.

Such a scheme leads to a fast algorithm: sparse networks take about 15 seconds for N = 500 nodes, and about five minutes for N = 5000 with a naive Matlab implementation.

Experiments on synthetic data

To assess the greedy optimization method, a simulation study was performed and the proposed solution was compared with available implementations of algorithms for SBM inference:

• vbmod, [START_REF] Hofman | A bayesian approach to network modularity[END_REF], a variational-based approach dedicated to the search of community structures, implemented in Matlab and C. The random graph model they considered can be seen as a constrained SBM where all terms on the diagonal of the connectivity matrix Π are set to a unique parameter λ and off-diagonal terms to another parameter ,

• mixer, [START_REF] Daudin | A mixture model for random graph[END_REF], another variational approach but one which can deal with all types of SBM models (not only communities structures) implemented in R and C,

• colsbm, [START_REF] Mc Daid | Improved bayesian inference for the stochastic block model with application to large networks[END_REF], a collapsed Gibbs sampler for SBM in C. The last version of the code is used in this experimental section. It involves an additional move type compared to the algorithm described in the associated publication. This move was found to greatly enhance the results.

Our goal here is to evaluate the ability of the different solutions to recover a simulated clustering without knowing the number of clusters. Only a reasonable upper bound K up on K will be provided to the algorithms when needed. We recall that the variational methods optimize a bound for various values of K and select K * such that a model selection criterion is maximized: ICL for mixer and ILvB for vbmod. Conversely, the collapsed Gibbs sampler automatically provides an estimate of K, since the posterior of K is made available.

The performances are assessed in terms of normalized mutual information (see for instance [START_REF] Vinh | Information theoretic measures for clusterings comparison: Variants, properties, normalization and correction for chance[END_REF]) between the estimated cluster membership matrix Z e and the simulated one Z s . The mutual information I(Z e , Z s ) between two partitions is to this end defined by:

I(Z e , Z s ) = K k,l p kl log p kl p e k p s l , (4)
with

p kl = 1 N N i,j Z e ik Z s jl , p e k = 1 N N i=1 Z e ik , p s l = 1 N N i=1 Z s il .
The measure I(Z e , Z s ) describes how much is learnt about the true partition if the estimated one is known, and vice versa. The mutual information is not an ideal similarity measure when the two partitions have a different number of clusters and it is therefore preferable to use a normalized version of the mutual information such as:

N I(Z e , Z s ) = I(Z e , Z s ) max (H(Z e ), H(Z s )) , (5)
with

H(Z) = -K k=1 p k log(p k ) and p k = 1 N N i Z ik .
The performances are evaluated on simulated clustering problems of varying complexity and with different settings, in order to give insights about the influence of the number K of clusters, of the number of vertices N and of the type of connectivity matrix Π.

Setting 1: small scale community structures

The first setting is a classical community simulation with N = 100 vertices and K = 5 clusters. The cluster proportions are set to α = (1/5, 1/5, 1/5, 1/5, 1/5) and the connectivity matrix takes a diagonal form with off-diagonal elements equal to 0.01: Π kl = 0.01, ∀k = l and diagonal elements given by Π kk = β, ∀k. β is a complexity tuning parameter which ranges from 0.45 to 0.01. When β reaches 0.01, the model is not identifiable (the connectivity matrix is constant) and the true cluster memberships cannot be recovered. The set of simulated problems is therefore of varying complexity: from problems with a clear structure when β = 0.45 to problems without any structure when β = 0.01. The experiments are performed twenty times for each value of β and the average of the normalized mutual information over these twenty simulated graphs is depicted in Figure 1 (left) for all the algorithms. In order to produce results as comparable as possible, the parameters of the different algorithms were set as follows: vbmod, mixer and greedy icl were all started ten times and for each method the best run was selected according to the corresponding model selection criterion. The variational methods were run with K between 2 and 20 and the best clustering kept as a final result. For greedy icl, the parameters of the prior η 0 , ζ 0 and n 0 k were set to 1 and K up fixed to twenty. Finally the collapsed Gibbs sampler was run for 250 000 iterations (more than twice the default value).

The results illustrated in Figure 1 show that greedy icl outperforms the other methods for complex problems, i.e. low values of β. The simulated clustering is recovered until β reaches 0.25. Above this value the different algorithms perform identically, but beyond this limit the results of greedy icl are a little bit better.

During the transition greedy icl gets slightly better results than the other algorithms, it is followed by colsbm and vbmod which give close results and mixer that deviates earlier from the planted clustering.

Setting 2: small scale community with a hub cluster

The second setting aims at exploring the performances of the methods when the latent structure does not correspond only to communities. To this end, graphs were generated using the stochastic block model with affiliation probability matrix Π of the form as in [START_REF] Latouche | Variational bayesian inference and complexity control for stochastic block models[END_REF]:

Π =            β β . . . . . . β β β . . . β β . . . β . . . β β . . . . . . β           
.

The clusters correspond therefore to communities, except one cluster of hubs which connects with probability β to all other clusters. Graphs with N = 100 vertices, K = 5 clusters and α = (1/5, 1/5, 1/5, 1/5, 1/5) were generated using this connection pattern. The parameter was set to 0.01 and β ranged as previously from 0.45 to 0.01. Eventually, the other simulation parameters did not change. The results are shown in Figure 1 (right).

As expected, the vbmod algorithm, which looks only for communities, is strongly affected by this change of setting and systematically misses the hub cluster. For the remaining methods, the best results are achieved by greedy icl which still uncovers the planted clustering when β > 0.25, whereas mixer starts to drop at β equals 0.4.

The collapsed Gibbs sampler achieves also good results in this setting, very close to those of greedy icl and out-performs mixer.

Setting 3: medium scale community structures

The third setting is a classical community simulation but with more nodes and clusters, in order to study the effect of these two parameters. Thus, the number of vertices was set to N = 500 and the number of clusters to K = 10. The cluster proportions were defined as α = (1/10, . . . , 1/10) and all the other parameters kept the same value as previously. For this third experiment, the results presented in Figure 2 are identical for greedy icl, colsbm and mixer. The vbmod algorithm seems to be more affected by the dimensionality of the problem, and did not recover exactly the true clusters when β is under 0.2. The results obtained by the different algorithms in this setting are better than those obtained previously. This can easily be explained by the increase in the number of nodes per cluster. The transitions between high and low values of the normalized mutual information were also sharper than in the previous experiments, for the same reasons.

Setting 4: large scale problem with complex structure

The final setting involves larger graphs with N = 10 000 vertices. The planted structure is also not a purely community pattern. Some interactions between clusters are activated randomly using a Bernoulli distribution as described by the following generative model:

Π kl =      ZU + (1 -Z) , if k = l U , if k = l (6)
with Z ∼ B(0.1), U ∼ U(0.45) and = 0.01. The size of the problem and the complex nature of the underlying structure, let only two algorithms able to deal with these graphs namely greedy icl and colsbm, since mixer cannot handle such large graphs and vbmod deals only with community structure. Both were used to cluster 20 simulated graphs generated using this scheme. The greedy algorithm was started using K up = 100 and same setting as previously for the prior distributions.

The results presented as boxplots in Figure 3 give a clear advantage to greedy icl over the collapsed sampler. Thus, greedy icl achieves an average normalized mutual information of 0.88 whereas colsbm reaches only 0.67. In fact, the greedy solution ended with around 80 clusters for all the simulations whereas the Gibbs sampler gives more than 240 clusters in average and therefore produces highly over segmented partitions of the graphs. To summarize the results we obtained in all the experiments we carried out, it appears that greedy icl compares favourably with the other existing solutions for SBM, in all the settings. The results obtained in complex setting, i.e. large graphs and a complex underlying structure (Setting 4) are particularly encouraging since greedy icl clearly outperforms the collapsed Gibbs sampler.

Real dataset: communities of blogs

The proposed algorithm was finally tested on a real network where vertices correspond to blogs and edges to known hyperlinks between the blogs. All the blogs considered are related to a common topic, i.e. illustrations and comics.

The network was built using a community extraction procedure [START_REF] Côme | The noise cluster model, a greedy solution to the network community extraction problem[END_REF] which starts from known seeds and expands them to find a dense core of nodes surrounding them.

It is made of 1360 blogs linked by 33 805 edges. The data set is expected to present specific patterns, namely communities, where two blogs of the same community more likely to be connected that nodes of different communities. To test this hypothesis we used the greedy ICL algorithm and did a qualitative comparison of the results with those obtained with the community discovery method of Blondel et al. [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF].

Starting with K up = 100 clusters, greedy ICL found K = 37 clusters. The corresponding clusters are illustrated in Figure 4 which is an image of the adjacency matrix with rows/columns sorted by cluster number. Thus, it appears that the clusters found correspond in their vast majority to small sub-communities. These sub-communities all correspond to known groups. For instance a group of blogs of illustrators for Disney was found. Other examples include clusters of blogs of students who went to the same illustration school such as the ECMA school of Angouleme or the "Gobelins École de l'image". However, some clusters have more complex connectivity structures and are made of hubs which highly connect to blogs of different clusters. They correspond to blogs of famous writers such as Boulet.

To give a qualitative idea of the interest of the found clustering, we also give the results obtained by the community discovery algorithm of Blondel et al. [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] in Figure 5. With this approach only 8 clusters are found, corresponding all to sub-communities. Clusters of hubs could not be recovered. The major difference between the number of clusters estimated by the two methods may be explained by two facts. Firstly, modularity is known to be prone to a resolution limit problem [START_REF] Fortunato | Resolution limit in community detection[END_REF] which prevents such a solution to extract small scale structures. This explains why the small sub-community extracted by greedy icl are not recovered using the modularity. For the time being, the behaviour of the ICL ex criterion with respect to the resolution limit problem is not clear and will deserve further investigations.

However, we notice that on this dataset finer structures than those obtained using modularity are recovered. Secondly, the difference in the way the two criteria use degree correction or not [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF] can also explain the disparity in the number of clusters.

While modularity is a degree-corrected criterion which downscales the weights of the edges between highly connected vertices, the ICL ex criterion for the basic stochastic block model used here is not. Using a degree correction or not is a modelling choice which deserves to be validated and investigated; however, it seems that even without degree correction the results obtained by greedy icl are meaningful, the hub clusters being interesting per se.

Conclusion

In this paper, we showed how an analytical expression of the integrated complete data log likelihood could be derived using conjugate priors for the model parameters, and that no asymptotic approximations were required. We then proposed a greedy optimization algorithm to maximize this exact quantity. Starting from an over segmented partition, the approach simplifies the model, while clustering the vertices, until a local maximum is reached. This greedy algorithm has a competitive complexity and may handle networks with tens of thousands of vertices. We illustrated on simulated data that the method improves over existing graph clustering algorithms, both in terms of model selection and clustering of the vertices. A quali- tative comparison between methods was also carried out on an original network we built from blogs related to illustration, comics, and animations.

We emphasize that the methodology we considered can be adapted to other mixture models. In particular, we will investigate the case of the degree corrected stochastic block model which have been shown to give promising results on real data.

well as Z test . Z describes the current partition of the vertices in the network, while Z test represents the partition after applying the swap g → h:

         Z test j = Z j , ∀j = i Z test ik = Z ik = 0, ∀k = g, h while          Z test ig = 0, Z ig = 1
Z test ih = 1, Z ih = 0 Thus ∆ g→h = ICL ex (Z test , K test ) -ICL ex (Z, K).

Note that ∆ g→h takes two forms whether cluster g is empty after removing i or not.

In the later scenario, the model dimensionality changes (K test = K -1) and this must be taken into account to evaluate the possible increase induced by the swap movement.

AppendixB. kl the changes in edges counter η kl induced by the label swap:

δ (i) kl = 1 {k=h} N j =i Z jl X ij + 1 {l=h} N j =i Z jk X ji -1 {k=g} N j =i Z jl X ij -1 {l=g} N j =i Z jk X ji .
Moreover, ρ

kl is defined in the following:

ρ (i) kl = 1 {k=h} -1 {k=g} (n l -n 0 l -Z il ) + 1 {l=h} -1 {l=g} (n k -n 0 k -Z ik ) -δ (i)
kl .

These update quantities can be computed in O(l i ) with l i the degree of i (total number of edges from and to i). Therefore the complexity of finding the best swap movement for a node is O(l + K 2 ), l for computing the δ

(i)
kl and K 2 to compute the ∆ swap with all the possible h labels and keep the best one.

AppendixB.2. Case 2 :

i Z test ig = 0, cluster g disappear

In this case the dimensionality of n 0 changes and we will denote by n 0 * = (n 0 , . . . , n 0 ) the corresponding vector of size K -1.

kl is defined in the following:

Algorithm 1 :

 1 Greedy ICL Set K = K up ; stop = 0 ; Initialize the N × K up matrix Z ; Compute η,ζ,n ; while stop = 1 do V = {1, . . . , N } ; stop = 1 ; while V not empty do Select a node i randomly in V ; Remove i from V ; If i is in cluster g, compute all terms ∆ g→h , ∀h = g ; if at least one ∆ g→h is positive then stop = 0 ; Find h such that ∆ g→h is maximum ; Swap labels of i : Z ig = 0 and Z ih = 1 ; if g is empty then Remove column g in Z ; Set K = K -1 ; end Update rows and columns (g, h) of the matrices η and ζ ; Update the components g and h of vector n; end end end Result: (Z, K)

Figure 1 :

 1 Figure 1: Mean of mutual information between estimated and true cluster membership matrices using 20 simulated graphs for each value of β in {0.45, 0.43, . . . , 0.03, 0.01}, and with N = 100, K = 5, = 0.01 for the different algorithms greedy icl, vbmod, colsbm and mixer.

Figure 2 :

 2 Figure 2: Mean of mutual information between estimated and true cluster membership matrices using 20 simulated graphs for each value of β in {0.45, 0.43, . . . , 0.03, 0.01}, and with N = 500, K = 10, = 0.01 for the different algorithms greedy icl, vbmod, colsbm and mixer.

Figure 3 :

 3 Figure 3: Mean of the mutual information between estimated and true cluster membership matrices using 20 simulated graphs with N = 10000 and K = 50.

Figure 4 :

 4 Figure 4: Adjacency matrix of the of blogs, the rows/columns are sorted by cluster number with clusters found by the greedy ICL algorithm. The cluster boundaries are depicted with white lines.

5 :

 5 Adjacency matrix of the network of blogs, the rows/columns are sorted by cluster number with clusters found by modularity optimization. The clusters boundaries are depicted with white lines.

∆

  g→h = log C(n 0 ) C(n) C(n test ) C(n 0 * ) + (k,l) =g k=h or l=h log B(η kl + δ(i) kl , ζ kl + ρ (i) kl ) B(η kl , ζ kl ) + k=g or l=g log B(η 0 kl , ζ 0 kl ) B(η kl , ζ kl ) = log n h n 0 Γ ((K -1)n 0 ) Γ(Kn 0 + N ) Γ(K n 0 )Γ((K -1)n 0 + N ) + (k,l) =g k=h or l=h log B(η kl + δ (i) kl , ζ kl + ρ (i) kl ) B(η kl , ζ kl ) + k=g or l=g log B(η 0 kl , ζ 0 kl ) B(η kl , ζ kl )the complexity in this case is the same as previously i.e. O(l + K 2 ).AppendixC. Change in ICL induced by a merge movement∆ g∪h = log C(n 0 ) kl , ζ kl ) + k=g or l=g log B(η 0 kl , ζ 0 kl ) B(η kl , ζ kl ) = log Γ(n 0 ) Γ ((K -1)n 0 ) Γ(Kn 0 + N ) Γ(K n 0 )Γ((K -1)n 0 + N ) Γ(n h + n g -n 0 ) Γ(n g )Γ(n h ) + (k,l) =g k=h or l=h log B(η kl + δ (i) kl , ζ kl + ρ (i) kl ) B(η kl , ζ kl ) + k=g or l=g log B(η 0 kl , ζ 0 kl ) B(η kl , ζ kl ) with δ (i)kl the changes in edges counter η kl induced by the merge: δ(i) kl = 1 {k=h} (η gl -η 0 gl ) + 1 {l=h} (η kg -η 0 kg ) + 1 {k=h and l=h} (η gg -η 0 gg ). (C.1) Moreover, ρ

  ρ

= 1

 1 {k=h} (ζ gl -ζ 0 gl ) + 1 {l=h} (ζ kg -ζ 0 kg ) + 1 {k=h and l=h} (ζ gg -ζ 0 gg ). (C.2)

  1. Case 1 : i Z test ig > 0. Cluster g not empty after removing i

	with δ	(i)
	∆ g→h = log	C(n test ) C(n)	+	K k,l	log	B(η test kl , ζ test kl ) B(η kl , ζ kl )
		= log	Γ(n test g )Γ(n test h ) Γ(n g )Γ(n h )	+	K l=1 k∈{g,h}	log	B(η test kl , ζ test kl ) B(η kl , ζ kl )
		+	k / ∈{g,h} l∈{g,h}	log			B(η test kl , ζ test kl ) B(η kl , ζ kl )
		= log	Γ(n g -1)Γ(n h + 1) Γ(n g )Γ(n h )	+	K l=1 k∈{g,h}	log	B(η kl + δ kl , ζ kl + ρ (i) (i) kl ) B(η kl , ζ kl )
		+	k / ∈{g,h} l∈{g,h}	log			B(η kl + δ kl , ζ kl + ρ (i) B(η kl , ζ kl )	(i) kl )
		= log	n h n g -1	+	K l=1 k∈{g,h}	log	B(η kl + δ kl , ζ kl + ρ (i) B(η kl , ζ kl )	(i) kl )
		+	k / ∈{g,h} l∈{g,h}	log			B(η kl + δ B(η kl , ζ kl ) (i)	(i) kl )	,

kl , ζ kl + ρ

APPENDIX

AppendixA. Marginal distributions Proposition AppendixA.1. marginal distribution p(Z|K) is given by:

where the components of the vector n are n k = n 0 k + N i=1 Z ik , for all k in {1, . . . , K} and the function

and we denote n the vector with components n k = n 0 k + N i=1 Z ik for all k in {1, . . . , K}. Thus

Proposition AppendixA.2. The marginal distribution p(X|Z, K) is given by:

,

Proof.

where

.

Proposition AppendixA.3. Using factorized and conjugate prior distributions over the model parameters, the integrated complete data log likelihood is given by:

where

• the components of the vector n are

AppendixB. Change in ICL induced by a swap movement i : g → h At each step of the greedy ICL algorithm, a single node i is considered. If i is currently in cluster g, the method tests every possible label swapping g → h, that is removing i from cluster g and assigning it to a cluster h = g. The corresponding changes in the ICL ex criterion are denoted ∆ g→h . In order to derive the calculation of each term ∆ g→h , for all h = g, we consider two cluster indicator matrices Z as