Towards Robust Co-Clustering

1 Introduction

Clustering is a fundamental topic in unsupervised machine learning. To group similar data points into clusters, numerous clustering algorithms have been developed, e.g., Kmeans, spectral clustering and NMF [Seung and Lee, 2001]. However, these methods only consider data side clustering and omit the cluster structure on feature side. Typically, the clustering of samples and features are closely related in real world applications.

Motivated by the duality and interdependence between data and feature clusters, several co-clustering approaches have been proposed to partition both dimensions simultaneously by modeling the relationship between samples and features, such as graph-based method [Dhillon, 2001], informationtheoretic based method [Dhillon et al., 2003] and matrix factorization based method [Ding et al., 2006]. More recently, it has been shown that the intersample and interfeature relationships also play important roles in co-clustering task [Gu and Zhou, 2009; Zhang et al., 2012b]. [Gu and Zhou, 2009; Shang et al., 2012] proposed dual regularized (semi)-GNMTF models, which impose graph regularization on both data and feature cluster assignment matrices, to explore the geometric structure of data manifold and feature manifold. [Gu et al., 2011] further imposed a normalized cut like constraint to make GNMTF model well-defined. [Wang et al., 2011] proposed a fast GNMTF method by constraining factors to be cluster indicator matrices. Besides, GNMTF method has also been developed to address other aspect of co-clustering [Long et al., 2012].

Although GNMTF based co-clustering methods have attracted increasing attention in recent years, they still surffer from the following problems in real-world applications. First, these methods adopt the squared loss, which is well-known to be unstable with respect to outliers and noises, to measure the reconstruction error incurred by modeling the relationship between sample and features. Accordingly, a few noisy entries with large errors may easily dominate the factorization. Second, the intersample and interfeature relationships, such as the regularization on data and feature graphs in GNMTF, are also essential for co-clustering task. When the data set is corrupted by outliers and noises, the constructed graphs may also be contaminated [Zhang et al., 2012a] and fail to capture the desired geometrical information, for instance, two samples may be mis-linked and incurred inappropriate weight in a nearest neighbor graph. As a result, the graph regularization may mislead the factorization and deteriorate the performance.

