Learning a Robust Consensus Matrix for Clustering Ensemble via Kullback-Leibler Divergence Minimization

1 Introduction

Clustering is a fundamental problem in machine learning. According to [Wang et al., 2009], traditional single clustering algorithms usually suffer from the robustness problems, because: (1) different clustering methods may discover very different structure in a given data set due to their different objective functions; (2) for a single clustering method, since no ground truth is available, we can hardly validate the clustering results; (3) some methods (such as K-means) are highly depend on their initializations. To improve the quality of clustering, the idea of ensemble has been proposed.

Clustering ensemble provides a framework for combining multiple base clusterings of a data set to generate a consensus clustering [Topchy et al., 2004]. In past decades, many clustering ensemble methods have been proposed. [Strehl and Ghosh, 2003] formalized clustering ensemble as a combinatorial optimization problem in terms of shared mutual information. [Topchy et al., 2003] also used information theoretic method to combine clusterings. [Fern and Brodley, 2004] proposed a graph cut method. [Li et al., 2007; Li and Ding, 2008] and [Du et al., 2011] applied non-negative matrix factorization (NMF) to clustering ensemble. [Wang et al., 2009] learned the consensus clustering results by minimizing the Bregman divergence over all input clusterings. [Wang et al., 2011] applied a Bayesian method to clustering ensemble. [Du et al., 2013] proposed a self-supervised framework for clustering ensemble. These methods try to learn the consensus clustering results by taking advantage of diversity between base clusterings and reducing the redundance in clustering ensemble. In this paper we treat clustering ensemble with a new perspective. The inputs of clustering ensemble task are several weak base clustering results. Since the clusters identified by these base clusterings are often imperfect, these base clustering results may not be fully reliable. Thus these imperfect base clusterings can be regarded as intrinsic clustering corrupted with "noises" and "outliers". As a result, it is necessary to recover these contaminated results for consensus clustering.

∗Corresponding author

Most existing clustering ensemble methods blindly combine multiple base clusterings of data sets without taking into account noises and outliers, thus incurring the robustness problem, i.e., their performance would be severely degraded by noises and outliers. In this paper, we propose a novel Robust Clustering Ensemble (RCE) method, which explicitly characterizes the noises in each clustering, and uses them to get a robust and consensus clustering. In detail, given a set of input clusterings for ensemble we first construct the connective matrices, where each entry indicates the probability of two instances belonging to the same class. Since the input clusterings may contain noises and outliers, these connective matrices may also be contaminated. We introduce a sparse and symmetric error matrix for each connective matrix to explicitly identify the noises, and integrate the error matrices and connective matrices into ensemble framework. Considering that the connective matrices have a clear probabilistic interpretation, we use the Kullback-Leibler divergence for consensus measuring. We further impose a low-rank constraint on the final consensus matrix to obtain a more clear cluster structure. The resulting optimization problem turns out to be hard to solve due to the involvement of the noise matrices, divergence function and the low-rank constraint. To solve the objective function, we develop a block coordinate descent algorithm which can be theoretically guaranteed to converge.

Main contributions of our work are summarized as follows

• To improve the robustness of clustering ensemble, we introduce the sparse and symmetric error matrices to char-

acterize the noises in each input clustering. To learn the final robust, low-rank consensus matrix, we minimize the disagreements among the connective matrices using the Kullback-Leibler divergence.

• We propose a block coordinate descent algorithm to solve the complex objective function which involves Kullback-Leibler divergence, sparse term and low-rank term.
• The experiments on several benchmark data sets show that our method outperforms other compared algorithms, which indicates the importance of robustness for clustering ensemble.

2 Clustering Ensemble

Let X = {x1, x2, ..., xn} be a set of n data points. Suppose we are given a set of m clusterings C = {C 1 , C 2 , ..., Cm} of the data in X, each clustering C i consisting of a set of clusters {πi1 , πi2 , ..., πik }, where k is the number of clusters and X = ∪kj=1πij . Note that the number of clusters k could be different for different clusterings.

From C, we can construct symmetric connective matrix A(i) for partition C i as:

More generally, sometimes we can also get a soft connective matrix A(i) , where A(i) pq denotes the possibility that xp and xq belong to the same cluster. Thus A(i) pq ∈ [0, 1] for i from 1 to m. The task of clustering ensemble is to learn a consensus matrix ˜A fromA(1) ,...,A(m) , where ˜ Apq denotes the consensus probability that xp and xq belong to the same cluster.

3 Robust Clustering Ensemble

In this section, we present the framework of our RCE method and then discuss how to solve it.

