Clique-based Method for Social Network Clustering
Guang Ouyang Dipak K. Dey Panpan Zhang




 
       
       
       



 Clique-based Method for Social Network Clustering
     Guang Ouyang, 
       
       
       
       Dipak K. Dey, 
       
       
       
       
       Panpan Zhang
    
=====================================================================================================================




In this article, we develop a clique-based method for social network clustering. We introduce a new index to evaluate the quality of clustering results, and propose an efficient algorithm based on recursive bipartition to maximize an objective function of the proposed index. The optimization problem is NP-hard, so we approximate the semi-optimal solution via an implicitly restarted Lanczos method. One of the advantages of our algorithm is that the proposed index of each community in the clustering result is guaranteed to be higher than some predetermined threshold, p, which is completely controlled by users. We also account for the situation that p is unknown. A statistical procedure of controlling both under-clustering and over-clustering errors simultaneously is carried out to select localized threshold for each subnetwork, such that the community detection accuracy is optimized. Accordingly, we propose a localized clustering algorithm based on binary tree structure. Finally, we exploit the stochastic blockmodels to conduct simulation studies and demonstrate the accuracy and efficiency of our algorithms, both numerically and graphically.



Keywords: Clique-score index, localized clustering algorithm, modularity, social network, spectral analysis,stochastic block model




§ INTRODUCTION


Networks are proliferating all around us, and they appear in different forms, such as (hardwired) electrical grids or (virtual) social relationships. Networked systems spread in various scientific and applied disciplines, for example, the Internet, the World Wide Web, metabolic networks, neural networks, food webs and social networks, etc. In this paper, we place our focus on social networks. Social network analysis is a branch of the social science which is an academic discipline studying a society and the behavior of entities therein. In sociometric or other quantitative studies, social networks are usually modeled by graph structures consisting of a set of nodes or vertices connected by directed arcs or undirected edges. More specifically, the actors in a social network are represented by nodes. For each pair of nodes, if there appears some pattern of ties, they will be connected by an edge. Edges can be directed or undirected, depending on the feature and interpretation of the network. There are a variety of social networks, such as the Facebook <cit.>, marriage networks <cit.> and smart business networks <cit.>.

Due to solid theoretical foundation of statistics, statistical methodologies have been a significant input to the studies of social networks. A plethora of statistical methods have been established and developed to uncover relational structure of social networks, and dyadic ties between actors therein. In this paper, we focus on a significant research topic in social networks—clustering. More precisely, our goal is to develop some statistical methods for accurately clustering actors in a social network into mutually exclusive communities. This process, in many literatures, is called “community detection.” 

Roughly speaking, the fundamental principle underlying social network clustering is that a group of actors who are excessively connected are more likely to form a community. Technically speaking, the formation of a community requires that the connections of actors within the community are significantly higher than the connections between actors from different communities. From the sociological point of view, the occurrence of high connection density between actors in a community is usually due to some kind of homology or homogeneity of actors. For instance, consider a friendship network on Facebook. Intrinsically, students from the same department of a college are more likely form a friendship community, as they have a very high probability to know and friend each other. The homogeneities are reflected in the location parameter (i.e., college) and the academic parameter (i.e., department). On the other hand, students with different education, social or geographic background are much less likely to be connected. 

Past research on social network clustering can be summarized into two categories. One approach is to propose a (parametric or nonparametric) probabilistic graphical model (PGM) which characterizes the community structures of a social network. Pioneering work of such model-based approach were the p_1 model <cit.> and the stochastic blockmodel (SBM) <cit.>. Several successful models were proposed for community detection in the last two decades, including an extension of the SBM <cit.>, a latent position model <cit.>, a latent position cluster model <cit.>, and a mixed membership SBM in <cit.>, etc. We refer the interested readers to <cit.> for a complete and comprehensive review of PGMs for social network clustering. Another approach is to consider a metric that can be used to quantitatively evaluate the quality of social network clustering. The task is to specify an objective function based on the proposed metric, and the ultimate goal is to design an efficient algorithm to optimize the objective function over all potential clustering strategies. In this article, we call this class of approaches metric-based methods.  Precursory work in this direction traced back to <cit.>, in which a measure called cluster coefficient was proposed to evaluate mutual acquaintance between actors in a social network. We refer the interested readers to <cit.> and <cit.> for extensive discussions about cluster coefficient. While cluster coefficient is the metric that is often used for random graphs or dynamic network models, we place our focus on clustering problems of static networks in this manuscript. Representative and popular methods for static network clustering are the spectral clustering method <cit.>, the normalized cut approach <cit.>, and the modularity maximization method <cit.>, etc.

The contribution of this paper is as follows:

	
  * We propose a novel clustering method based on graph cliques such that the clique score (a new index defined to evaluate clustering outcomes) of each identified community is higher than some threshold p.
	
  * When threshold p is unknown, we develop a systematic strategy to select localized thresholds p (for each subnetwork which requires further subdivision) by controlling over-clustering and under-clustering simultaneously.
	
  * The two algorithms proposed in this manuscript are efficient, easy to implement, and provide consistently reliable clustering results.


We would like to point out that even though undirected binary social networks are considered in this manuscript for the sake of interpretation, our method can be extended to directed or weighted networks in a similar manner, done mutatis mutandis. The rest of this manuscript is organized as follows. In Section <ref>, we introduce some notations that will be used throughout the paper. In Section <ref>, we briefly review two classical clustering methods. Section <ref> is divided into two subsections. In Section <ref>, we propose a new metric to quantitatively measure the quality of clustering results, and then establish an objective function based on it, followed by an algorithm which is designed upon the idea of recursive bipartition to optimize the objective function. In Section <ref>, we consider the situation that the global parameter p is not predetermined, and develop a method to compute localized parameter p so as to maximize the detection accuracy of our algorithm. We further propose a localized algorithm (corresponding to localized parameter p), with some modifications on the algorithm developed in Section <ref>. We present several simulation examples in Section <ref> to demonstrate the efficiency of our algorithms and accuracy and consistency of clustering results. In addition, we numerically compare the performance of the proposed algorithm and the traditional modularity maximization algorithm. Lastly, we address some concluding remarks and potential future studies in Section <ref>.