In this paper, a novel Robust Co-Clustering (RCC) method is presented by extending GNMTF towards algorithmic robustness to both data outliers and unreliable graphs. Specifically, to model the sample-feature relationship faithfully, we introduce a sparse error matrix to explicitly capture the grossly corrupted entries in the data reconstruction function. With this error matrix, a much cleaned data could be recovered and a robust factorization result could be expected. Besides, to robustly model the intersample and iterfeature relationships, we apply the `1-norm error function, instead of the traditional `2 function, to alleviate the influence of large and unexpected regularization incurred by unreliable graphs. In this way, the error matrix in data reconstruction function alleviates the influence of noisy data, and the `1-norm reduces the impact of unreliable graph regularization errors and leads to sparse and faithful regularization. Hence, RCC is robust to noisy data and unreliable graphs. To solve the optimization problem of our approach, we derive an alternating iterative algorithm and prove its convergence. We also show the connection between the error matrix in data reconstruction function and the Huber M-estimator in robust statistics [Hampel et al., 2011]. Thus, the robustness of our data reconstruction function can be interpreted as reducing the influence of large reconstruction errors which are often caused by outlier entries. Experimental results on real data sets show that RCC consistently outperforms the other algorithms in terms of clustering results, and possess robustness to noisy data and unreliable graphs.

2 Preliminaries

Throughout the paper, we use small letters (e.g., a) to denote scalars, capital letters (e.g., A) to denote matrices. We denote Ai·↪A·j↪ and Aij as the i-th row, the j-th column and the entry at the i-th row and j-th column of a matrix A. We use tr(A) to denote the trace of a square matrix A, ⊙ to denote the element-wise product and ↵ for element-wise division. We also define 1d 2 Rd as a column vector with all its elements equals to 1.

Given a nonnegative data matrix X 2 Rd⇥n where rows correspond to features and columns to samples, the goal of co-clustering is to group the features ¶X1·↪ ▷ ▷ ▷ ↪Xd·♦ into k1 clusters ¶C 0i ♦k1 i=1 , while group the data points ¶X·1↪ ▷ ▷ ▷ ↪ X·n♦ into k2 clusters ¶Ci♦k2 i=1 .

2.1 Graph Regularized NMTF

To co-cluster data points and features simultaneously, [Gu and Zhou, 2009; Shang et al., 2012] proposed graph regularized (semi-)NMTF models by imposing graph regularization on data and feature graphs under the manifold assumption, which can be formulated as minimizing the following objective function

where F and G are feature and data point cluster assignments matrices, H reflects the association between feature clusters and data clusters, LF = DF -WF is the graph Laplacian on the feature graph, DF and WF are the degree matrix and the adjacency matrix of the feature graph. The definition of LG for data graph is the same as above, λF and λG are regularization parameters.

3 Robust Co-Clustering

3.1 Formulation

Motivated by the recent development in robust principle component analysis [Cands et al., 2011], we propose a novel nonnegative matrix tri-factorization model to approximate the data matrix, where we assume some entries of the data matrix may be arbitrarily corrupted, but the corruption is sparse. We introduce an error matrix S 2 Rd⇥n to explicit capture the sparse corruption. Thus, the goal of robust tri-factorization is to approximate the nonnegative data matrix X as

where F 2 Rd⇥k1 ↪H 2 Rk1⇥k2 and G 2 Rn⇥k2 are constrained to be nonnegative, and S is required to be sparse. Due to the presence of S, a much clean data matrix can be recovered from sparse corruption, and more faithful factorization result could be expected.

Besides, the intersample and interfeature relationships, such as the dual graph regularizers in GNMTF, are also essential for clustering task [Zhang et al., 2012b]. However, when the data set is corrupted by outliers and noises, the constructed graph may be contaminated and not reliable to capture the desired geometrical information. For instance, two data points are mis-linked in a KNN graph due to outliers. As a result, the graph regularization may mislead the factorization and deteriorate the performance. Considering two far away samples or features are forced to be close according to unreliable graph regularization, we further propose to minimize the `1-norm regularization error which can suppress the large and unexpected regularization errors, inspired by recent developments in graph embedding [Zhang et al., 2012a]. As a result, the affect of large regularization errors caused by unreliable graphs can be reduced. By minimizing the `1-norm of graph regularization error, we can also re-estimate more faithful graph structure such as the weights between samples or features as shown later.

By combining robust factorization and robust dual graph regularization into a unified framework, our Robust CoClustering (RCC) method can be formulated as follows

where λS↪ λF and λG are regularization parameters. With the `1 normalization constraints on each row of F and G, and the `1 norm of regularization errors in the last two terms, the optimization problem in Eq. (2) is well defined and no longer suffers from either scale transfer problem or trivial solution in GNMTF [Gu et al., 2011].

3.2 Algorithm

The problem in Eq. (2) is not convex in all variables together, but convex in them separately. In the following, we introduce an iterative algorithm based on block coordinate descent to minimize the objective function in Eq. (2). We separately update the value of S↪ F↪H↪ and G while holding the other variables as constant. Thus, a local minima can be expected by solving a sequence of convex optimization problems. We also prove the convergence of the algorithm.

Update S

Optimizing Eq. (2) w.r.t. S is equivalent to minimizing

It is easy to show that the above problem is element-wise decoupled, and the unique solution of each subproblem admits a closed form called the soft-thresholding operator [Bach et al., 2012], i.e.,

where Eij = (X - FHGT )ij. By substituting Eq. (4) into Eq. (3), we get

It is interesting to note that the right part of Eq. (5) is often denoted as Huber M-estimator in robust statistics [Hampel et al., 2011]. Based on the duality between Eq. (3) and the Huber loss in Eq. (5), it can be seen that: if λS > 2maxij ♣Eij♣, Eq. (3) degenerates to ordinary NMTF; if λS ! 0, Eq. (3) behaves like `1-norm factorization; Eq. (3) behaves like `2-norm on small errors (Eij  λS◁2) and like `1-norm on large errors by a trade off parameter λS. Different from traditional methods [Shang et al., 2012; Gu et al., 2011] which use the squared loss on all reconstruction errors, our method automatically adapts the `1-norm on large errors, which are often caused by outliers. Thus, the influence of outliers with large errors can be alleviated. Besides, the reformulation in Eq. (5) can also be used to guide the selection of the regularization parameter λS.

Update F

Optimizing Eq. (2) w.r.t. F is equivalent to minimizing

To deal with the nonsmoothness of the `1-norm on graph regularization, we first apply the quadratic variational form for the `1-norm [Bach et al., 2012], so that we can rewrite Eq. (6) in a compact form:

Note that although the `1-norm constraints on latent factors (F and G) have been explored in NMTF based methods in previous works [Zhuang et al., 2011; Long et al., 2012], the constrained optimization problem is usually solved by a two-step strategy, which first minimizes an unconstrained problem without considering the constraints and then uses an

