Orthogonal and Nonnegative Graph Reconstruction for Large Scale Clustering

1 Introduction

Clustering is an important topic in machine learning and data mining. As a branch of clustering, graph clustering has been drawing growing attention which aims to group vertices of the graph into clusters with the expection that there are more edges within each cluster and fewer edges between the clusters [Schaeffer, 2007; Zhang et al., 2014; Wang et al., 2016]. As the field of graph clustering involves too many aspects, we mainly focus on spectral clustering which is considered as a way to solve relaxed versions of graph cut problems including Normalized Cut (Ncut) and RatioCut. Spectral clustering has limited applicability to process large scale data, since it uses eigenvectors of the Laplacian matrix. In general, it takes O(n3 ) for eigen-decomposition with n denoting the number of data points.

Generally, there are two major directions to solve the scalability issue. One is to reduce the computational cost of the eigen-decomposition problem. [Fowlkes et al., 2004] firstly adopted the classical Nystr¨om method to extrapolate the complete clustering results using only a small number of samples. The works in [Li et al., 2011] and [Choromanska et al., 2013] also attempted to alleviate the computational burden based on Nystr¨om method with improvements on complexities or with attached performance guarantees. Spectral sparsification can be seen in [Spielman and Teng, 2011] with the computational cost relying on the number of edges in original graph.

From the point of view of hardware platform, [Chen et al., 2011] sparsified the matrix via retaining k-nearest neighbors, paralleling both memory use and computation on distributed computers.

Another direction is to reduce the data size by sampling some representative points beforehand. [Yan et al., 2009] developed a general framework that spectral clustering only needs to run on a small subset. [Shinnou and Sasaki, 2008] adopted a committees-based way with a slightly different process to obtain the representatives. [Cai and Chen, 2014] selected a few landmarks for spectral embedding while [Peng et al., 2013] chose some in-sample data for sparse subspace clustering, and the remaining data are represented as a linear combination of the pre-selected points. [Liu et al., 2013] devoted to handling graph data, by generating supernodes to compress the original graph into a sparse bipartite graph.

Most previous works can not achieve good clustering results since they are limited by the performance of post-processing. For example, kmeans is a common way to obtain the final cluster labels, while kmeans itself is sensitive to the initialization. To get rid of this extra step, and at the same time to be able to handle the scalability issue, we propose a novel approach called orthogonal and nonnegative graph reconstruction (ONGR). It is worthwhile to highlight the main contributions of the paper as follows:

1. ONGR is proposed from the viewpoint of graph reconstruction, which is different from both two kinds of existing methods aiming at handling the scalability issue.
2. ONGR adds nonnegative constraint to the objective, thus offering interpretability that the final cluster labels can be directly obtained without post-processing.
3. Due to the orthogonal and nonnegative constraints, ONGR reconstructs the graph by a structured one which is good for clustering tasks.
4. ONGR has close relationship with spectral clustering and nonnegative matrix factorization. Comprehensive experiments on a variety of datasets show the effectiveness of the proposed method.

Notations: Suppose we have n data points belonging to k clusters with the dimensionality denoted by d. Ik ∈ Rk×k is an identity matrix and 0 is a zero matrix of proper size. The Frobenius norm is denoted by ‖ · ‖F . Let Gr = (Vs, E) be an undirected weighted graph with vertices Vs = {vi}ni=1 and edges E. The symmetric similarity matrix of the graph is W = {wij}ni,j=1 and the degree of vi is defined as di = ∑nj=1 wij. The diagonal degree matrix D has {di}ni=1 on the diagonal. The Laplacian matrix is L = D -W, and the normalized one is ˜L =D-1/2LD-1/2 .

2 Revisit of Spectral Clustering

We first review some concepts in graph cut [Schaeffer, 2007; Von Luxburg, 2007]. Given a subset A ⊂ Vs, the complement of which is denoted by ¯ A, and the size can be measured by the number of vertices denoted by |A|, or the sum of degrees of vertices denoted by vol(A) = ∑i∈A di.

Thus, we can define the cut between A and ¯ A as cut(A, ¯ A) = ∑i∈A,j∈ ¯ A wij. Further, the definitions of RatioCut [Hagen and Kahng, 1992] and Ncut [Shi and Malik, 2000] can be respectively represented as RatioCut(A1, . . . , Ak) = ∑ki=1 cut(Ai, ¯ Ai) |Ai | , NCut(A1, . . . ,Ak) = ∑ki=1 cut(Ai, ¯ Ai) vol(Ai) , where k nonempty subsets A1, . . . ,Ak satisfy Ai ⋂ Aj = ∅(i, j = 1, . . . , k), and A1 ⋃, . . . , ⋃ Ak = Vs. Minimizing these two problems is NP hard and it is natural to solve the relaxed problems.

