Scalable Normalized Cut with Improved Spectral Rotation

1 Introduction

Clustering is a hot topic in machine learning and data mining. Over the past decades, many clustering algorithms have been proposed for cluster analysis of high-dimensional data, such as spectral clustering [Von Luxburg, 2007], subspace clustering [Kriegel et al., 2009; Chen et al., 2012], multi-view clustering [Cai et al., 2011; Chen et al., 2013], etc. Among them, spectral clustering is a popular method because it is easy to implement and often shows good clustering performance due to the use of manifold information. Various spectral clustering algorithms have been proposed, such as Ratio Cut [Hagen and Kahng, 1992], k-way Ratio Cut [Samaria and Harter, 1995], Normalized Cut [Ng et al., 2002], Spectral Embedded Clustering [Nie et al., 2011] and MinMax Cut [Nie et al., 2010]. They have been successfully applied to many highdimensional applications, such as image segmentation [Shi and Malik, 2000; Yu and Shi, 2003], clustering gene expression data [de Souto et al., 2008] and power network decomposing [Sánchez-García et al., 2014].

*Corresponding Author.

Spectral clustering methods usually transform the data into a weighted, undirected graph based on pairwise similarities. To obtain the final discrete clustering assignments, they often perform eigendecomposition first, and then the final clustering assignments can be obtained from eigenvectors by kmeans or spectral rotation [Yu and Shi, 2003]. According to the analysis in [Huang et al., 2013], spectral rotation can obtain better clustering result than k-means. However, spectral rotation involves a two-stage process in which an approximate continuous cluster assignment matrix is first computed, and the final discrete solution is a nearby discrete solution obtained from the approximate continuous cluster assignment matrix. A disadvantage of this two-stage process is that the final clustering structures may deviate from the true discrete solution.

Moreover, since both graph construction as well as spectral analysis are time consuming, spectral clustering usually has a time complexity of O(n3 ) where n is the number of samples. In recent years, much effort has been devoted for accelerating the spectral clustering. There are mainly two ways to handle the scalability issue of spectral clustering. One way is to reduce the computational cost of the eigendecomposition step [Fowlkes et al., 2010; Li et al., 2010], and another way is to sample the original data and perform clustering on the reduced data [Yan et al., 2009; Shinnou and Sasaki, 2008]. However, these methods are based on sampling, and a lot of information of the data will be lost in the sampling step. Recently, Cai et al. proposed a landmarks-based spectral clustering (LSC) method [Cai and Chen, 2015]. Given a data set with n samples, LSC generates m ≪ n representative data points to compute a representation matrix and the eigendecomposition can be performed on the low-size representation matrix. The final discrete clustering result is obtained from eigenvectors by k-means. The overall time of LSC is O(ndmt + nm2 ) where t is the number of iterations of k-means for anchor generation, which is significant reduction from O(n3 ) considering m ≪ n. However, how to effectively construct a representation matrix and how to effectively obtain the clustering assignments are still two problems that need to be solved.

In this paper, we propose a Scalable Normalized Cut method (SNC) for large scale data. Given a data set with n samples, we first use k-means to find m ≪ n representative data points, a new method to construct a low-size n ×m representation matrix on which the eigendecomposition can be performed. We propose an Improved Spectral Rotation method to obtain the final clustering assignments. SNC has the same computational complexity as LSC for large scale data. The main contributions of our work include:

1. We propose an Improved Spectral Rotation (ISR) to obtain the solution of the final clustering assignments.
2. We propose an efficient method to construct a small representation matrix, which can be used to computed an affinity matrix. We futher prove that the resulting affinity matrix is symmetric and doubly stochastic.
3. Comprehensive experiments on 14 benchmark data sets show the efficiency and effectiveness of the proposed method.

The rest of this paper is organized as follows. Notations and preliminaries are given in Section 2. We review the background and related work in Section 3. The Improved Spectral Rotation (ISR) is given in Section 4 and the Scalable Normalized Cut (SNC) is given in Section 5. We present experimental results and analysis in Section 6. Conclusions and future work are given in Section 7.

2 Notations and Definitions

We summarize the notations and the definition of norms used in this paper. Matrices are written as boldface uppercase letters. Vectors are written as boldface lowercase letters. For matrixM= (mij), its i-th row is denoted as mi , and its j-th column is denoted by mj . The Frobenius norm of the matrix

3 Background and RelatedWork

In this section, we introduce the anchor-based similarity matrix construction and spectral rotation.

3.1 Anchor-based SimilarityMatrix Construction

