Fast Nonnegative Matrix Tri-Factorization for Large-Scale Data Co-Clustering

1 Introduction

Clustering, which partitions a data set into different groups unsupervisedly, is one of the most fundamental topic in statistical learning. Most traditional clustering algorithms are designed for one-side clustering [Nie et al., 2009], i.e. cluster either data points or features. However, in many real world applications, the clustering based analysis is interested in twoside clustering results, i.e. group the data points and features simultaneously, e.g., "documents" and "words" in document analysis, "users" and "items" in collaborative filtering, "samples" and "genes" in microarray data analysis, etc. Typically, instead of being independent, the different clustering tasks on data and features are closely correlated, and it is challenging for traditional clustering algorithms to utilize the data and features interdependence efficiently. Consequently, coclustering techniques, which aim to cluster both data and features simultaneously by leveraging the interrelations between them, have been proposed in recent researches. To name a few, Dhillon [Dhillon, 2001] introduced a bipartite spectral graph partition approach to co-cluster words and documents; Cho et al. [Cho et al., 2004] suggested to co-cluster the experimental conditions and genes for microarray data by minimizing the sum-squared-residue; Long et al. [Long et al., 2006] presented a relation summary network model to co-cluster the heterogeneous data on a k-partite graph, and so on.

More recently, Ding et al. [Ding et al., 2005] explored the relationships between Nonnegative Matrix Factorization (NMF) [Lee and Seung, 1999; 2001] and K-means/spectral clustering, and proposed to use Nonnegative Matrix Trifactorization (NMTF) [Ding et al., 2006] to co-cluster words and documents at the same time. Due to its mathematical elegance and encouraging empirical results, NMTF method has been further developed to address various aspects of coclustering [Wang et al., 2008; Li et al., 2010; Gu and Zhou, 2009; Ding et al., 2010]. However, a notorious bottleneck of NMTF based co-clustering approaches is the slow computational speed because of intensive matrix multiplications involved in each iteration step of the solution algorithms, which makes these approaches hard to be applied to largescale data in real world applications. In this paper, we propose a novel Fast Nonnegative Matrix Tri-factorization (FNMTF) approach to efficiently conduct co-clustering on largescale data. Our new algorithms are interesting from the following perspectives:

• Instead of enforcing traditional nonnegative constraints on the factor matrices of NMTF, we constrain them to be cluster indicator matrices, a special type of nonnegative matrices. As a result, the clustering results of our approach are readily stored in the resulted factor matrices. However, existing NMF based methods require an extra post-processing step to extract cluster structures from the factor matrices, which often leads to nonunique clustering results.
• Due to the nature of indicator matrices, the optimization problems of our approach can be decoupled into subproblems with much smaller sizes, and the decoupled subproblems involve much less matrix multiplications. Therefore, our approach is computationally efficient and

scale well to large-scale input data.

• Taking into account the manifold structures in both data and feature spaces, we further develop a Locality Preserved FNMTF (LP-FNMTF) approach to incorporate manifold regularizations. Efficient algorithm to optimize the objective with quick convergence is presented.
• Promising experimental results on five benchmark data sets show that our approaches not only are faster than state-of-the-art co-clustering methods but also have competitive clustering performance.

Notations and problem formalization. Throughout this paper, we write matrices as boldface uppercase letters and vectors as boldface lowercase letters. Given a matrix M = (mij), its i-th row and j-th column are denoted as mi· and m·j respectively.

Traditional clustering methods focus on one-side clustering, i.e., clustering the data side based on the similarities along the feature side. In the co-clustering problem, we cluster data points based on the distributions of features, meanwhile cluster features based on the distributions of the data points. Formally, given a data set X = {x·i ∈ Rd}ni=1 , we write X = [x·1, . . . , x·n] = [xT1·, . . . , xTd· ]T . Our goal is to group the data points {x·1, . . . , x·n} into c clusters {Cj}cj=1 , and simultaneously group the features {x1·, . . . , xd·} into m clusters {Wj}mj=1 .