§ NOTATIONS


In this section, we introduce the notations that will be used throughout the manuscript. Let G = G(V, ) be an undirected graph to model a social network of  |V| = n actors, where  is an n × n adjacency matrix such that

    = (a_ij)_n × n,

in which 

    a_ij = 
    1,        ,
    
    
    0,        ,

for i, j ∈ V. Let  be an n × n diagonal matrix such that 

    = [ (1)   0   ⋯   0;   0 (2)   ⋯   0;   ⋮   ⋮   ⋱   ⋮;   0   0   ⋯ (n) ],

where (i) represents the degree of node i ∈ V. Consider the matrix  defined as

    =  - .

The normalized Laplacian matrix of G <cit.> is then given by

    = ^-1/2^-1/2 =  - ^-1/2^-1/2,

where  denotes the identity matrix of rank n.

Two graph invariants on which many classical network clustering methods depend on are volume and cut. The volume of a graph G, denoted (G), is the total number of degrees of the nodes in G, i.e., 

    (G) = ∑_i = 1^n(i).
 
The cut, on the other hand, is defined on subgraphs of G. Let G_1 and G_2 be two disjoint subgraphs of G, the cut of G_1 and G_2, denoted (G_1, G_2), is the number of edges linking G_1 and G_2, i.e.,

    (G_1, G_2) = ∑_i ∈ G_1, j ∈ G_2 a_ij.
 



§ PREREQUISITES


In this section, we give brief reviews of the spectral clustering method with a concentration on the normalized cut and the modularity maximization approach as prerequisites. 


 §.§ Spectral network clustering


We first look at the spectral network clustering algorithm developed in <cit.> and <cit.>, inspired from the spectral graph theory <cit.>. Generally speaking, a measure called the normalized cut was adopted as an objective function to quantify the number of edges across different communities. Suppose that a social network G is split into h communities G_1, G_2, …, G_h. The normalized cut is given by

    (G_1, G_2, …, G_h) = ∑_k = 1^h(G_k, G_k)/(G_k),

where G_k denotes the complement of G_k in G, for k = 1, 2, …, h. Consider an n × h community indicating matrix = (p_ik)_n × h, in which each entry is p_ik = 1/√((G_k)) for node i in community G_k. The optimization problem corresponding to Equation (<ref>) is, in fact, a discrete trace minimization problem, which can be approximated by a standard trace minimization problem with a relaxation of the discreteness condition as follows:

    min_G_1, G_2, …, G_h(^⊤) ≈min_∈ℝ^n × h(^⊤),

subject to ^⊤ =. Let = ^1/2. The optimization problem on the right hand side of Equation (<ref>) is equivalent to

    min_∈ℝ^n × h(^⊤),

subject to ^⊤ =. This is a standard Rayleigh quotient problem, and the solution of  is composed of the first h eigenvectors of ; see <cit.> and <cit.>. Thus, the solution to the right hand side of Equation (<ref>) can be obtained by solving the first h eigenvectors in a generalized eigenvalue system = λ, which can be done via an adaptive algorithm based on bipartition proposed in <cit.>.



 §.§ Modularity maximization


Another popular algorithm for social network clustering, known as the modularity maximization, was first introduced in <cit.>. In essence, the underlying principle is to partition a social network into mutually exclusive communities such that the number of edges across different communities is significantly less than the expectation, whereas the number of edges within each community is significantly greater than the expectation. In <cit.>, a bipartition situation was considered, and the clustering outcome was evaluated by a measure—modularity—defined as

    Q = 1/2 (G)^⊤,

where  is an n × 1 column indication vector such that

    s_i = 1,        ;
    
     s_i = -1,        ,

and = (b_ij)_n × n is an n × n matrix with entires

    b_ij = a_ij - (i)(i)/(G),

for i = 1, 2, …, n. Subdivisions for existing communities are available by repeatedly implementing the proposed modularity algorithm. The decision of whether or not subdividing an existing community G_k of size n^G_k depends on an associated modularity matrix ^G_k = (b^G_k_ij)_n^G_k× n^G_k with entries

    b^G_k _ij = b_ij - δ(i, j) ∑_l ∈ G_k b_il,

where δ(·, ·) denotes the Kronecker delta function. Subdivision for G_k is terminated if the largest eigenvalue of ^G_k is zero. 



§ CLIQUE-BASED CLUSTERING ALGORITHM


In this section, we propose a novel algorithm for social network clustering. More specifically, we define a new measure called the p-clique index to quantitatively evaluate the quality of clustering outcome. The section is divided into two subsections. We introduce a clique-based clustering algorithm in Section <ref>, and then extend it to a localized clique-based clustering algorithm in Section <ref>. The localized algorithm enables us to update the threshold p for each subnetwork so as to well control the potential over-clustering or under-clustering problems. We will discuss the details in the sequel.



 §.§ Clique-based clustering algorithm


In this section, We give a clique-based clustering algorithm, inspired from the modularity maximization algorithm proposed in <cit.>. In graph theory, a clique is defined as a complete graph on a set of nodes, i.e., each pair of nodes is connected by an edge. A clique is called maximal if it cannot be extended to any larger-size clique by including any adjacent node. The ideal clustering outcome is that each community in a network is a maximal clique. In reality, this is hard to achieve. For many real world social networks, it is even difficult to guarantee that each community forms a clique. Therefore, an appropriate measure to assess the degree of connectivity of a community is needed. The next graph invariant measure can be used to gauge the internal link density of each community in a network.
[clique score]
	The clique score of a community (cluster, subnetwork) is the ratio of the number of observed ties to the number of edges in the clique over the same number of nodes.