To handle the scalability issue of spectral clustering, Liu et al. proposed an anchor-based strategy [Liu et al., 2010], which is also called landmarks-based method [Cai and Chen, 2015]. Given a data set X ∈ Rd×n with n objects {x1, x2, ..., xn}, anchor-based strategy first seeks m anchors, where m ≪ n, and then construct the affinity matrix by calculating the distance between anchors and original samples. There are mainly two methods for anchor generation, i.e., random selection and k-means generation. Since clustering centers have a stronger representation power than random selected data, it is preferred to use k-means for anchor generation [Liu et al., 2010; Cai and Chen, 2015].

After we have m anchors W ∈ Rd×m, the next step is to obtain a representation matrix B such that X ≈ WB. With B, we can obtain an affinity matrix A as [Liu et al., 2010]

where ∆ ∈ Rm×m is a diagonal matrix and the j-th entry is defined as ∆jj = ∑ni=1 bij. The most important property of this similarity matrix is that it can be represented as A = PPT where P ∈ Rn×m = B∆- 1 2 .

3.2 Spectral Rotation

In this subsection, we introduce the spectral rotation which is used in Multiclass Spectral Clustering (MSC) [Yu and Shi, 2003]. Given an affinity matrix A, we can compute the corresponding degree matrix DA, which is a diagonal matrix with the i-th diagonal element as dii = ∑nj=1 aij. The objective function of MSC is

where Z ∈ Rn×c is the scaled partition matrix. It is hard to directly solve problem (2). A well known way is to relax the matrix Z from the discrete values to the continuous ones, and form the new problemmax

According to Proposition 1 in [Yu and Shi, 2003], the optimal solution of Z is {Z∗R : RT R = Ic} where Z∗ ∈ Rn×c is the c column vectors of the eigenvectors of D-1 A A which correspond to the c biggest eigenvalues.

To obtain the discrete solution Y, we first compute an approximate Y∗ asY∗ = Diag(Z∗ (Z∗ )T )- 12 Z∗ (4)

Then we can learn suitable R and Y such that Y∗R is closest to Y by solving the following problem

4 Improved Spectral Rotation

In MSC, approximate Y ∗ is first computed and then a suitable R is learnt for the final cluster indicator matrix Y. However, the final clustering results may deviate from the true discrete solution since Y ∗ is an approximate solution. In this paper, we propose a new spectral rotation method to obtain better discrete solution of Y. We first rewrite problem (2) as follows

where F ∈ Rn×c is the cluster indicator matrix. We can relax F to continuous matrix and form the new problem

It can be verified that the optimal solution of F is {F∗R : RT R = Ic} where F∗ ∈ Rn×c is the c column vectors of the eigenvectors of D- 1 2 A AD- 12 A which correspond to the c biggest eigenvalues.

With F∗ , the next step is to obtain the discrete solution of Y. In this paper, we propose to directly obtain the discrete solution Y by solving the following problem

Note that ∥ ∥∥D 12AY(YT DAY)- 12 -F∗R∥∥∥2 F = 2n - 2T r((YT DAY)- 1 2 YT D 12AF∗R), problem (8) can be rewritten as

We can apply the alternative optimization approach to solve problem (9).

4.1 Update Rwith Y fixed

If Y is fixed, denote (YT DAY)- 12YT D 12A as M ∈ Rc×n . Suppose the SVD of MF∗ is MF∗ = UΣVT , then we have

where E = VT RU, λii and eii are the (i, i) - th element of matrix Σ and E respectively.

Since ET E = UT RVVT RT U = Ic, i.e, ∑cj=1 e2ji = 1, we know eii ≤ 1(1 ≤ i ≤ c). On the other hand, λii ≥ 0 since λii is singular value. Therefore, Tr(MF∗R) = ∑ci=1 λiieii ≤ ∑ci=1 λii, and the equality holds when eii = 1 (1 ≤ i ≤ c). That is to say, Tr(MF∗R) reaches its maximum when E = Ic. Then we obtain the optimal solution of R as

4.2 Update Y with R fixed

Let G = F∗R. According to problem (9), the optimal solution of Y can be obtained by solving the following problem

which can be rewritten as

Since √yT j DAyj involves all rows of Y, we propose to sequentially solve Y row by row and fix the other rows of Y

as constants. Suppose we have obtained the optimal solution Y, which has the objective function J old (Y). To solve the i-th row yi ∈ Bc , we only need to consider the increment of the objective function value from yij = 0 to yij = 1. Since yTj DAyj and ∑nt=1 √dttytjgtj can be computed once before we solve yi , we can compute the increment as follows (See Figure 1)

Then the optimal solution of yi can be obtained as

where < . > is 1 if the argument is true or 0 otherwise, and sij is defined in Eq. (14).

4.3 Initialization of Y

We can use the the mapping in [Yu and Shi, 2003] to obtain the initial Y. We first compute an approximate Y∗ as

Then the initial discrete solution of Y is given by

4.4 The Optimization Algorithm