We use a partition matrix G = [gT 1· , . . . , gTn· ]T ∈ {0, 1}n×c to represent clustering result of data points, such that gij = 1 if x·i belongs to cluster Cj and gij = 0 otherwise. Similarly, we use another partition matrix F = [f T 1· , . . . , f T d· ]T ∈ {0, 1}d×m to represent the clustering results of features. Here, we call F and G as cluster indicator matrices, because each row of them, i.e., fi·(1 ≤ i ≤ d) or gi·(1 ≤ i ≤ n), has one and only one element equal to 1 to indicate the cluster membership, while the rest elements are 0. We denote the set of all cluster indicator matrices as Ψ.

2 Fast nonnegative matrix tri-factorization (FNMTF) for co-clustering

In this section, we first briefly review the background of coclustering using NMTF, which motivates the optimization objective of our approach. After that, an efficient algorithm to solve our objective will be introduced.

2.1 Objective of FNMTF approach

K-means clustering is a standard clustering method in statistical learning, which minimizes the following objective:

where ‖·‖ denotes the Frobenius norm of a matrix, c·j is the j-th centroid of the data set. Because G ∈ Ψ is a cluster indicator matrix, minimizing J1 is a combinatorial optimization problem, which is hard to be resolved in general. Therefore, the minimization of J1 is often relaxed to maximize the following objective [Zha et al., 2001; Ding and He, 2004]:

Note that, G in J2 is no longer an indicator matrix, but an arbitrary orthonormal matrix.

Recently, Ding et al. [Ding et al., 2005] proved the equivalence between the relaxed objective of K-means clustering in Eq. (2) and the NMF objective when orthonormal constraints are enforced on the factor matrices, which minimizes:

where X ∈ Rd×n + , F ∈ Rd×c + and G ∈ Rn×c + , and J3 aims to approximate the nonnegative data matrix X by the product of F and G. The orthonormal constraints here ensure the uniqueness (up to a permutation) of the solution, and together with the nonnegative constraints make the resulted F and G approximate the K-means clustering results on both features and data points (called as "soft labels") [Ding et al., 2005; 2006]. The latter, simultaneously clustering the rows (features) and the columns (data points) of an input data matrix, is one of the main strength of NMF defined in Eq. (3).

Because the two-factor NMF in Eq. (3) is restrictive, which often gives a rather poor low-rank matrix approximation, one more factor S ∈ Rm×c + was introduced to absorb the different scales of X, F and G. This leads to NMTF [Ding et al., 2006]: X ≈ FSGT , where F ∈ Rd×m + and G ∈ Rn×c + . S provides increased degrees of freedom such that the lowrank matrix representation remains accurate, while F gives row clusters and G gives column clusters. In order to achieve additional flexibility, in clustering scenarios, the nonnegative constraint on X (thereby the nonnegative constraint on S) can be relaxed [Ding et al., 2010], which leads to the semi-NMTF problemminimizing the following objective:

Despite its mathematical elegance, Eq. (4) suffers from two problems that impede its practical use. First, similar to Eq. (2), the relaxations on F and Gmake the immediate outputs of Eq. (4) are not cluster labels, which require an additional post-processing step and often lead to non-unique solutions. Second, and more important, Eq. (4) is usually solved by alternately iterative algorithms, and in each iteration step the intensive matrix multiplications are involved [Ding et al., 2005; 2006; 2010; Wang et al., 2008; Li et al., 2010; Gu and Zhou, 2009]. As a result, it is infeasible to apply such algorithms to large-scale real world data due to the expensive computational cost.

In order to tackle these difficulties, instead of solving the relaxed clustering problems as in Eqs. (2-4), we solve the original clustering problem similar to Eq. (1). Specifically, we constrain the factor matrices of NMTF to be cluster indicator matrices and minimize the following objective:

We call Eq. (5) as the proposed Fast Nonnegative Matrix Trifactorization (FNMTF) approach.

Note that, the orthonormal constraints on F and G are removed in our objective, as their purposes (unique solution and labeling approximation) are automatically accomplished by the new constraints. Surprisingly, with these new constraints, though more stringent, as shown theoretically shortly in this section and empirically later in Section 4, the computational speed of our approach can be significantly improved.

2.2 An efficient optimization algorithm

Following the standard optimization procedures, we alternately solve the three variables F, S and G of J5 in Eq. (5).

First, fixing F, G, and setting the derivative of J5 with respect to S as zero, we have

Second, when F and S are fixed, the optimization problem to obtain G can be decoupled and we solve the following simpler problem for each i (1 ≤ i ≤ n):