3.1 Formulation

Since A(1) ,...,A(m) may be contaminated by noises and outliers, we denote E(i) as the sparse error matrix of the i-th connective matrix. The error matrix E(i) should be symmetric because the connective matrix is symmetric. Given E(i) , we obtain a cleaned connective matrix A(i) - E(i) . To learn the consensus matrix, we minimize the average disagreement over the cleaned connective matrices. Since both ˜ Apq and A(i) pq -E(i) pq have a probabilistic interpretation, we minimize the Kullback-Leibler divergence between ˜ Apq andA(i) pq -E(i) pq instead of the Euclidean distance. Thus, we obtain the following optimization problem:

where KL(·) is Kullback-Leibler divergence and λ1 is a balancing parameter. ‖ · ‖1 is the ℓ1-norm, which makes E(i) sparse as [Du and Shen, 2013; He et al., 2014] did. The constraints on ˜ Apq and A(i) pq -E(i) pq ensure that the cleaned connective matrices and consensus matrix have a probabilistic interpretation. The constraint on E(i) makes it symmetric.

Additionally, to make the final consensus matrix ˜A have a clear structure for clustering, we further require ˜A to be low-rank. Here we use the nuclear norm of ˜A to approximate the rank of ˜A inspired by [Liu et al., 2010; Luo et al., 2012; Liu et al., 2013]:

where λ2 is another balancing parameter, and ‖ · ‖∗ is the nuclear norm.

Since ˜ Apq and A(i) pq -E(i) pq are the possibility of xp and xq belonging to the same cluster, ˜ Apq and A(i) pq - E(i) pq follow the Bernoulli distribution. Thus we can expand the KL(·) and rewrite Eq.(2) to the following objective function (here 0log0 = 0):

Observe that due to the explicit characterization (E(i) ) of sparse, symmetric and bounded noises in clusterings, the above framework of clustering ensemble is robust. Moreover, to be more suitable for clustering ensemble task, we also introduce the Kullback-Leibler divergence and low-rank term in our approach.

3.2 Optimization

Eq.(3) involves two groups variables (˜A and E(i) ), thus we present a block coordinate descent scheme to optimize it. In particular, we optimize the objective with respect to one variable while fixing the other variables. This procedure repeats until convergence.

To handle the nuclear norm term ‖ ˜A‖∗, we borrow the following result from [Grave et al., 2011]:

Lemma 1. LetM∈ Rn×m. The nuclear normofMis equal to:

and the infimum is attained for S = (MMT )1/2 , where tr(·) is the trace of a matrix, and S ⪰ 0 means S is positive semidefinite.

According to Lemma 1, we rewrite Eq.(3) as:

Optimize S by fixing ˜A and E(i)

When fixing ˜A and E(i) , Eq.(4) is re-written as:

As suggested by [Grave et al., 2011], we need to add a term µtr(S-1 ). Otherwise, the infimum over S could be attained at a non-invertible S, leading to a non-convergent algorithm. Here µ is a small parameter and fixed to 0.001 for simplicity. According to Lemma 1, the infimum over S is then attained for:

where I is an identity matrix. Considering that when updating other variables, we just need S-1 instead of S, we calculate S-1 directly. More specifically, we compute VDiag(σk)VT as the eigenvalue decomposition of ˜A˜AT , and then calculate S-1 as S-1 = VDiag(1/√σk + µ)VT .

Optimize ˜A by fixing S and E(i)

By fixing E(i) and S, Eq.(4) can be simplified as:

For the constraint 0 ≤ ˜ Apq ≤ 1, we cannot get the closedform of ˜A. To solve Eq.(7), we use the auxiliary function approach [Lee and Seung, 2000]. We first introduce the definition and lemma of auxiliary function.

Definition 1. [Lee and Seung, 2000] Z(h, h′ ) is an auxiliary function for F(h) if the conditions

are satisfied.

Lemma 2. [Lee and Seung, 2000] If Z is an auxiliary function for F, the F is non-increasing under the update

To find an auxiliary function for Eq.(7), we first denote

C is a matrix whose (p, q)-th element is Cpq. We further introduce C+ = (|C| + C)/2, C- = (|C| - C)/2, D+ = (|D| + D)/2 and D- = (|D| - D)/2. Then following Theorem defines an auxiliary function of Eq.(7).

Theorem 1. Let J(˜A) be the objective function of Eq.(7), then the following function

is an auxiliary function of J( ˜A).

Proof. See Appendix A in the supplementary material.