For example, a community consisting of 10 nodes and 18 internal links has clique score 18/102 = 0.4. It is obvious that a clique always has clique score one. Our goal is to develop an adaptive algorithm to maximize the overall clique scores (e.g., weighted average) of all communities in a network, subject to the clique score of each community exceeding some predetermined threshold 0 ≤ p ≤ 1. The condition of “exceedance” in our algorithm is essential as it guarantees that none of the communities in our clustering outcome performs extremely bad or does not achieve the minimum standard. In addition to this, the choice of p is flexible, depending on the users' needs or the realistic features of communities in a social network. In Section <ref>, we will discuss how to choose an appropriate value of p to optimize the performance of the proposed algorithm when no prior information of p is available. In the next definition, we bridge the gap between clique graph and p.
[p-clique]
	
	A p-clique is a random graph of a set of nodes, of which each pair is connected by an edge independently with probability p, for 0 ≤ p ≤ 1.

The p-clique defined in this manuscript is not novel, and structurally it is equivalent to the  model proposed by <cit.>. The definition of p-clique (c.f. Definition <ref>) adopts an alternative interpretation of the  model given by <cit.>. Since the nodes are connected independently, there is no structure of communities or clusters in p-cliques theoretically. Hence, p-cliques appear to be a proper benchmark model in our study. The expected number of edges in a p-clique on n nodes is np(1 - p).  

Suppose that a social network G of size n is clustered into h communities, G_1, G_2, …, G_h respectively with community sizes n_1, n_2, …, n_h such that ∑_k = 1^h n_k = n. Our task is to search for a clustering rule such that the total degree of nodes in each community is significantly larger than the expected number of edges of a p-clique of the same size, whereas the total number of links across different communities is minimized. Based on this idea, we propose a measure called the p-clique index as follows, and our ultimate goal is to design an adaptive algorithm for an optimal clustering rule where the p-clique index is maximized.
[p-clique index]
	
	Let G be a social network that consists of n nodes, and C = [G_1, G_2, …, G_h] be a clustering rule which divides G into h communities. The p-clique index is given by
	
    D(, p)    = 1/n(n - 1)(∑_k = 1^h((G_k) - p n_k (n_k - 1)) . 
           - . ∑_1 ≤ k ≠ l ≤ h( Cut(G_k, G_l) - p n_k n_l ) ) 
    		
       = 1/n(n - 1)∑_1 ≤ i ≠ j ≤ n((a_ij - p)δ(c_i, c_j) + (p - a_ij)(1 - δ(c_i, c_j))),

	where = (c_1, c_2, …, c_n)^⊤ is the membership indication vector for nodes.

The essence of p-clique index is to reward connected nodes in the same community (with 1 - p) and disconnected nodes from different communities (with p), but penalize connected nodes from different communities (with -(1 - p)) and disconnected nodes in the same community (with -p). Different from model-based methods, this approach is not to fit data (i.e., observation) to any type of generative model or p-clique structures. Rather, our goal is to determine a clustering rule, in which each community has a higher clique score than a predetermined threshold p.

The algorithm for our clique-based clustering approach is based on the hierarchical clustering algorithm developed for modularity maximization in <cit.>. To begin with, we consider bipartition, i.e., clustering a social network into two communities 1 and 2. Define an alternative membership indication vector, = (s_1, s_2, …, s_n)^⊤, as follows:

    s_i = 1,        ,
    
     s_i = -1,        ,

for i = 1, 2, …, n. Then, the p-clique index (c.f. Equation (<ref>)) is equivalent to 

    D(, p) = 1/n(n - 1)∑_1 ≤ i ≠ j ≤ n (a_ij - p) s_i s_j = 1/n(n - 1)^T(p) ,

where the p-clique matrix (p) is given by

    (p) =  - p(_n × n - _n × n),

in which _n × n is an n × n matrix of all ones. In what follows, the network clustering is converted to an optimization problem

    max_{^⊤(p) }.

This is equivalent to

    max_{D(, p) - D(_n × 1, p)}
       = 1/n(n - 1)max_{^⊤(p)  - ^⊤_n × 1(p) _n × 1}
       = 1/n(n - 1)max_{^⊤( (p) - (1/n∑_i, j = 1^n C(p)_ij) _n × n) },

where C(p)_ij denotes the (i, j)th entry of p-clique matrix (p). Since = {-1, 1}^n is dyadic, the optimization problem is NP-hard. To relax the problem, our strategy is to allow  to be any normalized real-valued vector. One solution to Equation (<ref>) is the eigenvector = (v_1, v_2, …, v_n)^⊤ corresponding to the largest eigenvalue of (p) - (1/n∑_i, j = 1^n C(p)_ij) _n × n, denoted by λ_ max. We thus obtain a natural approximation of the solution; that is, we cluster the nodes with respect to the signs of the components of :

    s^ approx_i = 1,        v_i ≥ 0,
    
     s_i^ approx = -1,        v_i < 0,

for i = 1, 2, …, n.

We remark that we would rather work on the optimization problem in Equation (<ref>) than that in Equation (<ref>) analytically, albeit they are mathematically identical. The underlying reason is that Equation (<ref>) provides an instinctive insight of the decision rule of partitioning an existing community. Our rule is that we do not further subdivide the existing cluster if the eigenvalue λ_ max≤ 0 and D(_n × 1, p) ≥ 0. Both of the conditions are needed since λ_ max≤ 0 indicates that D(_n × 1, p) ≥max{max_{D(, p)}, 0 } and D(_n × 1, p) ≥ 0 implies that the link density of the existing community is already higher than p.