The detailed algorithm to solve problem (9), named Improved Spectral Rotation (ISR), is summarized in Algorithm 1. In the new algorithm, we need O(r1(c3 + r2nc)) time to iteratively solve R and Y where r1 is the number of iterations to update R and r2 is the average number of iterations to update Y. Considering that c ≪ n for large scale data, the computational complexity for obtaining Y is O(r1r2nc). If we use k-means to obtain Y, we need O(tnc2 ) time where t is the number of iterations. Here, the discrete solution Y converges very fast due to its limited solution space so r2 is usually very small. Therefore, ISR has almost similar computational complexity as k-means for large scale data.

Algorithm 1 Improved Spectral Rotation (ISR) to solve problem (9)

• 1: Input: F∗ .
• 2: Initialize Y according to Eq. (17).
• 3: repeat
• 4: Update R according to Eq. (11), and G= F∗R.
• 5: repeat
• 6: Update Y according to Eq. (15).
• 7: until Y does not change
• 8: until problem (9) converges
• 9: Output: the clustering result Y.

5 The Scalable Normalized Cut for Large Data

In this section, we propose a Scalable Normalized Cut (SNC) for large scale data.

5.1 An Efficient Method for Construction Representation Matrix

Assume that we have obtained m anchors W ∈ Rd×m with k-means, the next step is to construct a representation matrix B ∈ Rn×m. Inspired from the work in [Nie et al., 2016], we assume that bij should be larger if xi is closer to wj and propose an efficient method to construct B. For the i-th sample xi ∈ X, we propose to obtain bi ∈ Rm by solving the following problem

According to the analysis in [Nie et al., 2016], the optimal solution bi to problem (18) is

where di,h is the square of Euclidean distance between xi and its h-th nearest neighbor, and Nk(xi) contains the k nearest neighbors of xi.

After obtaining the representation matrix B, we can compute the affinity matrix A according to Eq. (1). The following theorem ensures that A is symmetric and doubly stochastic.

Theorem 1. Given the representation matrix B computed according to Eq. (19), A computed from Eq. (1) is symmetric and doubly stochastic.

Proof. According to Eq. (1), we have

It can be easily verified that aij = aji, which indicates that A computed from Eq. (1) is symmetric.

We can also verify that

which implies that ∑ni=1 aij = ∑nj=1 aij = 1. Therefore, A is doubly stochastic.

5.2 The OptimizationModel

According to Theorem 1, it can be verified that the degree matrix of A should be an identity matrix. Then problem (6) can be rewritten as

We also relax F to continuous matrix, and obtain the optimal solution of F from the following problem

Note that A can be rewritten as A = PPT where P ∈ Rn×m = B∆- 12 , we can perform SVD on P instead of A. Suppose the SVD of P is P = UPΣPVTP , we have A = PPT = UPΣ2P UTP , which can be rewritten as

then we know that the optimal solution of F to problem (23) is the c column vectors in UP corresponding to c biggest diagonal entries in diagonal matrix Σ2P .

With the learnt optimal solution of F∗ , we can use Algorithm 1 to obtain the final solution of Y. Since the size of P is n × m, we can obtain F∗ within O(nm2 ).

Following the same analysis in [Nie et al., 2011], it can be verified that although 1n is a trivial vector in UP , it should be retained in order to generate the whole set of optima.

5.3 The Optimization Algorithm

The detailed algorithm to solve problem (22), named Scalable Normalized Cut (SNC), is summarized in Algorithm 2. Given a data matrix X ∈ Rd×n , we need O(ndmt) time to obtain m anchors by k-means where t is the number of iterations, O(ndm + nmlog(m)) time to construct P, O(nm2 ) time to obtain F∗ , and O(r1(c3 + r2nc)) time to iteratively solve R and Y where r1 is the number of iterations to update R and r2 is the average number of iterations to update Y. Here, the discrete solution Y converges very fast due to its limited solution space so r2 is usually very small. Considering that m ≪ n, d ≪ n and c < m for large scale data, the overall computational complexity is O(ndmt + nm2 ). Therefore, SNC has the same computational complexity as LSC.

Algorithm 2 Scalable Normalized Cut (SNC) to solve problem (22)

• 1: Input: DatamatrixX ∈ Rd×n , number of nearest neighbors k, number of anchors m.
• 2: Find m anchors W using k-means, and construct a sparse representation matrix B ∈ Rn×m with the method introduced in Section 5.1.
• 3: Obtain P ∈ Rn×m = B∆- 12 , where ∆ ∈ Rm×m is the degree matrix of B.
• 4: Perform SVD on P such that P = UPΣPVTP , then form F∗ by selecting c column vectors in UP which corresponds to c biggest diagonal entries in diagonal matrix Σ2P .
• 5: Call Algorithm 1 with input F∗ to obtain the optimal solution of Y.
• 6: Output: the clustering result Y.

6 Experimental results and analysis

