A Joint Optimization Framework of Sparse Coding and Discriminative Clustering

1 Introduction

Clustering [Yang et al., 2014a] plays an important role in many real world data mining applications. To learn the hidden patterns of the dataset in an unsupervised way, existing clustering algorithms can be described as either generative or discriminative in nature. Generative clustering algorithms model categories in terms of their geometric properties in feature spaces, or as statistical processes of data. Examples include K-means and Gaussian mixture model (GMM) clustering [Biernacki et al., 2000], which assume a parametric form of the underlying category distributions. Discriminative clustering techniques search for the boundaries or distinctions between categories. With fewer assumptions being made, these methods are powerful and flexible in practice. For example, maximum-margin clustering [Xu et al., 2004], [Xu and Schuurmans, 2005], [Zhao et al., 2008] aims to find the hyperplane, that can separate the data from different classes with maximum margins. Information theoretic clustering [Li et al., 2004], [Barber and Agakov, 2005] minimize the conditional entropy of all samples. Many recent discriminative clustering methods have achievedsatisfactory performances [Zhao et al., 2008].

Moreover, many clustering methods extract discriminative features from data prior to clustering. The Principal Component Analysis (PCA) feature is a common choice but not necessarily discriminative [Zheng et al., 2011]. In [Roth and Lange, 2003], the features are selected for optimizing the discriminativity by Linear Discriminant Analysis (LDA). More recently, sparse codes prove to be both robust to noise and scalable to high dimensional data [Wright et al., 2009]. Furthermore, ℓ1-graph [Cheng et al., 2010] builds the graph by reconstructing each data point sparsely and locally with other data. A spectral clustering [Ng et al., 2002] is followed based on the constructed graph matrix. In [Sprechmann and Sapiro, 2010], dictionary learning is combined with the clustering process, which makes the use of Lloyds-type algorithms that iteratively re-assign data to clusters and then optimize the dictionary associated with each cluster. In [Zheng et al., 2011], the authors learned the sparse codes that explicitly consider the local data manifold structures. Their results indicate that encoding geometrical information will significantly enhance the learning performance. However, their clustering step is neither discriminative nor jointly optimized.

In this paper, we propose to jointly optimize feature extraction and discriminative clustering, in which way they mutually reinforce each other. We focus on sparse codes as the extracted features, and develop our loss functions based on two representative discriminative clustering methods, the entropy-minimization [Li et al., 2004] and maximum-margin [Xu et al., 2004] clustering, respectively. A task-driven bilevel optimization model [Mairal et al., 2012], [Wang et al., 2015] is then built upon the proposed framework. The sparse coding step is formulated as the lower-level constraint, where a graph regularization is enforced to preserve the local manifold structure [Zheng et al., 2011]. The clustering-oriented cost functions are considered as the upper-level objectives to be minimized. Stochastic gradient descent algorithms are developed to solve both bi-level models. Experiments on several popular real datasets verify the noticeable performance improvement led by such a joint optimization framework.

2 Model Formulation

2.1 Sparse Coding with Graph Regularization

Sparse codes prove to be an effective feature for clustering. In [Cheng et al., 2010], the authors suggested that the contribution of one sample to the reconstruction of another sample was a good indicator of similarity between these two samples. Therefore, the reconstruction coefficients (sparse codes) can be used to constitute the similarity graph for spectral clustering. ℓ1-graph performs sparse representation for each data point separately without considering the geometric information and manifold structure of the entire data. Further research shows that the graph regularized sparse representations produce superior results in various clustering and classification tasks [Zheng et al., 2011], [Yang et al., 2014c]. In this paper, we adopt the graph regularized sparse codes as the features for clustering.

We assume that all the data samples X = [x1, x2, · · · , xn], xi ∈ Rm×1 , i = 1, 2, · · · , n, are encoded into their corresponding sparse codes A = [a1, a2, · · · , an], ai ∈ Rp×1 , i = 1, 2, · · · , n, using a learned dictionary D = [d1, d2, · · · , dp], where di ∈ Rm×1 , i = 1, 2, · · · , p are the learned atoms. Moreover, given a pairwise similarity matrix W, the sparse representations that capture the geometric structure of the data according to the manifold assumption should minimize the following objective: 1 2 ∑ni=1 ∑nj=1 Wij||ai - aj||22 = Tr(ALAT), where L is the graph Laplacian matrix constructed from W. In this paper, W is chosen as the Gaussian Kernel: Wij = exp(- ||xi-xj||22 δ2 ), where δ is the controlling parameter selected by cross-validation.