We next propose an algorithm based on recursive bipartition for our p-clique approach. In practice, we first determine the p-clique matrix, ^(G)(p) for a given network G of size n. We partition G into two clusters G_1 and G_2 with respect to the signs of the entries of the eigenvector corresponding to the largest eigenvalue of ^(G)(p), and then compute the additional contribution to the p-clique index due to the division, denoted by 

    Δ D(p) = 1/n(n - 1)∑_i, j ∈ G(s_i C^(G)(p)_ij s_j - C^(G)(p)_ij).

We continue to apply the bipartition algorithm respectively to G_1 (according to the associated p-clique matrix ^(G_1)(p), where ^(G_1)(p) is a submatrix extracted directly from ^(G)(p) only for the nodes contained in G_1) and G_2 (according to the associated p-clique matrix ^(G_2)(p) obtained in a similar manner), if at least one of the two regularity criteria is met. We terminate the algorithm until no further partition is needed. The algorithmic procedure is presented in Algorithm <ref>.



1em

1em



To conclude this section, we address some remarks on the proposed algorithm. The essence of Algorithm <ref> is to maximize the additional contributions to the overall p-clique index, whenever a division of an existing cluster is implemented (versus its current state) . As mentioned, the optimization problem is NP-hard, and our strategy is to relax the problem by allowing  to be any real vector with ℓ_2 norm equal to one. Our goal is to determine the eigenvector associated with the largest eigenvalue of the p-clique matrix, and specifically our strategy is to exploit an implicitly restarted Lanczos method developed in <cit.>. It is worthy of noting that our algorithm will continue to be executed if ∑_i,j ∈ G^' C^(G^')_ij(p) < 0 for some G^' even when Δ D(p) ≤ 0, since, under such circumstance, the clique index of G^' is less than p. We take the quality (i.e., having clique score higher than predetermined threshold p) of each resulting community as a higher priority, bearing with non-optimal grade for the overall p-clique index. We view this as a major feature, as well as an advantage, of our algorithm, since the choice of threshold p is absolutely flexible, depending on the users' needs. 

On the other hand, a natural question arising from our algorithm is “ what if there is no specific requirement or prior information for p?” An arbitrary choice of p may lead to over-clustering or under-clustering problems. Having this concern in mind, we propose a modified algorithm based on Algorithm <ref> in Section <ref>.In addition, we carry out a statistical procedure for “reasonable” selection of unspecified p so as to simultaneously minimize the errors of over-clustering and under-clustering.  



 §.§ Localized clustering algorithm


We have demonstrated that the algorithm proposed in Section <ref> is advantageous for social network clustering, as it allows users to control the quality of each resulting community. In many cases, the threshold parameter p is predetermined or preselected, based on users' needs or past research experience. For the cases that p is not specified, it is needed to develop a systematic method to choose an optimal assignment of p. Probably a question that needs to be answered in advance is whether or not there exists such optimum p. An analogous issue occurs for the modularity maximization algorithm, and was discussed in <cit.> and <cit.>, where the conclusion from both was negative for the existence of such overall optimal value of p. In this section, we consider a localized clustering strategy, the core idea of which is to select different values of p for different subnetworks. 

In essence, our strategy for the selection of p is a process of balancing over-clustering and under-clustering. Intrinsically, a large value of p usually results in small sizes of clusters (i.e., over-clustering), whereas a small value of p sometimes fails to guarantee the quality of clusters (i.e., under-clustering). The procedure that balances between over-clustering and under-clustering is analogous to dealing with type I errors and type II errors in statistical hypothesis testings. We borrow these two terminologies in our study. Let us call the error of over-clustering as type I error and the error of under-clustering as type II error. Our goal is determine a distinct value of p for each existing cluster such that both types of errors are well controlled.

Recall the p-clique (i.e., the  model) introduced in Definition <ref>. The model is constructed completely at random without any cluster structure, suggesting that we can exploit it to control type I error, and accordingly determine an upper bound of p. Consider an  graph, ER(n, p_0), where n and p_0 are given. Intuitively, we need to set the threshold parameter p significantly less than p_0 so that there is a small probability to divide ER(n, p_0) into two subgraphs or more by our algorithm. Assume that α is the maximal percentage of the number of nodes to be split from the  graph under tolerance. Our goal turns to find the largest value of p such that at most α n nodes are split off of the majority. Let X denote the observed inter-across link density between the split group (i.e., α n nodes) and the remainder (i.e., (1 - α)n nodes), and F_X be its distribution. An instinctive and reasonable choice of p is thus the αth percentile of F_X. Although the expected number of inter-across edges is (1 - α)np_0 in theory, the inter-across link density is anticipated to be much smaller than the internal density within the large group of (1 - α) n nodes. In this regard, we approximate F_X by a truncated normal distribution, where the top α is curtailed; that is, X is normally distributed with mean

    [X] = p_0 - ϕ(z_α)/1 - α√(p_0(1 - p_0)/(1 - α)n),

and variance

    [X] = p_0 (1 - p_0)/(1 - α)n[1 - z_αϕ(z_α)/1 - α - (ϕ(z_α)/1 - α)^2],

where ϕ(·) denote the density function of the standard normal distribution, and z_α is the (1 - α)th percentile of standard normal. Thus, we obtain an upper bound for p as a function of α, which is given by

    p_ U(α)    = max{0, [X] - z_α√([X])}
       = max{0, p_0 - ξ(α) √(p_0 (1 - p_0)/(1 - α)n)},

where ξ(α) is a constant depending on α only.