In this section, we present the experiments conducted on 14 real-life data sets to demonstrate the efficiency and effectiveness of the proposed method.

6.1 Experiments on ISR

We first compared ISR with k-means and the original spectral rotation for spectral clustering.

Benchmark data sets

8 benchmark data sets were selected from the UCI Machine Learning Repository and Feiping Nie's page 1 . Table 1 summarizes the characteristics of these 8 data sets.

Results and Analysis

We compared ISR with k-means (KM) and the original spectral rotation (SR) for normalized cut. For each data set, we set five neighborhood parameters k = {10, 20, . . . , 50} to construct five affinity matrices with the method in [Nie et al., 2016], and used these matrices to run three methods in order to perform fair comparison. For each algorithm on each data set, we computed the average accuracy and NMI and show them in Table 2. From these figures, we can see that ISR outperformed other methods in both accuracy and NMI on almost all data sets. Especially on D1, D4 and D8, ISR has over 5% improvement compared to the second best method. This indicates that ISR improves the original spectral rotation.

We selected D8 to show the convergence curves of the objective function value and the number of iterations for obtaining Y in each main loop. The results are drawn in Figure 2. From Figure 2(a), we can see that the objective function value drops very fast, indicating Algorithm 1 converges very fast. From Figure 2(b), we can see that the average number of iterations for obtaining Y is around 4. Therefore, ISR can quickly obtain the final clustering assignments.

6.2 Experiments on SNC

In this subsection, we compare SNC with the original normalized cut and other scalable spectral clustering methods.

1 http://www.escience.cn/people/fpnie/ index.html#

Figure 2: Objective function values of problem (9) and no. of iterations for obtaining Y in each main loop of ISR on D8.

Benchmark data sets

6 large scale benchmark data sets were selected from the UCI Machine Learning Repository and Feiping Nie's page 1 . Table 3 summarizes the characteristics of these 6 data sets.

Results and Analysis

We compared SNC with six spectral clustering methods, including NCut with k-means (NCut+KM) [Ng et al., 2002], multiclass spectral clustering (MSC) [Yu and Shi, 2003], kmeans-based approximate spectral clustering (KASP) [Yan et al., 2009], committees-based spectral clustering (CSC) [Shinnou and Sasaki, 2008], parallel spectral clustering (PSC) [Chen et al., 2011] and LSC [Cai and Chen, 2015]. For each data set, we used the same clustering result for anchors generation in KASP, CSC, LSC and SNC where 10 numbers were selected for m. The neighborhood parameters were set as {10, 20, . . . , 50} for all data sets. We used the Gaussian kernel to compute similarities for all methods excluding SNC, where the parameter h was set as the average distance between two points in the data set (used in [Cai and Chen, 2015]). The average clustering performance of seven spectral clustering algorithms are shown in Figure 3. From these figures, we can see that SNC outperformed other methods in accuracy and NMI on almost all data sets. Especially on D1, SNC has nearly 10% improvement compared to the second best method NCut+KM in terms of both accuracy and NMI. On five data sets, SNC outperformed both NCut+KM and MNC which perform clustering with the similarity matrix computed from the original data. From Figure 3(c), we can see that the time costs of SNC are much smaller than NCut+KM and MNC, especially on D2, D3 and D4. The time costs of SNC are similar as those of LSC. Although SNC spent more time than KASP, CSC and PSC, it produced better results than these methods.

Parameter study

(c) Time costs of similarity con-

struction versus m and k.

(d) Time costs of clustering ver-

sus m and k.

Figure 4: Accuracy, nmi and running time of SNC versus the no. of anchors m and neighborhood parameter k.

We select D6 to show the relationship between the clustering performance and two parameters m, k in SNC. The results are drawn in Figure 4. From these figures, it can be seen that SNC can achieve better clustering results (in terms of both accuracy and nmi) as both m and k increase. From Figure 4(c), we can see that the time cost of similarity construction grows linearly as m increases, and the time cost of clustering does not change too much as k increases. From Figure 4(d), we can see that the time cost of clustering is insensitive to both m and k. Since the time cost of similarity construction is much smaller than the time cost of clustering, we can say that the total time cost of SNC is nearly insensitive to both m and k.

7 Conclusions

In this paper, we have proposed a Scalable Normalized Cut (SNC) method for large scale data, in which a parameter-free method is proposed to construct a representation matrix, and an Improved Spectral Rotation (ISR) method is proposed to obtain the final clustering assignments. Experimental results show that ISR can obtain better clustering results than kmeans and the original spectral rotation. Comparison results with other scalable spectral clustering methods show that our method can obtain better results without increasing running time too much. Therefore, the new method is effective and efficient for large scale data. In the future work, we will study new anchor generation method.

Acknowledgement

This research was supported by NSFC under Grant no.61305059, 61473194 and U1636202.