The graph regularized sparse codes are obtained by solving the following convex optimization

Note λ2 > 0 is necessary for proving the differentiability of the objective function (see [7.1] in the Appendix). However, setting λ2 = 0 proves to work well in practice, and thus the term λ2||A||2F will be omitted by default hereinafter (except for the differentiability proof).

Obviously, the effect of sparse codes A largely depends on the quality of dictionary D. Dictionary learning methods, such as K-SVD algorithm [Elad and Aharon, 2006], are widely used in sparse coding literature. In regard to clustering, the authors in [Cheng et al., 2010], [Yang et al., 2014c], [Yang et al., 2014b] constructed the dictionary by directly selecting atoms from data samples. [Zheng et al., 2011] learned the dictionary that can reconstruct input data well. However, it does not necessarily lead to discriminative features. In contrast, we will optimize D together with the clustering task.

2.2 Bi-level Optimization Formulation

The objective cost function for the joint framework can be expressed by the following bi-level optimization:

where C(A,w) is a cost function evaluating the loss of clustering, with A as its input and w as the parameter. It can be formulated differently based on various clustering principles, two of which will be discussed and solved in Section 3.

Bilevel optimization [Yang et al., 2012] has been investigated in both theory and application sides. In [Yang et al., 2012], the authors proposed a general bilevel sparse coding model for learning dictionaries across coupled signal spaces. Another similar formulation has been studied in [Mairal et al., 2012] for general regression tasks.

3 Clustering-oriented Cost Functions

3.1 Entropy-Minimization Loss

Maximization of the mutual information with respect to parameters of the encoder model effectively defines a discriminative unsupervised optimization framework. The model is parameterized similarly to a conditionally trained classifier, but the cluster allocations are unknown [Barber and Agakov, 2005]. In [Dai and Hu, 2010], [Li et al., 2004], the authors adopted an information-theoretic framework as an implementation of the low-density separation assumption by minimizing the conditional entropy. By substituting the logistic posterior probability into the minimum conditional entropy principle, the authors got the logistics clustering algorithm, which is equivalent to find a labelling strategy so that the total entropy of data clustering is minimized.

Assuming K clusters, since the true cluster label of each xi is unknown, we introduce the predicted confidence probability pij that sample xi belongs to cluster j, i = 1, 2, · · · ,N, j = 1, 2, · · · ,K, which is set as the likelihood of multinomial logistic (softmax) regression:

The loss function for all data could be defined accordingly in a entropy-like form:

The predicted cluster label of ai is the cluster j where it achieves the largest likelihood probability pij. The logistics regression can deal with multi-class problems more easily compared with the support vector machine (SVM). The next important thing we need to study is the differentiability of (2).

Theorem 3.1. The objective C(A, w) defined in (4) is differentiable on D × w.

Proof: Denote X ∈ X, and D ∈ D. Also let the objective function C(A,w) in (4) be denoted as C for short. The differentiability of C with respect to w is easy to show, assuming the compactness of X, as well as the fact that C is twice differentiable.

We will therefore focus on showing that C is differentiable with respect to D, which is more difficult since A, and thus ai, is not differentiable everywhere. Without loss of generality, we use a vector a instead of A for simplifying the derivations hereinafter. In some cases, we may equivalently express a as a(D,w) in order to emphasize the functional dependence. Based on [7.1] in Appendix, and given a small perturbation E ∈ Rm×p , it follows that

where the term O(||E||2F ) is based on the fact that a(D, x) is uniformly Lipschitz and X ×D is compact. It is then possible

where g has the form given in Algorithm I. This shows that C is differentiable on D.□

Building on the differentiability proof, we are able to solve (1) using a projected first order stochastic gradient descent (SGD) algorithm, whose detailed steps are outlined in Algorithm 1. At a high level overview, it consists of an outer SGD loop that incrementally samples the training data. It uses each sample to approximate gradients with respect to w and D, which are then used to update them.