The next task is to determine the lower bound for p. Suppose that we have two independent  graphs which are respectively denoted by ER(n_1, p_1) and ER(n_2, p_2), and p_12 is the link density between the two graphs. Let β be the probability that the two graphs are merged. We once again consider normal approximation, and conclude that the two graphs are both statistically significant if we have

    min{p_1, p_2} > p_12 + z_β√(p_12(1 - p_12)/n_1 n_2),

the right-hand side of which in fact forms the lower bound for the parameter p. However, this lower bound is intractable in general as it requires the information of p_12, which is usually unknown for most real-world networks. Therefore, we place focus on the upper bound for p, and our strategy is to set the parameter p as large as possible under tolerance so as to minimize the probability of merging significant subnetworks.

Finally, we propose a new algorithm for unspecified p based on Algorithm <ref> proposed in Section <ref>. The new algorithm is designed for controlling type I and II errors simultaneously. As mentioned, however, it seems to be unreasonable to have a global threshold p for our algorithm. We elaborate the reasonings and demonstrate our solution via the following illustrative example. 

Suppose that we cluster a network G of size n into two subnetworks G_1 of size n_1 and G_2 of size n_2. Initially, we adopt a threshold p^(G) which depends on the observed link density of G. However, if we continue to use p^(G) as the threshold when G_1 or G_2 or both need to be further clustered, we may have to bear with over-clustering or under-clustering risks in the subdivision processes. Alternatively, we suggest to reset the threshold(s) before subdividing G_1 or G_2 or both. In other words, we choose

    p^(G_1) = max{0, p_1 - ξ(α)√(p_1 (1 - p_1)/(1 - α)n_1)}

and

    p^(G_2) = max{0, p_2 - ξ(α)√(p_2 (1 - p_2)/(1 - α)n_2)}

as updated thresholds for G_1 and G_2, respectively. From then on, the threshold parameters are refreshed in this manner for all subnetworks that need to be further subdivided. As for each subnetwork to which a fresh threshold parameter is assigned, we call the new algorithm localized clustering algorithm. 

One of the most effective ways to illustrate this heirachical clustering process is probably to exploit binary tree. We start with a root node that represents the original network which needs to be clustered. At the first level, two child nodes are attached to the root node. The child nodes represent either subnetworks or communities. If a child node is a subnetwork which requires further subdivision, it will carry over two higher-level child nodes and itself turns to an internal node. If a child node is a community which does not need further subdivision, it becomes a terminal node in the tree. We use rectangle for internal nodes and circle for terminal nodes. When the algorithm is terminated, the clustering result is reflected in all the terminal nodes in the binary tree. An example of binary tree is given in Figure <ref>.


	



We are now in the position of defining a new clique index, the local clique index, modified from the p-clique index given in Definition <ref>.
[local clique index]
	Consider a network G consisting of n nodes. Let T be the binary tree that describes a hierarchical clustering procedure on G, and {v} be the collection of all internal nodes of T. The local clique index is given by
	
    LD(T) = 1/n(n - 1)∑_v ∈ T∑_i ≠ j ∈ v((a_ij - p_v)δ_ij^(v) + (p_v - a_ij)(1 - δ_ij^(v))),

	where p_v is the threshold of parameter p for internal node (subnetwork) v.


As shown in Figure <ref>, our new algorithm is again based on recursive bipartition procedures, and the goal is to maximize the score function given in Equation (<ref>). Analogous to Algorithm <ref>, we need to determine the local clique matrix for each v in T, i.e.,

    ^(v) = ^(v) - p_v (_n_v × n_v - _n_v × n_v),

where ^(v) is the updated adjacency matrix for subnetwork v and n_v is the size of v. Besides, we need to compute the additional contribution from each bipartition of v to the overall score LD(T), i.e.,

    Δ LD(v) = 1/n(n - 1)∑_i ≠ j ∈ v((a_ij - p_v)δ_ij^(v) + (p_v - a_ij)(1 - δ_ij^(v))).

The algorithm for hierarchical clustering process is terminated as Δ LD(v) ≤1/n(n - 1)∑_i ≠ j ∈ v (a_ij - p_v) for all existing v in T, and then all remaining internal nodes v turn out to be terminal nodes, i.e., added to the community list. We summarize our strategy in Algorithm <ref>, slightly modified from Algorithm <ref>.



1em

1em





§ SIMULATIONS


In this section, we conduct some simulation studies to evaluate the performance of two clique-based clustering algorithms proposed in the previous sections. We also show that the proposed algorithm outperforms the modularity maximization approach <cit.> for network clustering. The SBM in <cit.> is adopted to generate networks, as they allow us to predetermine the structure of simulated networks, which can be used as the ground truth for comparisons.

To begin with, we briefly review the SBM. Given a network of n nodes which belong to h nonempty communities, let n_1, n_2, …, n_h be respectively the size of each community. For i = 1, 2, …, n, a mapping c: {1, 2, …, n}↦{1, 2, …, h} preserves the membership information for each node labeled with i. An h × h probability matrix  describes the link densities within every community, as well as those between different communities, i.e., = (B_c_(i), c_(j)). Notice that we have to require B_c_(i), c_(j)≤min{B_c_(i), c_(i), B_c_(j), c_(j)} for all i, j, when predefining matrix ; otherwise, the simulated network would be against the clustering structure obtained by our algorithm(s). Suppose that  (the probabilistic structure of network) is specified, we are able to simulate the entries in the adjacency matrix 

    a_ij =  Bernoulli(B_c_(i), c_(j))

for i < j. Suppose that the network is undirected, we have a_ji equal a_ij by symmetry. We also assume that  all the entries on the diagonal equal 0, as loops are not considered in our study.