Because gi· (1 ≤ i ≤ n) ∈ Ψ1×c is a cluster indicator vector in which one and only one element is 1 and the rest are zeros, the solution to Eq. (7) can be easily obtained by:

where ˜F = FS and ˜f·k is the k-th column of ˜F. Note that, Eq. (8) simply enumerates the c vector norms and seeks the maximum one, without involving any matrix multiplication.

Finally, when G and S are fixed, the optimization problem to obtain F can be similarly decoupled and we solve the following simpler problem for each j (1 ≤ j ≤ d):

Again, since fj· (1 ≤ i ≤ d) ∈ Ψ1×m is a cluster indicator vector for feature side, the solution to Eq. (9) is:

where ˜GT = SGT and ˜gl· is the l-th row of ˜GT .

The procedures to solve J5 are summarized in Algorithm 1. Due to the nature of alternating optimization, Algorithm 1 is guaranteed to converge to a local minima (existing NMF algorithms [Ding et al., 2005; 2006; 2010] also converges to a local minima because the objectives J3 and J4 are not convex in both variables F and G), and the proof is skipped due to space limit. As can be seen, in step 2 and step 3, the solutions are obtained by enumerating the vector norms, which is definitely much faster than the matrix multiplication used in the existing NMF methods. Upon solution, G gives the clustering results of data points, and F gives the clustering results of features directly.

Algorithm 1: Algorithm to solve J5 in Eq. (5).

Input: Data matrix X = [x·1, . . . , x·n] ∈ Rd×n. Initialize G ∈ Ψn×c and F ∈ Ψd×m with arbitrary class indicator matrices;

repeat

1. calculate S by Eq. (6) ;

2. calculate G by Eq. (8) ;

3. calculate F by Eq. (10) ;

until converges;

Output: Indicator matrices G for data point

clustering and F for feature clustering.

3 Locality preserved FNMTF (LP-FNMTF)

Recent researches showed that many real world data are sampled from the nonlinear manifolds which are embedded in the high dimensional ambient space [Belkin and Niyogi, 2002]. However, similar to traditional NMF and NMTF, the proposed FNMTF approach assumes that the data points and features are sampled from Euclidean spaces, and fails to discover the intrinsic geometrical and discriminative data and feature structures. Therefore, we further develop our FNMTF approach and propose the Locality Preserved FNMTF (LP-FNMTF) approach to enforce two geometrically based regularizers from both data and feature sides.

3.1 Manifold regularization and the optimization objective of LP-FNMTF method

Because we co-cluster an input data matrix on both data and feature dimensions, we consider two undirected graphs, one constructed from data points, denoted as Gd, and the other one from features, denoted as Gf . The corresponding affinity matrices Wd and Wf could be either computed from the input data matrix X (e.g., as in [Gu and Zhou, 2009]), or obtained from prior knowledge. According to manifold assumption [Belkin and Niyogi, 2002], the regularization terms to measure the smoothness with respect to the intrinsic manifolds of data points and features are given by [Cai et al., 2008; Gu and Zhou, 2009]:

where Ld = I - D12 d WdD12 d is the the normalized graph Laplacian of Gd, and Dd is the degree matrix of Gd; similarly, Lf = I - D12 f WfD12 f , andDf is the degree matrix of Gf . Incorporating Eq. (11) into Eq. (5), the objective of LPFNMTF approach is to minimize:

where α and β are regularization parameters to balance the reconstruction error of co-clustering in the first term and labeling smoothness in the data point space and feature space in the second and third terms, respectively.

Because F and G are constrained to be cluster indicator matrices, it is difficult to solve Eq. (12) in general. Hence we simplify this problem using the following proposition.

Proposition 1 Given a symmetric matrix A and its eigendecomposition A = PΣPT , where Σ ∈ Rc×c is a diagonal matrix with diagonal elements as the c largest eigenvalues, and P is the corresponding eigenvector matrix, the following two optimization problems are equivalent:

where Q is an arbitrary orthonormal matrix and

Proof. (P1) is equivalent to maxC∈Ψ tr (CT AC) that is further equivalent to minC∈Ψ ‖CCT - A‖2 . By definition the low-rank approximation of A is given by A = BQ(BQ)T , thus (P1) becomes minC∈Ψ,QT Q=I ‖CCT - BQ(BQ)T ‖2 . C approximating BQ is equivalent to CCT approximatingBQ(BQ)T . Hence, solving (P1) in Eq. (13) can be reasonably transformed to solve (P2) in Eq. (14), which completes the proof of Proposition 1. ■