Convergence andComplexityAnalysis SGDconverges to stationary points under a few stricter assumptions than ones satisfied in this paper. A non-convex convergence proof assumes three times differentiable cost functions [Mairal et al., 2012]. As a typical case in machine learning, we use SGD in a setting where it is not guaranteed to converge in theory, but behaves well in practice.

Assuming n samples and dictionary size p, in each iteration of Algorithm 1, step 8 takes O(n) time. Step 4 is solved by the feature-sign algorithm [Lee et al., 2006], which is reduced to a series of quadratic programming (QP) problems. The computational bottleneck lies in solving the inverse of matrix DT D of size p × p, where applying the Gauss-Jordan elimination method takes O(p3 ) time per sample. Thus, Algorithm 1 takes O(np3 ) time per iteration, and O(Cnp3 ) in total (C is a constant absorbing epoch numbers, etc.). Further, if p is a constant, Algorithm I reaches O(Cn) time complexity.

3.2 Maximum-Margin Loss

Xu et al. [Xu et al., 2004] proposed maximum margin clustering (MMC), which borrows the idea from the SVM theory. Their experimental results showed that the MMC technique could often obtain more accurate results than conventional clustering methods. Technically, what MMC does is just to find a way to label the samples by running an SVM implicitly, and the SVM margin obtained would be maximized over all possible labelings [Zhao et al., 2008]. However, unlike supervised large margin methods which are usually formulated as convex optimization problems, maximum margin clustering is a non-convex integer optimization problem, which is much more difficult to solve. [Li et al., 2009] made several relaxations to the original MMC problem and reformulated it as a semi-definite programming (SDP) problem. The cutting plane maximum margin clustering (CPMMC) algorithm was presented in [Zhao et al., 2008] to solve MMC with a much improved efficiency.

To develop the multi-class max-margin loss of clustering, we refer to the classical multi-class SVM formulation in [Crammer and Singer, 2002]. Given the sparse code ai are the features to be clustered, we define the multi-class model:

where f j is the prototype for the j-th cluster and wj is its corresponding weight vector. The predicted cluster label of ai is the cluster of the weight vector that achieves the maximum value wT j ai. Let w = [w1, ...,wK], the multi-class maxmargin loss for ai could be defined as:

Note that different from training a multi-class SVM classier, where yi is given as a training label, the clustering scenario requires us to jointly estimate yi as a variable. The overall max-margin loss to be minimized is (λ as the coefficient):

But to solve (8) or (9) with respect to the same framework as logistic loss will involve two additional concerns, which needs to be handled specifically.

First, the hinge loss of the form (8) is non-differentiable, with only subgradients existing. That makes the objective function C(A,w) non-differentiable on D×w, and further the analysis in Theorem [3.1] proof can not be applied. We could have used the squared hinge loss or modified Huber loss for a quadratically smoothed loss function [Lee and Lin, 2013]. However, as we checked in the experiments, the quadratically smoothed loss is not as good as hinge loss in training time and sparsity. Also, though not theoretically guaranteed, using the subgradient of C(A,w) works well in our case.

Second, given that w is fixed, it should be noted that yi and ri are both functions of ai. Therefore, calculating the derivative of (8) over ai would involve expanding both ri and yi, and become quite complicated. Instead, we borrow ideas from the regularity of the elastic net solution [Mairal et al., 2012], that the set of non-zero coefficients of the elastic net solution should not change for small perturbations. Similarly, due to the continuity of the objective, it is assumed that a sufficiently small perturbation over the current ai will not change yi and ri. Therefore in each iteration, we could directly pre-calculate yi and ri using the current w and ai and fix them for ai updates 1 .

Given the above two handling, for a single sample ai, if the hinge loss is above 0, the derivative of (8) over w is:

where ∆ji denote the j-th element of the derivative for the sample ai. If the hinge loss is less than 0, then ∆ji = λwj i . The derivative of (8) over ai is wri -wyi if the hinge loss is over 0, and 0 otherwise. Note the above deduction can be conducted in a batch mode. It is then similarly solved using a projected SGD algorithm, whose steps are outlined in Algorithm 2. The convergence and complexity analysis is similar to Algorithm 1.

Algorithm 1 Stochastic gradient descent algorithm for solving (2), with C(A, w) as defined in (4)

Require: X, σ; λ; D0 and w0 (initial dictionary and classifier parameter); ITER (number of iterations); t0, ρ (learning rate)