To quantitatively evaluate the performance of algorithms, we adopt two well-defined robust measures: the normalized mutual information (NMI) and the Adjusted Rand Index (ARI), respectively proposed in <cit.> and in <cit.>. Suppose that T is the ground truth of community structure, and S is the clustering result of an algorithm, the NMI of T and S is given by

    NMI(T, S) = -2 ∑_k = 1^C_T∑_l = 1^C_S N_kllog(N_kl N_··/N_k ·N_· l)/∑_k = 1^C_T N_k ·log(N_k ·/N_··) + ∑_l = 1^C_S N_· llog(N_· l/N_··),

where C_T and C_S are the number of communities for T and S, respectively; (N_kl) is a C_T × C_L confusion matrix, in which N_kl denotes the number of nodes that should be in community k according to the truth, but are mis-clustered into community l according to algorithm S; N_k ·, N_l ·, and N_·· are standard definitions of the sum of the kth row, the sum of the lth column, and the overall sum of the confusion matrix, respectively. We borrow the same notations and give the definition of ARI as follows:


    ARI(T, S) = ∑_k = 1^C_T∑_l = 1^C_SN_kl2 - [∑_k = 1^C_TN_k ·2∑_l = 1^C_SN_l ·2] /N_··2/1/2[∑_k = 1^C_TN_k ·2 + ∑_l = 1^C_SN_l ·2] - [∑_k = 1^C_TN_k ·2∑_l = 1^C_SN_l ·2] /N_··2.


In addition, we propose a measure analogous to NMI.  Once again, assuming that T is the true structure of the simulated network based on SBMs, and S is the analogy from our clustering algorithm, we consider two n × n binary matrices, (M^(T)_ij) and (M^(S)_ij), defined respectively as follows:

    M_ij^(T) = 
    1,        ,
    
     0,        ;

and 

    M_ij^(S) = 
    1,        ,
    
     0,        .

Our measure is defined on a cluster level, and in the form of an h × h matrix. The fundamental principle of our measure is simple; that is, we compute the proportion of mis-clustered nodes, from cluster to cluster. More precisely, for all nodes in different communities k and l (k ≠ l) under T, the error measure is given by

    ϵ_kl = ∑_i, j 
     c^(T)(i) = k, c^(T)(j) = l|M_ij^(T) - M_ij^(S)|/n_k n_l,

where n_k and n_l are the sizes of communities k and l, respectively. The error measure for a single community k is defined analogously, i.e.,

    ϵ_kk = ∑_i, j 
     c^(T)(i) = c^(T)(j) = k|M_ij^(T) - M_ij^(S)|/n_k (n_k - 1).
  
Ultimately, the error matrix is given by

    = (ϵ_kl)_h × h.

We would like to point out that our error matrix measure can be more appealing in some comprehensive analysis, as it shows where (i.e., in which community or communities) exactly a large amount of mis-clustering occurs when the clustering result far deviates from ground truth.



 §.§ Effect of threshold p


In the first experiment, we show that a correct choice of threshold p has a significant impact on the p-clique index clustering algorithm (c.f. Algorithm <ref>). We simulate an SBM consisting of 120 nodes which are clustered into three communities of sizes 100, 10 and 10, respectively. The parameters of in-group and cross-group link densities are summarized in Table <ref>.


The expected clique score of the simulated SBM1 is approximately 0.1597, which can be used as an estimate of the overall link density of the network. Let the error parameter α be equal to 0.025. According to Equation (<ref>), the upper bound for the threshold p is 0.09541. On the other hand, the lower bound for p (c.f. Equation (<ref>)) is max{0.0635, 0.0927} = 0.0927. Suppose that we choose p = 0.09, simulate 100 independent SBM1s, and compute the distance-based error matrix . Albeit the threshold p is just slightly less than the lower bound, we have a large under-clustering error for clusters 2 and 3, i.e., ϵ_23 = ϵ_32 = 0.2218 with standard error SE(ϵ_23) =  SE(ϵ_32) = 0.0376, indicating that the probability of misclassifying nodes from these two clusters is about 22% in average. Suppose that the value of p is increased to 0.11 (even though p = 0.11 is greater than the upper bound), the errors ϵ_23 = ϵ_32 drop dramatically to 0.0614 with standard error SE(ϵ_23) =  SE(ϵ_32) = 0.0172. As ϵ_23 = ϵ_32≈ 6 % is the largest entry in the error matrix  for p = 0.11, it seems that over-clustering does not bring too much trouble in this experiment; the reason is that all predefined in-group densities are significantly larger than cross-group densities. Due to the limit of space, we refer the interested readers to <cit.> for more analogous examples.

Nevertheless, as long as the parameter p is fairly close to the proposed threshold selection interval from the above, the clustering outcomes for this experiment (c.f. Table <ref>) are under satisfactory. We choose the optimal value of p = 0.9527, and depict the result in Figure <ref>, where the three communities are clearly identified and colored by blue, red and green.




 §.§ Global threshold v.s. localized thresholds


The next experiment is designed to compare the performance of the algorithms proposed in this manuscript. We show that the localized algorithm provides more reliable outcomes when the threshold p is unknown. We consider SBM2 with 140 nodes which are clustered into three communities of sizes 100, 20 and 20, respectively. The parameters of in-group and cross-group link densities are summarized in Table <ref>.


The expected clique score of the simulated SBM2 is approximately 0.1546. We again set the error parameter α at 0.025. The associated threshold parameter p according to Equation (<ref>) is 0.0959. As the parameter p is less than the cross-group link density between communities 2 and 3, it seems to be difficult to separate the two clusters with the overall threshold p via Algorithm <ref>. To verify our conjecture, we simulate 100 independent SBM2s, compute the NMI for each replication, and take the average as an estimate. We obtain NMI = 0.8488 with standard error 0.0025. In addition, we compute the proposed block-wise distance-based measure, summarized in Table <ref>.