Applying Proposition 1 in Eq. (12), the objective of our LP-FNMTF approach is transformed to minimize:

where Bd and Bf are computed from Ld and Lf following the procedures described in Proposition 1.

3.2 Optimization algorithm

Again, we use the alternating iterative method to solve Eq. (16). We first introduce the following theorem.

Theorem 1 Let H = BT C and the Singular Value Decomposition (SVD) of H be given by H = UΛVT . Fixing C and B, the optimumQ to problem (P2) defined in Eq. (14) is given by Q = UVT .

Proof. When C is fixed, problem (P2) in Eq. (14) is equivalent to maxQT Q=I tr (QT H).

We have tr (QT H) = tr (QUΛVT ) = tr (ΛVT QU) = tr (ΛZ) = ∑i λiizii, where Z = VT QU and λii and zii are the (i, i)-th entry of Λ and Z, respectively.

Note that Z is orthonormal, i.e., ZT Z = I, thus zii ≤ 1. On the other hand, λii ≥ 0 as λii is singular value of H. Therefore, tr (QT H) = ∑i λiizii ≤ ∑i λii, and when zii = 1 (1 ≤ i ≤ c), the equality holds. That is to say, tr (QT H) reaches its maximum when Z = I. Recall that Z = VT QU, the solution to maxQT Q=I tr (QT H) or (P2) is Q = UZT VT = UVT . Theorem 1 is proved. ■

Now, we solve Eq. (16). First, fixing F and G, by setting the derivative of J7 with respect to S as 0, we obtain

Second, fixing F, G and S, we can decouple Eq. (16) into two following subproblems:

Applying Theorem 1, Qd = UdVT d where Ud and Vd are obtained by SVD on BTd G; Qf = UfVT f whereUf andVf are obtained by SVD on BTf F.

Third, we fix S, F and Qd to update G. Because G is a cluster indicator matrix, Eq. (16) is decoupled to the following simpler problems for each 1 ≤ i ≤ n:

where ˜F = FS, ( ˜bd)i· denotes the i-th row of ˜Bd = BdQd. Thus, the solution can be obtained by

(20) Finally, whenfixingS, GandQf , let ˜Bf = BfQf , ˜GT = SGT and ˜gl· is the l-th row of ˜GT , we similarly obtain F as:

The procedures to solve J7 are summarized in Algorithm 2. Again, because step 4 and step 5 only involve vector norm enumeration without matrix multiplication, our algorithm is more computationally efficient. Empirical results show that the convergence of our algorithm is fast, which make it feasible to solve the large-scale real world problems via our approach in practice.

Algorithm 2: Algorithm to solve J7 in Eq. (16).

Input: Data matrix X = [x·1, . . . , x·n] ∈ Rd×n. 1. Initialize G ∈ Ψn×c and F ∈ Ψd×m with arbitrary class indicator matrices;

1. Calculate Bd and Bf fromLd and Lf following the description of Proposition 1;

repeat

1. Calculate S by Eq. (17) ;
2. Calculate Qd = UdVT d where Ud andVd are obtained by SVD on BTd G ;
3. Calculate Qf = UfVT f whereUf andVf are obtained by SVD on BTf F ;
4. Calculate G by Eq. (20) ;
5. Calculate F by Eq. (21) ;

until Converges;

Output: Class indicator matrices G and F for data and feature clustering tasks, respectively.

4 Experiments

In this section, we evaluate the proposed FNMTF and LPFNMTF approaches, and compare them against state-of-theart (co-)clustering methods, including Semi-NMF (SNMF) [Ding et al., 2010], Orthogonal NMTF (ONMTF) [Ding et al., 2006], Graph regularized NMF (GNMF) [Cai et al., 2008]

and Dual Regularized Co-clustering (DRCC) [Gu and Zhou, 2009] methods. We also report the clustering results by Kmeans and NMF [Lee and Seung, 2001] methods as baselines.

