Clustering nodes in a graph is a useful general technique in data mining of large network data sets. In this context, Newman and Girvan [9] recently proposed an objective function for graph clustering called the Q function which allows automatic selection of the number of clusters. Empirically, higher values of the Q function have been shown to correlate well with good graph clusterings. In this paper we show how optimizing the Q function can be reformulated as a spectral relaxation problem and propose two new spectral clustering algorithms that seek to maximize Q. Experimental results indicate that the new algorithms are efficient and effective at finding both good clusterings and the appropriate number of clusters across a variety of real-world graph data sets. In addition, the spectral algorithms are much faster for large sparse graphs, scaling roughly linearly with the number of nodes n in the graph, compared to O(n 2 ) for previous clustering algorithms using the Q function. Introduction Large complex graphs representing relationships among sets of entities are an increasingly common focus of scientific inquiry. Examples include social networks, Web graphs, telecommunication networks, semantic networks, and biological networks. One of the key questions in understanding such data is "How many communities are there and what are the community memberships"? Algorithms for finding such communities, or automatically grouping nodes in a graph into clusters, have been developed in a variety of different areas, including VLSI design, parallel computing, computer vision, social networks, and more recently in machine learning. Good algorithms for graph clustering hinge on the quality of the objective function being used. A variety of different objective functions and clustering algorithms have been proposed for this problem, ranging from hierarchical clustering to max-flow/min-cut methods to methods based on truncating the eigenspace of a suitably-defined matrix. In recent years, much attention has been paid to spectral clustering algorithms (e.g., [11] , [12] , [14] ) that, explicitly or implicitly, attempt to * The research in this paper was supported by the National Science Foundation under Grant IRI-9703120 as part of the Knowledge Discovery and Dissemination program. SW was also supported by a National Defense Science and Engineering Graduate Fellowship. † Department of Computer Science, University of California, Irvine globally optimize cost functions such as the Normalized Cut measure [12] . The majority of these approaches attempt to balance the size of the clusters while minimizing the interaction between dissimilar nodes. However, for the types of complex heterogeneous networks that arise naturally in many domains, the bias that these approaches have towards clusters of equal size can be seen as a drawback. Furthermore, many of these measures, such as Normalized Cut, can not be used directly for selecting the number of clusters, k, since they increase (or decrease) monotonically as k is varied. Recently, a new approach was developed by Newman and Girvan [9] to overcome limitations of previous measures for measuring community structure. They proposed the "modularity function" Q, which directly measures the quality of a particular clustering of nodes in a graph. It can also be used to automatically select the optimal number of clusters k, by finding the value of k for which Q is maximized, in contrast to most other objective functions used for graph clustering. Let G(V, E, W ) be an undirected graph consisting of the set of nodes V , the set of edges E, and a symmetric weight matrix W ∈ n×n , where n is the number of vertices. The weights w ij = w ji = [W ] ij are positive if there is an edge between vertices v i and v j , and 0 otherwise. The modularity function Q can be defined as where P k is a partition of the vertices into k groups and where A(V , V ) = i∈V ,j∈V w(i, j). Thus, A(V c , V c ) measures the within-cluster sum of weights, A(V c , V ) measures the sum of weights over all edges attached to nodes in cluster c, and A(V, V ) measures the sum of all edge weights in the graph. Considering binary weights for simplicity, the term is the empirical probabilityp c,c that both ends of a randomly selected edge from G lie in cluster c. Similarly, Newman [10] and Newman and Girvan [9] showed across a wide variety of simulated and real-world graphs that larger Q values are correlated with better graph clusterings. In addition, they found that real-world unweighted networks with high community structure generally have Q values within a range from 0.3 to 0.7. Figure 1 shows an example of a simple toy graph with binary weights, where the structure of the graph visually suggests 3 clusters. Also shown are the maximum values of the Q function for different numbers of clusters k, and indeed Q is maximized for k = 3. As pointed out in [10] , if no edges exist that connect vertices across clusters then Q = 1, and conversely if the number of inter-cluster edges is no better than random then Q = 0. We have found empirically that the Q measure works well in practice in terms of both (a) finding good clusterings of nodes in graphs where community structure is evident, and (b) indicating what the appropriate number of clusters k is for such a graph. In this paper we show how Newman's Q measure can be related to the broader family of spectral clustering methods. Specifically: • We show how the problem of maximizing the modularity measure Q can be reformulated as an eigenvector problem involving a matrix we call the "Q-Laplacian." In this manner we link work on graph clustering using the Q measure to relevant work on spectral clustering (e.g., • We use the eigenvector formulation of maximizing Q to derive two new spectral graph clustering algorithms. One of these algorithms directly seeks a global optimum of the Q function. The other algorithm is similar to Newman's agglomerative clustering algorithm [10] , in that it attempts to maximize Q via local iterative improvement. • We compare the new algorithms with Newman's algorithm on different graph data sets and empirically illustrate that: -the spectral approach to maximizing Q produces results that, in terms of cluster quality, are comparable or better than results from Newman's hierarchical algorithm, and -the proposed algorithms are linear per iteration in the number of nodes and edges in the graph, compared to quadratic complexity in the number of nodes for the original algorithm proposed by Newman [10] . 2 Spectral Approaches to Maximizing the Q Function Consider for the moment that the number of clusters, k, is fixed. We use the following strategy to address the problem of finding a partitioning that maximizes Q(P k ) as follows: 1. Reformulate the problem of maximizing Newman's Q function as a discrete quadratic assignment problem. 2. Approximate the resulting assignment problem by relaxing it to a continuous one which can be solved analytically using eigen-decomposition techniques. 3. Compute the top k − 1 eigenvectors of this solution to form a k−1-dimensional embedding of the graph into a Euclidean space. Use "hard-assignment" geometric clustering (the k-means algorithm) on this embedding to generate a clustering P k . Copyright © by SIAM. Unauthorized reproduction of this article is prohibited 275 Below we outline each of these steps. In the section which follows, Section 3, we then describe computational details for two proposed clustering algorithms based on this approach. We assume G is simple, i.e., G contains no self-loops nor parallel edges, and is connected, i.e., there is a path from any vertex to any other vertex. Let D ∈ n×n be a diagonal matrix having d i in the ith diagonal entry and 0 everywhere else, where d i = j w ij . We denote the diagonal matrix derived from an n × n matrix X as diag(X) ∈ n×n where [diag(X)] ij = [X] ij if i = j and 0 otherwise. Finally, let tr(X) = i [X] ii be the trace of matrix X. To simplify notation, we rewrite Q as follows: , and x ic = 0 otherwise, for 1 ≤ c ≤ k. Observe that since each vertex can only be in one cluster, X1 k = 1 n . We can reformulate Q in terms of the assignment matrix X as follows: where d ∈ n×1 such that component d i equals the weighted degree of vertex i and we rewrite A(V, V ) as vol G , the volume of graph G. Thus, where D = dd T . Since for any matrix A and assignment matrix X, tr( Ax c , we can further reduce Q as follows: where W = vol G W . The problem of maximizing Q can then be expressed as: where M ∈ k×k is a diagonal matrix with diagonal entry [M ] ii = |V i |, where |V i | is the number of nodes in cluster V i . Finding an assignment matrix X which maximizes (2.3) is NP-complete. To address this we can attempt to derive a good approximation by relaxing the discreteness constraints that the X ij ∈ {0, 1}, so that instead the X ij ∈ 1 . This transforms the discrete optimization problem into one that is continuous. To find the optimal relaxed X, take the derivative of the following expression with respect to X: where Λ ∈ k×k is the diagonal matrix of Lagrangian multipliers. Setting this equal to 0 and rearranging terms we have: where L Q = D −W and we refer to this diagonal matrix as the "Q-Laplacian". Aside from noting its similarity in form to the standard Laplacian, we observe that this is a standard matrix eigenvalue problem which can be solved using standard eigendecomposition methods. Furthermore, had we normalized the original matrix W so that all rows sum to one and had we also added back in the normalization constant that we took out from equation (2.2) then we would have the following eigenvalue equation: n W where E is a matrix of all ones and W is the matrix W normalized so all rows sum to one. The first term of this equation can be seen as a damping term that ensures that there are edges between all of the nodes of very small weight and the second term is the original weight matrix after scaling and normalization. As n → ∞, the first term will approach 0 much faster than the second term, and hence will play a negligible role in determining the eigenspace of the matrix. Thus, for even moderately large values of n, it seems reasonable that W will provide a close approximation to L Q , which we refer to as the "normalized Q Laplacian." 1 In this paper, we will adopt the simplest method Copyright © by SIAM. Unauthorized reproduction of this article is prohibited 276 for normalizing a matrix so its rows sum to one, namely, to left multiply the matrix W by D −1 . The advantage of using W = D −1 W as an approximation to L Q is that it is easy to compute, it is well studied, especially in relation to Markov chains where it is known as the transition matrix, and it retains sparsity so we can use fast methods for eigendecomposing a sparse matrix. The final step in this framework is to iterate over different values of k, to search for the best clusterings (highest Q(P k ) scores). For each k, we try to find the optimal partitioning, i.e., a "hard-assignment" of the nodes to k clusters, based on clustering the rows of the matrix X. In this section, we propose two new algorithms for clustering graphs that build on insights developed in the previous section. Assume that we are seeking up to a maximum of K clusters and that we have a weight matrix W ∈ nxn . Both of our proposed algorithms below begin by computing the top K − 1 eigenvectors (ignoring the trivial all-ones eigenvector) corresponding to Equation 2.6. Specifically: ] from M using a sparse eigenvector decomposition method such as a variant of the Lanczos method or subspace iteration. In the experimental results in this paper we compute the K − 1 eigenvectors using the Implicitly Restarted Lanczos Method (IRLM) [2] . If one makes the conservative assumption that there are O(K) extra Lanczos steps, then the IRLM has worse-case time complexity of O(mKh + nK 2 h + K 3 h) where m is the number of edges in the graph, and h is the number of iterations required until convergence. For sparse graphs, where m ∼ n, and where K n, we found the IRLM to be extremely fast, taking near linear time with respect to the number of nodes n. In the algorithms below, we initialized k-means so that the starting centroids were chosen to be as close to orthogonal as possible. Initializing k-means in this way does not change the time-complexity but can significantly help to improve the quality of the clusterings, as discussed in [11] , while at the same reducing the need for multiple random restarts. In addition, both algorithms below can be run for any range of k values between a lower bound k min and an upper bound k max . When not stated otherwise, we will assume in what follows that we have k min = 1 and k max = K. In the case where k = 1, Q = 0 and the cluster is just all the vertices in the graph. 3.2 Algorithm Spectral-1 This algorithm takes as input an eigenvector matrix U K , and consists of the following steps: (a) Form the matrix U k from the first k − 1 columns of U K . (b) Scale the rows of U k using the l 2 -norm so they all have unit length (c) Cluster the row vectors of U k using k-means or any other fast vector-based clustering algorithm. For k = 1, the cluster is just the graph itself. This algorithm is similar in spirit to the one developed in [11] . Both algorithms embed the input graph into a Euclidean space by eigendecomposing a suitable matrix and then cluster the embedding using a geometric clustering algorithm. We experimentally validated the claim made in [11] that row-normalizing the matrix of eigenvectors, so that the row vectors are projected onto the unit hypersphere, gives much higher quality results. The Spectral-1 algorithm is different in three key respects to this earlier work (in addition to the minor ontological point that our framework is designed to cluster graphs while theirs is designed to cluster realvalued points): 1. Whereas in [11] the matrix that is eigendecom- 2 , was implicitly chosen to optimize the Normalized Cut, our algorithm is explicitly designed to optimize the modularity Q. 2. Our algorithm has a natural method for modelselection, the Q measure, which is the same objective function our embedding is based on. Since Normalized Cut is biased by the size of k, it can not be used for choosing the best k. 3. Our algorithm does not require an extra step of model selection to ensure: a) the edge weights are scaled correctly and b) the graph is sparsified. If links are not sparsified in the algorithm in [11] , the time complexity is O(n 3 ). Copyright © by SIAM. Unauthorized reproduction of this article is prohibited 277 3.3 Algorithm Spectral-2 The second algorithm we propose is a faster version of the algorithm in the previous section (Spectral-1). It uses a greedy strategy of recursive bisection to search over different values of k. Because of this strategy it need not find as high quality a clustering (as high a Q value) as the other approach, but it will be faster since in going from k clusters to k+1 only a portion of the data needs to be clustered rather than all of the data. The algorithm again takes as input the eigenvector matrix U K as before and consists of the following steps: 1. Initialize k, the current number of clusters, to k min . Initialize P , the best clustering assignment seen so far, to the clustering produced by running k-means with k set to k min clusters. If k min = 1, then simply initialize P to be one cluster containing all nodes in the graph. 2. Repeat until k > K or no more splits are possible: i. If not already formed, form the matrix U k from the first k − 1 columns of U K and scale the rows using the l 2 -norm so they all have unit length ii. Form the matrix U k,c from U k by keeping only rows corresponding to nodes in V c iii. Run k-means with 2 clusters on U k,c to get two new sub-clusters,V c,1 and V c,2 . iv. Form a new partition, P by setting P = P and replacing the corresponding V c with V c,1 and V c,2 v. If Q(P ) > Q(P ), accept the split by replacing the corresponding V c in P new with V c,1 and V c,2 , otherwise reject it and leave P new unchanged. vi. Assign k to be the (possibly new) number of clusters in P new The idea behind this algorithm is to start with k min clusters and instead of rerunning k-means on the entire graph for subsequent values of k as we did in the previous algorithm, we instead try recursively splitting each cluster into two child clusters if the split produces a higher value of Q. By continuing this procedure until no more splits are possible or until K clusters have been found, we end up with a clustering with the highest value of Q encountered along the way. The particular algorithm above is order-sensitive in the sense that cycling through the clusters in a different order could produce different results-however, we have not noticed any particular sensitivity to order in the data sets described later in the paper. Unlike many other recursive bisection methods, model selection here is natural and straightforward. We choose to accept a split if the split results in a higher value of Q. Of course, the drawback with this algorithm is that we tradeoff speed with accuracy. The algorithm uses a greedy strategy, where a split can never be revoked and so one bad choice negatively affects all other choices further down the same branch. Thus, the quality of the results are not necessarily as good as with Spectral-1, although they are still generally competitive as we will see below. There is a subtlety in this algorithm in that, if left unchecked, each cluster for which a split attempt was made but failed would be retried again the next time through the loop. For this reason, we only allow a single attempt to split a given cluster. O(ndke) where n is the number of data points, d is the dimensionality of each point, and e is the number of iterations required for k-means to converge. Elkan [5] proposed a much faster version of k-means. It produces exactly the same results as standard k-means but uses various geometric inequalities to significantly reduce the number of distance computations required. Elkan found that with his proposed algorithm, the overall time complexity is roughly O(nke) where e is the number of iterations. We use this algorithm for our implementation of k-means in both Spectral-1 and Spectral-2. The resulting complexity for clustering in Spectral-1, using Elkan's fast k-means algorithm, is roughly O(nK 2 e). For Spectral-2 the computational complexity is not as easy to estimate. In the worst case (completely imbalanced clusters where the largest cluster is split at each iteration) it will have the same complexity as Spectral-1. However, in practice we have found that it is considerably faster than Spectral-1 and will show experimental results later in the paper that illustrate this. In addition, for both algorithms there is the additional complexity of O(mKh + nK 2 h + K 3 h) for computing the matrix of eigenvectors U K , using the IRLM. Thus, Spectral-1 and Spectral-2 have an overall worstcase time complexity of O(mKh+nK 2 h+K 3 h+nK 2 e). Thus, for sparse graphs, where m ∼ n, the algorithms will scale roughly linearly as a function of the number of nodes n. This is in contrast to Newman's algorithm which has complexity O(n 2 ) even for sparse graphs, and thus does not scale up as well to large values of n. Experimental Results 4.1 Clustering words from WordNet We first illustrate how different choices for the graph embedding can affect the quality of clustering. We use a relatively small unweighted graph extracted from the WordNet database of word meanings and semantic relationships [8] . The reason we chose this data set was because one can immediately judge the quality of the clusters since intuitively clusters should contain words that share common semantic features which are recognizable. We created an unweighted undirected graph where nodes represent words and an edge exists between two nodes if any of the following semantic relationships exist between them: synonymy, antonymy, hypernymy, hyponomy, meronymy, troponymy, causality, and entailment. The entire graph contains 82670 nodes. We extracted a subgraph comprised of all nodes whose shortest path distance away from the word "Science" is no more than 3. We also removed all nodes with degree 1 so that the graph layout would not be too cluttered. Adding this constraint had little effect on the quality of the clusters. The resulting subgraph contains 230 nodes and 389 edges. Figure 2 shows the best clustering found by the Spectral-1 algorithm. For this clustering there were 12 clusters and Q = 0.696. Table 1 shows ten representative words from five random clusters. Figure 3 shows how the modularity Q varies with k when each of the following three types of matrices is used in step 1 of the first algorithm and the eigenvectors are computed exactly: standard Q Laplacian L Q , normalized Q Laplacian L Q , and the transition matrix D −1 W . We can see that using the transition matrix in step 1 provides a very good approximation to the normalized Q matrix. We can also see that the standard Q Laplacian slightly underperforms both of these matrices. This agrees with observations by other authors that the normalized Laplacian gives better results than the standard Laplacian (e.g., [12] , [14] ). Our next example demonstrates the ability of both of our algorithms to identify known clusters. The unweighted network was drawn from the schedule of games played between 115 NCAA Division I-A American college football teams in the year 2000. Each node in this network represents a college football team and each edge represents the fact that two teams played together. Because the true conference to which each team belongs is known a-priori, and because inter-conference games are played more often than intra-conference games, the groups of teams that form conferences correspond to natural clusters that should be identifiable in the graph. Figure 4 shows that this is indeed the case where the Spectral-1 algorithm identified the correct number of clusters and, furthermore, each team assignment to a cluster made by the algorithm was correct. 3 Our second algorithm, Spectral-2, did almost as well making very few mistakes. The mistakes were:
White and Smyth (WS) propose a fast spectral clustering algorithm in REF .
A spectral clustering approach to finding communities in graphs
In SIAM International Conference on Data Mining
Computer Science