Although community 1 is successfully identified, we observe that the estimate of error rate between communities 2 and 3 is ϵ_23 = 0.9890 = 98.90 %, which suggests that the BiPartition algorithm with the global threshold p = 0.0959 fails to separate these two communities almost surely. This is also reflected in the clustering result (via Algorithm <ref>) of the simulated SBM (c.f. Table <ref>) given in Figure <ref>. The entire network is divided into two cluster (rather than three), which are colored by red and green.


We remedy this problem by applying the localized clustering algorithm, i.e., Algorithm <ref>. We start with the root node (the original simulated SBM) in the binary tree structure. We use threshold p = 0.0959 to bipartition the root node, and obtain community 1, and subnetwork 2, where the latter requires for further clustering. Subnetwork 2 consists of 40 nodes and the expected clique score is approximately 0.4026. With the same α = 0.025, the updated threshold parameter for subnetwork 2 is 0.2536. We would like to mention that in practice we also compute the threshold parameter for subnetwork 1 (community 1), and obtain 0.1231 for implementing Algorithm <ref>, where we find that no further subdivision is needed. After subdividing subnetwork 2, we continue to compute new thresholds for both subsequent subnetworks, and ultimately find that they form the other two communities in the network.

Another 100 independent SBMs are simulated, and the NMI for each simulated network is computed. The mean estimate is NMI_ local = 0.9677 with standard error 0.0034. We also compute the block-wise error rates, and present them in Table <ref>, where we find a significant improvement of ϵ_23 = ϵ_32 = 0.0019 versus 0.9890. The small error rate implies that the communities 2 and 3 are significantly identifiable according to the updated threshold for subnetwork 2. Analogously, we depict the clustering result of a simulated SBM via Algorithm <ref>, shown in Figure <ref>. The network is successfully divided into three communities as predefined in Table <ref>. The communities are colored by blue, red and green.






 §.§ The propose algorithm v.s. modularity maximization


Notice that the algorithm proposed in this manuscript is inspired from the modularity maximization algorithm. In this section, we would like to compare the proposed algorithm with the modularity approach. In the literature, there is an extensive discussion about the limitations and drawbacks of the modularity-based algorithms. One of the most significant problems of the modularity method is that the maximization of modularity score is not always consistent with the optimized clustering outcome, which was addressed by <cit.>. This limitation is related to under-clustering, as it seems that the modularity maximization algorithm intends to merge small clusters, especially in large-size networks. It is believed that the modularity-based methods undergo over-clustering problems as well. Both of these issues were stressed and discussed in <cit.> and <cit.>. To the best of our knowledge, neither of them has been systematically solved up to date. We refer the readers to <cit.> for extensive discussions about inconsistency of clustering results from the modularity maximization algorithm. On the other hand, both over-clustering and under-clustering concerns are considered and well controlled in our algorithm, which empirically ensures consistency. We present the following three simulation examples to show resistance of the proposed algorithm in this manuscript.

We reconsider SBM1, a simulated network contains two extremely small communities which are extremely loosely linked; SBM2, a simulated network contains two relatively small communities which are fairly loosely linked; and SBM3, a simulated network which is partitioned into communities of approximately even size. The structural parameters are summarized in Table <ref>.



We respectively apply the modularity maximization algorithm and the localized clique-based algorithm to these three networks, and evaluate the clustering results by using the mean estimates of both NMI and ARI, based on 100 independent copies of each simulated network. The results are presented in Table <ref>.



Based on the experimental results, we observe that the modularity maximization algorithm performs extremely poorly for SBM1. As a matter of fact, the modularity algorithm even fails to distinguish the two smaller communities (c.f. clusters 2 and 3). As the sizes of small communities and the link density there-between increase, the modularity maximization algorithm recovers (as reflected in the results for SBM2), but apparently still does not appear as well-performed as our algorithm. For networks comprising similar-size communities (e.g., SBM3), the modularity maximization algorithm seems to outperform our algorithm, but not significantly. Having the randomness of simulated networks in mind, we thus conclude that the proposed algorithm in this manuscript is preferred due to its consistency and reliability.



 §.§ Time complexity

In the era of big data, researchers always concern about the efficiency of algorithms, especially the newly proposed algorithm, when the network size is large. In the last part of this section, we look into this issue numerically by implementing the localized algorithm in Python for several networks of different sizes and structure. Both mean estimates of NMI and ARI are computed. Our testing results are summarized in Table <ref>. We are convinced that the proposed algorithm is able to provide accurate clustering results in a relatively short amount of time.
 



§ CONCLUDING REMARKS


In this section, we address several remarks and discuss some possible future work. To conclude, we propose a new clique-based measure in this manuscript to evaluate the quality of network clustering, and design an algorithm to maximize an objective function based on the proposed measure. The clique score of each community in our clustering result is guaranteed to be higher than some predetermined threshold p. We also consider the situation at which the threshold p is unspecified. We develop an approach which accounts for over-clustering and under-clustering problems simultaneously to select localized p for the network and its subsequent networks. An associated localized algorithm is proposed and discussed.

Studies of networks or network models usually involve big data problems. When network size or parameter space or both are large, it is always challenging to use model-based clustering methods, as many of them depend on accurate but slow Bayesian MCMC algorithms, for example <cit.>. However, the methods proposed in this manuscript attempt to convert clustering analysis to optimization problems. Therefore, many sophisticated machine learning techniques and well-developed approximation methods are ready to use to deal with big data issues.

It is worthy of mentioning that the methods considered in this manuscript also applicable to sparse social network, since the computation of p-clique index is primarily based on the p-clique matrix (p) =  - p(_n × n - _n × n). This matrix is not sparse even if the adjacency matrix  is sparse. Numerical methods, such as the Lanczos method in <cit.>, promise that the eigenvector corresponding to the largest eigenvalue of (p) always can be determined very fast.