Define k cluster indicator vectors fj = (f1j, . . . , fnj)′ by

where i = 1, . . . , n, j = 1, . . . , k. Let F ∈ Rn×k be the indicator matrix that consists of the indicator vectors as columns. It is not difficult to verify F T DF = I, f ′ i Lfi = cut(Ai, ¯ Ai) vol(Aj (i = 1, . . . , k). So the minimization problem of NCut can be rewritten as

Substituting H = D1/2F to Eq.(2), the relaxation of minimizing NCut can be obtained by allowing the entries of indicator matrix to be arbitrary real values, i.e.,

Eq.(3) is exactly the objective of normalized spectral clustering[Ng et al., 2002]. Similarly, to approximate RatioCut, the indicator vectors are defined the same as Eq.(1) with |Aj| replacing vol(Aj). We have F T F = I, and the relaxation of minimizing RatioCut can be reformulated as

which is the objective of unnormalized spectral clustering. Note that when W is a doubly-stochastic matrix (a nonnegative square matrix satisfies that row sum and column sum all equal to 1), we have ˜L = I - W = L,H = F. Then Eq.(3) becomes the same as Eq.(4), which means the equivalence between RatioCut and Ncut under the condition.

Utilizing a doubly-stochastic similarity matrix is usually good for clustering tasks [Zass and Shashua, 2006; Wang et al.,

2010; 2016]. In the subsequent sections, W is by default to be doubly-stochastic, and we try to reconstruct it by a structured graph.

3 The ProposedMethod

3.1 Formulation

In spectral clustering, the eigenvectors of Laplacian matrix can be considered as the relaxed indicator vectors, but it lacks interpretability and thus relies on the post-processing. To get rid of the extra step, we add additional nonnegative constraint to get discrete indicator vectors. The objective with two constraints is as follows:

The Ncut problem in Eq.(2) has discrete indicator vectors whose entries are also nonnegative, while spectral clustering relaxes the entries to be any real values. Intuitively, compared to spectral clustering, the objective in Eq.(5) is much closer to Ncut, thus our model tends to get better performance.

Due to the nonnegativity, the objective in Eq.(5) offers the interpretability that entries in the indicator matrix directly correspond to relationship between data points and clusters. The constraints lead F to be a matrix that there is only one non-zero entry in each row, and the ℓ2-norm of each column is 1. We conceptually illustrate the constraints in Figure 1(a). The cluster labels can then be obtained by finding the column index of the non-zero entry in each row of F . Eq.(5) is hard to tackle. Considering the computational cost, we take a circuitous way by proposing an approximated model. We first have the following proposition to transform Eq.(5).

Proposition 1. Solving Eq.(5) is equivalent to solve:

Proof. Note that W is symetric and doubly-stochastic. Tr(W T W) and Tr (FF T FF T ) are two constant terms.

Eq.(6) is important for developing the newmodel. It can be considered as a kind of graph reconstructoin with orthogonal and nonnegative constraints. The original constructed graph does not have clear structure since there usually exists noise in the data, while the graph reconstruction is an optimization process to learn a structured graph. We illustrate the reconstruction proecss in Figure 1(b). As we can see, the objective in Eq.(6) tries to reconstruct the similarity matrix by a block diagonal matrix. Such a clear structure [Nie et al., 2014] contains more accurate information about the clusters thus it is good for clustering tasks. Our new model is formulated as follows:

The first term is the reconstruction term and the second one is the regularization term which forces F and G to be close to each other. λ is the trade-off parameter that balances the two terms, and G ∈ Rn×k can be called the label matrix that gives the final cluster labels. It can be seen that Eq.(7) is equivalent to Eq.(6) for a large enough value of λ.

F is relaxed to have continuous values, but it differs from spectral clustering since F in our new model has the restriction of being close to a nonnegative matrix G. Compared to Eq.(6), Eq.(7) is easier to solve with less computational cost but still provides interpretability. In a sense, the interpretability of F is passed on to G. Entries in the label matrix G can be regarded as a soft relationship between data points and clusters, and cluster labels can be obtained by finding the column index of the largest entry in each row of G.

3.2 Optimization

To solve the non-convex problem in Eq.(7), We divide it into two subproblems that are optimized iteratively. The whole procedure is summarized in Algorithm 1 which converges fast, and in the experimental part we will see that the local optimal solution is good enough due to the effective initializaiton.

Step 1: Update G with F fixed An usual way to solve Eq.(7) with F fixed may use Lagrange multiplier method and the KKT condition. Note that the fixed F still meets the condition of F T F = I. Instead, we take a simpler way to obtain the optimal solution. By adding or removing some constant terms, we can rewrite the problem as follows:

The optimal solution of the label matrix can then be written as

where (·)+ is to make each entry of the matrix its absolute value.

Step 2: Update F with G fixed Suppose the SVD of (WG + λG) is UΛV T , and Q = V T F T U ∈ Rk×n . Eq.(7)

Algorithm 1 ONGR Algorithm

Input: X = [x1, . . . , xn] ∈ Rd×n , cluster number k, anchor number m, nearest anchor number s and trade-off parameter λ;

Output: Cluster labels Y .

• 1: Select m anchors using kmeans or random selection;
• 2: Construct a sparse regression matrix Z by Eq.(12);
• 3: Initialize F ∈ Rn×k by left singular vectors of ˆ Z;
• 4: repeat
• 5: Update G by Eq.(9);
• 6: Update F by Eq.(11);
• 7: until converges.
• 8: Find the column index of the largest entry in each row of the label matrix G.

with G fixed can be rewritten as:

Λii, Qii are the (i, i)-th entry of Λ and Q, respectively. Note that Λii is singular value and QQT = Ik, thus we have Λii ≥ 0 and Qii ≤ 1. Therefore, ∑i ΛiiQii ≤ ∑i Λii and the equality holds when {Qii = 1}ki=1 , i.e., Q = [Ik, 0]. Recall that Q = V T F T U, so the optimal solution to Eq.(10) is

Graph Construction The traditional KNN-based way to construct graph is time consuming. A normalized step is also needed to insure Proposition 1. Here we adopt an efficient way to compute a doubly-stochastic similarity matrix.

Given data matrix X = [x1, . . . , xn] ∈ Rd×n , and the subset Xsub = [u1, . . . , um] ∈ Rd×m in which each point plays a role as an anchor. According to [Liu et al., 2010], we can design a sparse regression matrix Z ∈ Rn×m as follows:

where 〈i〉 is the set containing the indexes of s nearest anchors of xi. We choose Kh(·) as the common used Gaussian kernel with self-tuning bandwidth h [Chen et al., 2011]. The doubly-stochastic similarity matrix is then computed by W = ZΣ-1Z T , where Σ is a diagonal matrix with the entry Σkk = ∑ni=1 Zik.

Important tips: We do not need to calculate W explicitly since it would be both time and memory consuming. We can utilize the low-rank matrix multiplication to speed up. Concretely, we can compute (WF ) and (WG) by Z (Σ-1 (Z T F )) and Z (Σ-1 (Z T G)), respectively.

Initialization The indicator matrix F can be randomly initialized with orthogonal constraint. An alternative way is to use the k eigenvectors corresponding to the k smallest eigenvalues of ˜L, or the k eigenvectors corresponding to the k largest eigenvalues of W.

Let ˆ Z = ZΣ-1/2 , we have W = ˆ Z ˆ Z T . Supposing the SVD of ˆ Z is ˆ Z = U1Λ1V T 1 , it's easy to find that the left singular vectors of ˆ Z are the eigenvectors of W. Therefore, we can initialize F by taking the largest k singular vectors in U1 as columns.

3.3 Computational Complexity

Suppose the number of anchor points is m. It takes O(mndt1) to find the anchors and takes O(m2n +m3 ) to initialize F by using left singular vectors of ˆ Z, where t1 denotes the number of iterations in kmeans. The construction of Z needs O(mnd) while the alternating and iterative procedure needs O ((mnk +m2k + nk2 + k3 )t2 ). t2 denotes the number of iterations, and O(nk2 + k3 ) is the cost of SVD of M. For large scale problem, there usually exist m > k and n ≫ k, m. Thus the total complexity is about O(mndt1+m2n+mnkt2) that scales linearly with the data size n. See Table 1 for a brief summary of time complexity of different methods (references can be seen in the section of comparison algorithms).

Except for KNN-SC that costs too much on constructing the similarity matrix, all other methods are linear time complexity while they all depend on the post-processing step. SSSC is more sensitive to the dimensionality d, which makes it unsuitable to high dimensional data such as images. The total complexities of Nystr ¨om and LSC are O(mnd+m2n+nk2 t0), O(mndt1 +m2n+nk2 t0), respectively.

3.4 Discussions

Our proposed model has similar formulation with non-negative matrix factorization (NMF) [Lee and Seung, 2001; Li and Ding, 2006; Ding et al., 2005]. In this subsection, we compare ONGR with NMF methods and show the essential differences between them.

Given a nonnegative matrix X ∈ Rd×n , and a reduced rank k, the problem of NMF is formulated as follows:

where ˜ F ∈ Rd×k , ˜G ∈ Rn×k . Due to its nonnegativity, NMF provides interpretability that the cluster labels can be obtained by finding the index of the largest entry in each row of ˜G.

Compared to our models in Eq.(6) and Eq.(7) that are based on graph reconstruction, NMF in Eq.(13) is based on data reconstruction and requires the input data to be nonnegative. In many cases, the data points lies on a nonlinear manifold. It is more suitable to reconstruct the graph which describes the relationship between data points, rather than reconstructing the data matrix directly.

Recently, a nonnegative symmetric factorization (SymNMF) [Kuang et al., 2012] of the similarity matrix is proposed as follows:

where ˜H ∈ Rn×k , and A ∈ Rn×n can be any symmetric matrix representing similarity values. Again, due to the nonnegativity, the column index of the largest entry in each row of ˜H indicates the cluster label.

Comparing our model in Eq.(6), the normalized spectral clustering in Eq.(3) and SymNMF in Eq.(14), it is easy to find they share the common objective, while the constraints are totally different. Specifically, our model in Eq.(6) is a combination of the other two by introducing the orthogonal and nonnegative constraints simultaneously, thus the reconstructed graph has more clear structure. Our new model in Eq.(7) naturally inherits the good property for clustering, while it has low complexity and can be applied to large scale clustering.

There are many other NMF methods, such as Semi-NMF, Convex-NMF and Tri-factorization. They have different constraints or assumptions. These NMF methods are based on matrix multiplication updating rules and they often take a long time to converge. To seek a trade-off between convergence rate and computational cost, SymNMF developes a Newtonlike algorithm that still has complexity of O(n3k) in each iteration.

4 Experiments

4.1 Datasets and EvaluationMetrics

Datasets A variety of datasets are adopted to evaluate the proposed method. They can be downloaded from the UCI Machine Learning Repository 1 , the LibSVM Data page 2 , and three webpages 3 4 5 . Most of these datasets have more than 10,000 samples.

We obtained the MINIST-extend dataset by translating the original images in MINIST by one pixel in each direction like [Liu et al., 2010]. We scaled CoverType to [0, 1] by feature. The letter symbols in Connect-4 were replaced by digits. For USPS, COIL-20, COIL-100, MINIST and its extended version, we firstly scaled them to [0, 1], then we conducted dimensionality reduction by PCA with 95% information reserved, except for MINIST-extend which reserves 90% information. Details can be seen in Table 2.

1 https://archive.ics.uci.edu/ml/datasets.html

2 https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/

3 http://alumni.cs.ucsb.edu/ wychen/sc.html

4 http://yann.lecun.com/exdb/mnist/index.html

5 http://www.cs.columbia.edu/CAVE/software/

EvaluationMetrics The clustering quality is measured by Normalized Mutual Information (NMI) and Accuracy (Acc) [Cai et al., 2005]. The values of NMI and Acc range from 0 to 1 with higher score corresponding to better performance. Note that the clustering accuracy is the average performance of label matching results between ground truth labels and predicted labels, which is different from the classification accuracy.

4.2 Comparison Algorithms

We provide two versions of ONGR in which ONGR-K uses kmeans to find anchors and ONGR-R adopts random selection. The comparison algorithms are (1) Nystr¨om [Fowlkes et al., 2004]: Spectral clustering that uses Nystr¨om method with orthogonalization. (2) KNN-SC [Chen et al., 2011]: Spectral clustering that uses KNN to construct similarity matrix. (3) LSC [Cai and Chen, 2014]: LSC-K uses kmeans for landmark selection, and LSC-R randomly selects landmarks. (4) SSSC [Peng et al., 2013]: Scalable sparse subspace clustering. Due to the requirement for nonnegativity of data matrix or high computational cost, NMF methods are not considered as the comparison methods.

4.3 Experimental Setting

There are three parameters in ONGR, namely m, s and λ. LSC needs to tune parameters p, r. Nystr¨om has parameter n1, and KNN-SC has parameter k. The three parameters n2, λ1 and δ are from SSSC.

For fair comparison, we took the same kmeans centroids as anchors or landmarks in ONGR-K and LSC-K, and also took the same random selection points in ONGR-R, LSCR, Nystr¨om and SSSC. For m, p, n1, n2, we searched their value in the range of [100,1200] with step size 100, while searching s and r in the range of [2,8] with step size 1. For KNN-SC, we chose k in the range of [5,20] with step size 5. As suggested by the authors of SSSC, we searched λ1 and δ among {10-7 , 10-6 , 10-5 } and {10-3 , 10-2 , 10-1 }, respectively. For our trade-off parameter λ, we searched log λ in the range of [-6,3] with step size 1.

We ran each algorithm except for ONGR for 20 times under each parameter setting. Note that given the pre-selected points, ONGR is the only method that has stable performance. For the last three datasets in Table 2, we did not test ONGR-K and LSC-K since kmeans becomes costly for finding anchors or landmarks. We also did not test KNN-SC due to the unbearable running time.

Convergence judgement: With data growing in scale, it is a heavy burden for us to compute the objective value of Eq.(7). Recall that we do not calculate W explicitly. We adopted another criterion, in which the (i+1)-th iteration is terminated when just a very small percentage of data points change their predicted labels comparing to the i-th iteration. Empirically, we set the threshold to be 0.001.

4.4 Experimental Results

Table 3 records the running time corresponding to the best performance on five datasets. The results are consistent with the complexity analysis. Table 4 reports the clustering performance along with the standard deviation (std), and we can see that ONGR achieves the best or second best results no matter what the metric is. The average results of each algorithm on all datasets is also reported, showing the superiority of ONGR more clearly.

From Table 3 and Table 4, it can be concluded that ONGR achieves stable and much better performance in the shortest or relatively less time. Specifically, within about 20 and 200 seconds, ONGR-R gains 8.14% and 9.24% increment of accuracy over the second best results on PenDigits and CoverType, respectively. Comparing LSC and ONGR, the two methods with the common graph construction, we see that ONGR-R exceeds 5.98% than LSC-R and ONGR-K exceeds 7.41% than LSC-K with respect to accuracy, which shows the advantage of removing post-processing.

To verify the convergence of the proposed method, in the upper part of Figure 2(a), we plotted the objective value of Eq.(7) on the relatively small dataset COIL-20. The curve of the squared Frobenius norm of (F -G) demonstrates that F

and G do get close to each other in each iteration. Thus the model in Eq.(7) gets closer and closer to the one in Eq.(6), and the inheritance of good clustering quality from Eq.(6) to Eq.(7) is guaranteed. As can be seen, ONGR does converge fast within about 30 iterations.

In Figure 2(b), we plotted the curves of clustering performance on MINIST-extend. It can be seen that the performance is getting better as the iteration number increases, showing the effectiveness of the optimization process. Note that in order to completely present the curves, we set the maximum iteration number to be 80 and did not terminate the iteration.

As λ is the only parameter introduced by the model in Eq.(7), in Figure 2(c), we examined the effect of λ to clustering performance on USPS. Parameters m and s were set to be 800 and 4, respectively. Under the condition, we see that ONGR is pretty robust to λ.

Moreover, in Figure 3, we tested ONGR-R on COIL-20 to give a practical illustration of graph reconstruction. As we can see, compared to the constructed graph, the reconstructed graph has clear structure which contains more accurate information about the clusters. However, the zero rows in the label matrix lead to the missing block in the reconstructed graph, i.e., zeros rows will destroy the structure. This is a topic that we will explore to further improve the clustering performance.

5 Conclusion

In this paper, we have proposed a novel approach called ONGR for large scale clustering, which is based on the viewpoint of graph reconstruction. With orthogonal and nonnegative constraints, the reconstructed graph naturally has clear structure about the clusters. Due to the nonnegativity, the interpretability is provided and the post-processing is no longer needed. Given the anchors, ONGR has stable and much better performance than other state-of-the-art methods, demonstrated by extensive experiments.

In the future, we plan to study the sparsity of label matrix to avoid zero rows. We are also going to propose a robust version to better deal with outliers. The ℓ1-norm and ℓ2,1-norm may be considered to achieve robustness and sparsity.

Acknowledgments

This work was supported in part by the National Science Foundation of China under Grants 61522207 and 61473231.