According to Lemma 2, we minimize Z( ˜A, A′ ) instead of J(˜A). Let ∂Z(˜A,A′ ) ∂ ˜A = 0. We get the updating rule of ˜A:

where ˜ A(t+1) pq expresses the value of ˜ Apq in the (t + 1)-th iteration, and

It is obvious that ˜ Apq ≥ 0 always holds when updating ˜ Apq by Eq.(10). When ˜ Apq is greater than 1, since Z(˜A,A′ ) is convex w.r.t. ˜A, Z( ˜A, A′ ) is a monotone decreasing function in interval [0, 1]. Thus the optima of ˜ Apq is at ˜ Apq = 1. To sum up, we update ˜ Apq as follows:

Theorem 2. Updating ˜A by Eq.(11) will monotonically decrease the value of the objective of Eq.(7); hence it converges.

Proof. By Lemma 2 and Theorem 1, we can get that J(˜A(1) ) = Z(˜A(1) , ˜A(1) ) ≥ Z(˜A(2) , ˜A(1) ) ≥ J(˜A(2) ) ≥ .... So J(˜A) is monotonically decreasing. Since J(˜A) is obviously bounded below, we prove this theorem.

Optimize E(i) by fixing S and ˜A

When S and ˜A are fixed, the remaining problem can be further decomposed into n × n × m sub-problems, where only E(i) pq is involved. Considering the symmetry of A(i) and E(i) , there are two cases, i.e. off-diagonal and diagonal. For the i-th connective matrix, the sub-problem of the off-diagonal element E(i) pq (p ̸= q) is:

The sub-problem of the diagonal element E(i) jj is similar:

Since Eq.(13) can be regarded as half of Eq.(12), we only need to solve Eq.(12); and the solution of Eq.(13) is similar.

For simplicity, we denote a = ˜ Apq, b = ˜ Aqp, c = A(i) pq and x = E(i) pq . Theorem 3 gives the solution to Eq.(12).

Theorem 3. Eq.(12) is equivalent to

and the solution of Eq.(14) is

where

Proof. It is obvious that Eq.(14) is equal to Eq.(12). Nowwe solve Eq.(14). Since there is an absolute value function in Eq.(14), we consider two cases: x < 0 and x ≥ 0.

When x < 0, setting the derivative of Eq.(14) w.r.t. x to zero, we get:

To handle the constraint in Eq.(14), we introduce Lemma 3:

Lemma 3.

and

satisfies at any value of λ1 > 0 when a, b ∈ [0, 1].

Proof. See Appendix B in the supplementary material.

According to Lemma 3, to satisfy the constraint in Eq.(14), we can only take the second solution. Although the first solution sometimes can reach the boundary value c-1, it is easy to verify that when x = c-1, the objective function becomes infinitely large. Moreover, if 2λ1c-λ1-1+√(λ1+1)2-2λ1(a+b) 2λ1 > 0, which means the objective function is monotonically decreasing in [c - 1, 0], thus the optima is 0. To sum up, if x < 0, then the optima is:

Similarly, when x ≥ 0, the optima is:

Therefore, to optimize x, we compute two candidate values according to Eq.(17) and Eq.(18), then we choose the one which makes objective function of Eq.(14) smaller.

Algorithm 1 summarizes the whole optimization process. After getting ˜A, we set W = ( ˜A+ ˜AT )/2, and apply spectral clustering [Ng et al., 2001] onW to get the final clustering.

Algorithm 1 Robust Clustering Ensemble

Input: multiple connective matrices A(1) , ...,A(m) , parameters λ1, λ2 and µ.

Output: consensus matrix ˜A, error matrices E(i) .

• 1: Initialize ˜A= 1m ∑mi=1 A(i) .
• 2: while not converge do
• 3: Compute E(i) according to Theorem 3.
• 4: Compute VDiag(σk)VT as the eigenvalue decomposition of ˜A˜AT .
• 5: Set S-1 =VDiag(1/√σk +µ)VT .
• 6: Update ˜A by Eq.(11).
• 7: end while

3.3 Convergence and Complexity Analysis

According to Theorem 2, when updating ˜A, the objective function decreases monotonically. When updating E(i) , we find the global solution of this sub-problem which also makes the objective function decrease. In addition, the objective function has a lower bound. Thus Algorithm 1 converges.

In each iteration, when updating S-1 , the time complexity is O(n3 ) due to the eigenvalue decomposition, where n is the number of instances. When updating ˜A, since there is a matrix multiplication of two n×n matrices (D+ ˜A), the complexity is O(n3 ). When computing E(i) , which only contains element-wise operation, the complexity is O(n2m), where m is the number of connective matrices. Thus the overall complexity is O((n3+n2m)l), where l is the number of iterations.