The development of the localized clustering algorithm in this paper depends on a strong assumption that the null network is an  model. This may not be true for many real-world networks. It is suggested in <cit.> that many networks around us follow power law, which is not the case for the  graphs. We conjecture that it may be more accurate to consider a scale-free network as null for the studies of social networks possessing power-law property, and we will consider future research in this direction. 













	
	
	
	
	Airoldi
	Airoldi, E.M., Blei, D.M., Fienberg, S.E. and Xing, E.P.
	Mixed membership stochastic blockmodels. 
	Journal of Machine Learning Research, 9, 1981–2014. (2008)
	
	Barabasi
	Barabási, A.-L. and Albert, R. Emergence of scaling in random networks. Science, 286, 509–512. (1999)
	
	Bickle
	Bickel, P.J. and Chen, A. A nonparametric view of network models and Newman-Girvan and other modularities Proceedings of the National Academy of Sciences of the United States of America, 106, 21068–21073. (2009)
	
	Calvetti
	Calvetti, D., Reichel, L. and Sorensen, D.
	An implicitly restarted Lanczos method for large symmetric eigenvalue problems. Electronic Transactions
		on Numerical Analysis, 2, 1–21. (1994)
	
	
	Chung
	Chung, F.R.K. Spectral Graph Theory. American Mathematical Society, Providence, RI. (1997)
	
	Erdos
	Erdös, P. and Rényi, A. On random graphs I. Publicationes Mathematicae, 6, 290–297. (1959)
	
	Fortunato
	Fortunato, S. and Barthélemy, M. Resolution limit in community detection. Proceedings of the National Academy of Sciences of the United States of America, 104, 36–41. (2007)
	
	Fred
	Fred, A. and Jain, A. Robust data clustering, in the Proceedings of IEEE Computer Society Conference on Computer Vision and Pattern Recognition, 2, 128–133. (2003)
	
	Gilbert
	Gilbert, E. N. Random graphs. Annals of Mathematical Statistics, 30, 1141–1144. (1959) 
	
	Goldenberg
	Goldenberg, A., Zheng, A.X., Fienberg, S.E. and Airoldi, E.M.  A survey of statistical network models. Foundations and Trends in Machine Learning, 2, 129–233. (2010)
	
	Handcock
	Handcock, M.S., Raftery, A.E. and Tantrum, J.M. Model-based clustering for social networks. 
	Journal of the Royal Statistical Society, Series A, 170, 301–354. (2007)
	
	Hoff
	Hoff, P.D., Raftery, A.E. and Handcock, M.S. Latent space approaches to social network analysis. 
	Journal of the American Statistical Association, 97, 1090–1098. (2002)
	
	Holland1981
	Holland, P. W. and Leinhardt, S. An exponential family of probability distributions for directed graphs. Journal of the American Statistical Association, 76, 33–50. (1981)
	
	Holland1983
	Holland, P.W., Laskey, K.B. and Leinhardt, S. Stochastic blockmodels: First steps. Social Networks, 5, 109–137. (1983)
	
	Horn
	Horn, R.A. and Johnson, C.R. Matrix Analysis. Cambridge University Press, New York. (1985)
	
	Hubert
	Hubert, L. and Abrabie, P. Comparing partitions. Journal of Classification, 2, 193–218. (1985
	
	Lancichinetti
	Lancichinetti, A. and Fortunato, S. Limits of modularity maximization in community detection. Physical Review E, 84, 066122. (2011)
	
	Newman2001a
	Newman, M.E.J. The structure of scientific collaboration networks. Proceedings of the National Academy of Sciences of the United States of America, 98, 404–409. (2001)
	
	Newman2001b
	Newman, M.E.J., Strogatz, S.H. and Watts, D.J.. Random graphs with arbitrary degree distributions and their applications. Physical Review E, 64, 026118. (2001)
	
	Newman2006
	Newman, M.E.J. Modularity and community structure in networks. Proceedings of the National Academy of Sciences of the United States of America, 103, 8577–8582. (2006)
	
	Ng
	Ng, A.Y., Jordan, M.I. and Weiss, Y. On spectral clustering: Analysis and an algorithm. Advances in Neural Information Processing Systems, 14, 849–856. (2001)
	
	Snijders
	Snijders, T.A.B. and Nowicki, K.
	Estimation and prediction for stochastic blockmodels
	for graphs with latent block structure. Journal of Classification, 14, 75–100. (1997)
	
	Ouyang
	Ouyang, G.
	Social Network Community Detection. Ph.D. dissertation.
	University of Connecticut. (2015)
	
	Pao
	Pao, L.-F. Discovering the dynamics of smart business networks. Computational Management Science, 1, 445–458. (2014)
	
	Pei
	Pei, X., Zhan, X.-X. and Jin, Z. Application of pair approximation method to modeling and analysis of a marriage network. Applied Mathematics and Computation, 294, 280–293. (2017)
	
	Reichardt
	Reichardt, J. and Bornholdt, S. Statistical mechanics of community detection. Physical Review E, 74, 016110. (2006)
	
	Shi
	Shi, J. and Malik, J. Normalized cuts and image segmentation. IEEE Transaction on Pattern Analysis and Machine Intelligence, 22, 888–905. (2000)
	
	Watts
	Watts, D.J. and Strogatz, S.H. Collective dynamics of “small-world” networks. Nature, 440–442. (1998)
	
	Wohlgemuth
	
	Wohlgemuth, J. and Matache, M.T. Small-wold properties of Facebook group networks. Complex Systems, 23, 197–225. (2014) 
	
	