where f W F ii0 = W F ii0 2♣♣Fi·-Fi0· ♣♣2 , eDF ii = Pi0 f W F ii0 , P = (X - S)GHT , Q = HGT GHT , and we introduce P + ij = (♣Pij♣ + Pij)◁2, P - ij = (♣Pij♣ -Pij)◁2.

additional `1 normalization operator on each row of F and G in each iteration. However, this is not a principled way to solve the problem, and the price of normalization is potential performance decrease [Sun and Hamme, 2012].

The auxiliary function [Gu and Zhou, 2009] of the above objective function can be constructed as follows

where the solution of F derived from e F always satisfies the `1 normalization constraints.

Let A = (P - + F tQ + λF eDF F t ) ↵ F t , C = (P + + λF f W F F t ) ⊙F t , then Eq. (10) is equivalent to

Note that only F is involved in the Eq. (13). Summing Eq. (13) over k, we have

It is obvious that the solution of Eq (14) always satisfies the `1 normalization constraint for almost all λF .

In this paper, we use Algorithm 1 developed in [Sun and Hamme, 2012] to compute the fixed point of Eq. (13), where only element-wise operations are involved and a few iterations are needed to converge (i.e, nIter = 10).

Algorithm 1 Algorithm to compute the fixed point of Eq.(13)

Update H

Optimizing Eq. (2) w.r.t. H is equivalent to minimizing

Introducing the Lagrangian multiplier ⇤H 2 Rk1⇥k2 and setting the partial derivative to zero, we obtain

Using the Karush-Kuhn-Tucker condition (⇤H)ijH2 ij = 0, we get the following updating rule

Update G

Optimizing Eq. (2) w.r.t. G is equivalent to minimizing

Considering the duality between F and G, the minimization of Eq. (18) w.r.t. G can also be solved by the algorithm developed for optimizing F in Eq. (6).

In summary, we present the alternating iterative algorithm for optimizing Eq. (2) in Algorithm 2.

Algorithm 2 Algorithm to solve the problem in Eq. (2)

Input: The data matrix X, the initial values of F↪H↪G, data and feature graphs WF ↪WG, regularization parameters λS↪ λF ↪ λG.

Output: Partition matrices F and G.

repeat

Update S according to Eq. (4);

3.3 Convergence Analysis

In this subsection, we prove the convergence of the proposed algorithm by using the properties of auxiliary function developed in [Seung and Lee, 2001] and the following Lemma provided in [Nie et al., 2010].

Lemma 1. For any non-zeros vectors a↪ b 2 Rd , the following results follows [Nie et al., 2010]

Theorem 1. The Algorithm 2 monotonically decreases the objective function of RCC in Eq. (2) in each iteration and converges to an optimal solution.

Proof. Since the objective function of RCC in Eq. (2) is bounded below, we only need to prove the objective function monotonically decreases under the updates of St+1 ↪ F t+1 ↪Ht+1 ↪ and Gt+1 in each iteration. The solution of St+1 in Eq. (4) is the global minima of the subproblem in Eq. (3). Based on the properties of the auxiliary function [Seung and Lee, 2001], we have J4(F t ) ≥ Z(F t+1 ↪ F t ) ≥ Z(F t+1 ↪ F t+1 ) = J4(F t+1 ). Thus the objective function in Eq. (7) will monotonically decrease under the update of F t+1 , that is J4(F t+1 )  J4(F t ), i.e.,

So, we only need to prove the update of F t+1 satisfies J4(F t+1 )  J3(F t ). And we have the following inequal-

Recalling the result in Lemma 1, we know that

Combining Eq. (20) with Eq. (21), we have the following results

This inequality indicates that the objective function J3 in Eq. (6) will monotonically decrease in each iteration. Similarly, we can also prove the objective function J6 in Eq. (18) monotonically decreases under the update of Gt+1 . The nonincreasing of the update of Ht+1 in Eq. (17) can be proved according to the properties of the auxiliary function [Seung and Lee, 2001].

4 Experiments

In this section, we evaluate the effectiveness and robustness of our proposed RCC method for clustering. We compared RCC with the following clustering and co-clustering methods: Kmeans, NMF [Seung and Lee, 2001], ONMTF [Ding et al., 2006], GNMF [Cai et al., 2011], DRCC [Gu and Zhou, 2009], DNMTF [Shang et al., 2012], LDCC [Zhang et al., 2012b].

4.1 Parameters Setting

Since each clustering algorithm has one or more parameters to be tuned, in order to compare these algorithms fairly, we run these algorithms under different parameter settings and select the best average result for comparison. We set the number of clusters equal to the true number of classes for all the data sets and clustering algorithms. For graph regularized matrix factorization methods, the neighborhood size p to construct the graph is searched over the grid ¶1↪ 2↪ ▷ ▷ ▷ ↪ 10♦ and the binary weighting scheme is used for its simplicity. The regularization parameters are searched over the grid ¶0▷1↪ 1↪ 10↪ 100↪ 500↪ 1000♦ according to [Gu and Zhou, 2009]. For dual graph regularized methods, the regularization parameters on data and feature graphs are set to the same value. For tri-factorization methods, the number of feature clusters is set to be the same as the number of data clusters. For LDCC, the local regularization parameter is set to be 1 and the embedding dimension is searched over ¶5↪ 10↪ ▷ ▷ ▷ ↪ 50♦ according to [Zhang et al., 2012b]. Based on the duality between Eq. (3) and the Huber loss in Eq. (5), the regularization parameter λS is empirically set to λS = 2medianij(Eij). Since the reconstruction error Eij is changed in each iteration and we also adjust λS accordingly.

4.2 Data Sets

In our experiment, we have collected five public datasets to show the effectiveness of the proposed method. These data sets include one face image dataset, i.e., JAFFE, two hand written digit image datasets from UCI, i.e, MFEA and OPTDIGIT [Shang et al., 2012], two gene expression datasets, i.e., LUNG and GLIOMA, we use the processed data from [Nie et al., 2010]. Table 1 summarizes the characteristics of the data sets used in this experiment.

4.3 EvaluationMetrics

The result is evaluated by comparing the cluster label of each sample with the label provided by the data set. We use three metrics to evaluate the clustering performance. One metric is accuracy (ACC), which is used to measure the percentage of correct labels. The second metric is the normalized mutual information (NMI). In clustering applications, mutual information is used to measure how similar two sets of clusters are. The cluster purity is also used to measure the extent to which each cluster contained data points from primarily one class. The definition of these metrics can be found in [Ding et al., 2006; Gu and Zhou, 2009].

4.4 Clustering Results

Under each parameter setting of each method mentioned above, we independently repeat the experiments for 10 times and report the best average result in Table 2. We can see that RCC consistently outperforms the other algorithms in terms of ACC, NMI and Purity on all datasets. This may imply that real world datasets may contain noises and outliers to some extent due to the conditions under which the data sets are captured. By explicitly considering these unexpected corruptions, our method can improve the clustering results in general cases.

5 Conclusion

In this paper, we propose a novel Robust Co-Clustering (RCC) method by extending GNMTF towards robustness to noisy data and unreliable graphs. Specifically, we introduce a sparse outlier matrix into the data reconstruction function and adopt the `1 norm to measure the graph regularization errors. As a result, the proposed method can achieve robust factorization by approximating the cleaned data recovered from sparse outliers, and achieve robust regularization by suppressing the large regularization errors of unreliable graphs via `1 norm. One possible extension is to devise robust GNMTF model under the half-quadratic minimization framework which is helpful in general cases.

6 Acknowledgments

This work is supported in part by NSFC grant 60970045 and China National 973 project 2013CB329305.