In our experiments, we use 5 data sets to evaluate the compared methods, which are summarized in Table 1. The first four are widely used as benchmarks in clustering literatures [Ding et al., 2006; Cai et al., 2008; Gu and Zhou, 2009; Ding et al., 2010]. The last one has very large sample size and feature size [Chen et al., 2010]. In order to run the experiments on contemporary computers, for RCV1 data set, following previous studies, we remove the keywords (features) appearing less than 100 times in the corpus, which results in 1979 (out of 47236) keywords in our experiments.

To evaluate the clustering results, we adopt three widely used standard metrics: clustering accuracy [Cai et al., 2008], normalized mutual information (NMI) [Cai et al., 2008] and cluster purity [Ding et al., 2006].

4.1 Clustering results

Experiments setup. Because each clustering algorithm has one or more parameters to be tuned, in order to compare fairly, we run these algorithms under different parameter settings, and select the best average result for each one. We set the number of clusters as the true number of classes for all clustering algorithms on all data sets.

For co-clustering methods, including ONMTF, DRCC and our two methods, the number of feature clusters is set to be the same as that of data clusters, i.e., m = c.

For manifold regularized methods, including GNMF, DRCC and our LP-FNMTF methods, we construct nearestneighbor graph following [Gu and Zhou, 2009], where the neighborhood size for graph construction is set by searching the grid of {1, 2, . . . , 10}, and the regularization parameters (i.e., α and β in Eq. (16)) are set by searching the grid of {0.1, 1, 10, 100, 500, 1000}.

Given the number of clusters, no parameter selection is needed for K-means, NMF and SNMF methods.

Results. Under each parameter setting of every method in comparison, we repeat clustering 50 times, and the average result is computed. We report the best average result for each method on five data sets in Table 2.

From the results in Table 2, we can see that the proposed FNMTF and LP-FNMTF methods consistently outperform the other compared methods, sometimes very significantly, which demonstrate the advantage of our approaches in terms of clustering performance. A more careful examination on the results shows that, the co-clustering methods, including ONMTF, DDRC and our proposed FNMTF and LP-FNMTF methods, generally achieve better clustering results, which are consistent with the widely accepted hypothesis that clustering of features can help clustering of data points. Finally, LP-FNMTF method is superior to FNMTF method on all the data sets except CSTR. This indicates that exploiting the geometric structures in data and feature spaces indeed can improve the cluster performance, which verifies manifold assumption and confirms the correctness of our algorithms.

4.2 Studies of computational speeds

In this subsection, we evaluate the computational speeds of the compared (co-)clustering methods. All our experiments are performed on a Dell PowerEdge 2900 server, which has two quad-core Intel Xeon 5300 sequence CPU processors at 3.0 GHz and 48G bytes memory.

We first examine the convergence rates of the compared methods. We repeat clustering 50 times by each method with its optimal parameters on each data set. The average iteration numbers of each method on each data set are reported in Table 3. The results show that the proposed FNMTF and LPFNMTF approaches require much less iterations to converge, therefore they are more computationally efficient.

In addition, we also report the average convergence time of the compared methods on WebKB data and RCV1 data sets as in Figure 1. From the results, we can see that our FNMTF and LP-FNMTF methods are only slower than K-means method while much faster than all other state-of-the-art clustering methods. These results are consistent with our theoretical analysis that our methods are implemented on subproblems with much smaller sizes and use much less matrix multiplications. Therefore, our approaches are suitable for clustering on large-scale data. The results on three other smaller data sets are not shown due to space limit, from which the same observations can be seen.

5 Conclusions

In this work, we proposed a novel Fast Nonnegative Matrix Tri-factorization (FNMTF) method to simultaneously cluster both data side and feature side of an input data matrix. We adopt the idea of NMF/NMTF based co-clustering methods, but constrain the factor matrices to be cluster indicator matrices, a special type of nonnegative matrices. Through the new constraints, the optimization problem of our method is decoupled into a number of much smaller subproblems that require much less matrix multiplication than the existing NMF based co-clustering algorithms, which makes our approaches of particular use for real world large-scale data. We further developed our method to incorporate manifold information and proposed Locality Preserved FNMTF (LP-FNMTF) method. We conducted extensive experiments on benchmark data sets that demonstrate promising results of our methods, which is consistent with our theoretical analysis.

Acknowledgments

This research was supported by NSF-CNS 0923494, NSF-IIS 1041637, NSF-CNS 1035913.