4 Experiments

In this section, we evaluate the effectiveness of the proposed RCE method by comparing with several state-of-the-art clustering ensemble methods on benchmark data sets.

4.1 Data Sets

We use totally 11 data sets to evaluate the effectiveness of our proposed RCE, including images, texts, and UCI data sets. Data sets from different areas serve as a good test bed for a comprehensive evaluation. The basic information of these data sets are summarized in Table 1.

4.2 ComparedMethods

we compare RCE with the following algorithms:

• K-means, which is randomly initialized and the results are averaged over 200 independent runs.
• KC, which represents the results of applying K-means to a consensus similarity matrix, which is often used as a baseline in clustering ensemble methods such as [Li and Ding, 2008].
• Cluster-based Similarity Partitioning Algorithm (CSPA)[Strehl and Ghosh, 2003], which signifies a relationship between objects in the same cluster and can thus be used to establish a measure of pairwise similarity.
• HyperGraph Partitioning Algorithm (HGPA)[Strehl and Ghosh, 2003], which approximates the maximum

mutual information objective with a constrained minimum cut objective.

• Meta-CLustering Algorithm (MCLA)[Strehl and Ghosh, 2003], the objective of integration is viewed as a cluster correspondence problem.
• Nonnegative Matrix Factorization based Consensus clustering (NMFC)[Li and Ding, 2008], which uses NMF to aggregate clustering results.
• Generalized Weighted Cluster Aggregation (GWCA)[Wang et al., 2009], which learns the consensus clustering results by minimizing the Bregman divergence over all the input clusterings.
• Bayesian Clustering Ensemble (BCE)[Wang et al., 2011], which is a Bayesian model for ensemble.
• Robust Multiview Spectral Clustering (RMSC)[Xia et al., 2014]. Although it is not designed for clustering ensemble task, it also contains mechanism to handle noises and can be easily adopted to deal with clustering ensemble. We conduct the probabilistic transition matrices with the connective matrices, and then apply RMSC on them to obtain final clustering results.

4.3 Experiment Setup

Following the similar experimental protocol in [Wang et al., 2011], we run k-means 200 times with different initializations to obtain 200 base clustering results, which are divided evenly into 10 subsets, with 20 base results in each of them. Clustering ensemble methods are then applied on each subset.

The number of clusters is set to be the true number of classes for all data sets and algorithms. We fix the parameter µ to 0.001 when updating S and tune both the parameters λ1 and λ2 from [10-4 , 104 ] by grid search. We tune the parameters in compared methods as suggested in their papers. To measure the performance of clustering, the average clustering Accuracy (ACC) and Normalized Mutual Information (NMI) on the 10 subsets are reported. To validate the statistic significance of results, we also calculate the p-value of t-test.

4.4 Experimental Results

Table 2 shows the clustering results. Bold font expresses that the difference is statistically significant (p-value of t-test is less than 0.05). From Table 2, it can be seen that most clustering ensemble methods perform better than K-means, which indicates the benefit of ensemble methods. It can also be seen that, our proposed method RCE shows superior performance gains over the baselines w.r.t. ACC and NMI on most of the 11 data sets, which demonstrates that our robust ensemble can improve the performance of clustering.

RMSC was introduced for multivew spectral clustering. It has a mechanism to handle noises in multiple views, but differs significantly from ours. First, it is based on probabilistic transition matrix instead of connective matrix; second, it considers only a sparse error matrix, while we consider a symmetric, sparse and bounded matrix; it seems that our characterization of noises is more suitable for ensemble clustering tasks, as demonstrated in our experiments.

4.5 Parameter Study

We explore the effect of the parameters on clustering performance. There are two parameters in our method: λ1 and λ2. We tune these two parameters from [10-4 , 104 ]. We show the ACC and NMI on JAFFE and ORL data sets and the results are similar on other data sets. Figure 1 shows the results, from which we see that the performance of our method is stable across a wide range of the parameters.

5 Conclusion

In this paper we proposed a unified framework for robust clustering ensemble. We introduced symmetric and sparse error matrices to characterize noises and integrated them into a robust framework to learn a low-rank consensus matrix. We presented a block coordinate descent algorithm to solve the induced hard optimization problem and proved its convergence. Finally, experiments on benchmark data sets demonstrated the effectiveness of our method.

As ongoing work we are considering methods for handling connective matrices with missing values.

Acknowledgments

This work is supported in part by China National 973 program 2014CB340301 and NSFC grant 61379043, 61273291.