• 1: Construct the matrix L fromX and σ.
• 2: FOR t=1 to ITER DO
• 3: Draw a subset (Xt,Yt) from (X,Y)
• 4: Graph-regularized sparse coding: computer A∗ :
• A∗ = argminA 1 2 ||X - DA||2F + λ ∑i ||ai||1 + T r(ALAT).
• 5: Compute the active set S (the nonzero support of A∗ )
• 6: Compute β∗ : Set β∗SC = 0 and β∗S = (DTS DS + λ2I)-1∇AS [C(A, w)]
• 7: Choose the learning rate ρt = min(ρ, ρ t0 t )
• 8: Update D and W by a projected gradient step:

where ∏w and ∏D are respectively orthogonal projections on the embedding spaces of w and D.

• 9: END FOR

Ensure: D and w

Algorithm 2 Stochastic gradient descent algorithm for solving (2), with C(A, w) as defined in (9)

Require: X, σ; λ; D0 and w0 (initial dictionary and classifier parameter); ITER (number of iterations); t0, ρ (learning rate)

• 1: Construct the matrix L fromX and σ.
• 2: Estimate the initialization of yi and ri by pre-clustering, i = 1, 2, ..., N
• 3: FOR t=1 to ITER DO
• 4: Conduct the same step 4-7 in Algorithm 1.
• 5: Update D and W by a projected gradient step, based on the derivates of (9) over ai and w (10).
• 6: Update yi and ri using the current w and ai, i = 1, 2, ..., N.
• 7: END FOR

Ensure: D and w

4 Experiments

4.1 Dataset and Evaluation

We conduct our clustering experiments on four popular real datasets, which are summarized in Table 1. We apply two widely-used measures to evaluate the performance of the clustering methods: the accuracy and the Normalized Mutual Information(NMI) [Zheng et al., 2011], [Cheng et al., 2010]. Suppose the predicted label of the xi is ˆ yi which is produced by the clustering method, and yi is the ground truth label. The accuracy is defined as:Acc = ✶Φ(ˆ

1 To avoid ambiguity, if yi and ri are the same, i.e., the max value is reached by two cluster prototypes simultaneously in current iteration, then we ignore the gradient update corresponding to ai.

where ✶ is the indicator function, and Φ is the best permutation mapping function [Lov´asz and Plummer, 2009]. On the other hand, suppose the clusters obtained from the predicted labels {ˆ yi}ni=1 and {yi}ni=1 as ˆ C and C, respectively. The mutual information between ˆ C and C is defined as:

where p(ˆ c) and p(c) are the probabilities that a data point belongs to the clusters ˆ C and C, respectively, and p(ˆ c, c) is the probability that a data point jointly belongs to ˆ C and C.

The normalized mutual information(NMI) is defined as:

where H( ˆ C) and H(C) are the entropies of ˆ C and C, respectively. NMI takes values between [0,1].

4.2 Comparison Experiments

Comparison Methods We compare the following eight methods on all four datasets:

• KM: K-Means clustering on the input data.
• KM + SC: A dictionary D is first learned from the input data by K-SVD [Elad and Aharon, 2006]. Then KM is performed on the graph-regularized sparse code features (1) over D
• EMC: Entropy-minimization clustering, byminimizing (4) on the input data.
• EMC + SC: EMC performed on the graph-regularized sparse codes over the pre-learned K-SVD dictionary D.
• MMC: Maximum-margin clustering [Xu and Schuurmans, 2005].
• MMC + SC: MMC performed on the graph-regularized sparse codes over the pre-learned K-SVD dictionary D.
• Joint EMC: The proposed joint optimization (2), with C(A,w) as defined in (4).
• Joint MMC: The proposed joint optimization (2), with C(A,w) as defined in (9).

All images are first reshaped into vectors, and PCA is then applied to reducing the data dimensionality by keeping 98% information, which is also used in [Zheng et al., 2011] to improving efficiency. The multi-class MMC algorithm is implemented based on the publicly available CPMMC code for two-class clustering [Zhao et al., 2008], following the multiclass case descriptions in the original paper. For all algorithms that involve graph-regularized sparse coding, the graph regularization parameter α is fixed to be 1, and the dictionary size p is 128 by default. For joint EMC and joint MMC, we set ITER as 30, ρ as 0.9, and t0 as 5. Other parameters in competing methods are tuned in cross-validation experiments.

Comparison Analysis All the comparison results (accuracy and NMI) are listed in Table. 2, from which we could conclude the following:

• 1: The joint EMC and joint MMC methods each outperform their "non-joint" counterparts, e.g., EMC + SC and MMC + SC, respectively. For example, on the ORL dataset, joint MMC surpasses MMC + SC by around 5% in accuracy and 7% in NMI. Those demonstrate that the key contribution of this paper, i.e., jointly optimizing the sparse coding and clustering steps, indeed leads to improved performances.
• 2: KM + SC, EMC + SC, and MMC + SC all outperform their counterparts using raw input data, which verifies that sparse codes are effective features that help improve the clustering discriminability.
• 3: The joint MMC obtains the best performances in all cases, outperforming the others, including joint EMC, with significant margins. The MMC + SC obtains the second best performance for the last three datasets (for ORL, it is joint EMC that ranks the second). The above facts reveal the power of the max-margin loss (9).

Varying the number of clusters On the COIL20 dataset, We re-conduct the clustering experiments with the cluster number K ranging from 2 to 20, using EMC + SC, MMC + SC, joint EMC, and joint MMC. For each K except for 20, 10 test runs are conducted on different randomly chosen clusters, and the final scores are obtained by averaging over the 10 tests. Fig. 1 shows the clustering accuracy and NMI measurements versus the number of clusters. It is revealed that the two joint methods consistently outperforms their non-joint counterparts. When K goes up, the performances of joint methods seem to degrade less slowly.

Initialization and Parameters As observed in our experiments, a good initialization of D and w can affect the final results notably. We initialize Joint EMC by the D and w solved from EMC + SC, and Joint MMC by the solutions from MMC + SC, respectively.

There are two parameters that we need to set empirically: the graph regularization parameter α, and the dictionary size p. The regularization term imposes stronger smoothness constraints on the sparse codes when α grows larger. Also, while a compact dictionary is more desirable computationally, more redundant dictionaries may lead to less cluttered features that can be better discriminated. We investigate how the clustering performances EMC + SC, MMC + SC, joint EMC, and joint MMC change on the ORL dataset, with various α and p values. As depicted in Fig. 2 and 3, we observe that:

• 1: While α goes up, the accuracy result will first grow up then go down (the peak is around α =1). That could be interpreted as when α is too small, the local manifold information is not sufficiently encoded. On the other hand, when α turns overly large, the sparse codes are "over-smoothened" with a reduced discriminability.
• 2: Increasing dictionary size p will first improve the accuracy sharply, which however soon reaches a plateau. Thus in practice, we keep a medium dictionary size p =128 for all experiments.

5 Conclusion

We propose a joint framework to optimize sparse coding and discriminative clustering simultaneously. We adopt graph regularized sparse codes as the feature to be learned, and design two clustering-oriented cost functions, by entropyminimization and maximum-margin principles, respectively. The task-driven bi-level optimization mutually reinforces both sparse coding and clustering steps. Experiments on several benchmark datasets verify the remarkable performance improvements led by the proposed joint optimization.

Figure 3: The clustering accuracy and NMI measurements versus the parameter choices of p.

6 Acknowledgement

We sincerely appreciate the great help from Associate Professor Linli Xu, affiliated with University of Science and Technology of China (USTC), for sharing with us her original implementation of multi-class MMC algorithm [Xu and Schuurmans, 2005].

7 Appendix

We recall the following lemma [7.1] in [Mairal et al., 2012]:

Theorem 7.1 (Regularity of the elastic net solution). Consider the formulation in (1) (we may drop the last term to obtain the exact elastic net form, without affecting the differentiability conclusions). Let λ2 > 0 , and X is assumed to be compact. Then,

• a is uniformly Lipschitz on X ×D
• Let D ∈ D, σ be a positive scalar and s be a vector in {-1, 0, 1}p . Define Ks(D, σ) as the set of vectors x satisfying for all j in {1, ..., p},

Then there exists κ > 0 independent of s,D and σ so that for all x ∈ Ks(D, σ), the function a is twice continuously differentiable on Bκσ(x) × Bκσ(D), where Bκσ(x) and Bκσ(D) denote the open balls of radius κσ respectively centered on x and D.