Pseudo Supervised Matrix Factorization in Discriminative Subspace

1 Introduction

Data clustering is always a hot research topic that has been widely studied in various areas, such as document clustering [Cai et al., 2011a], gene selection [Jiang et al., 2004] and image segmentation [Shi and Malik, 2000]. Among them, Non-negative Matrix Factorization (NMF) and spectral clustering are two widely-used methods. NMF aims to factorize a matrix into two non-negative matrices whose product reconstructs the original data matrix. According to [Ding et al., 2010], the two matrices correspond to the cluster centroid and indicator, respectively. Thus we obtain the clustering result directly by the cluster indicator matrix without doing extra post-processing. Spectral clustering can adapt to a wider range of geometries and detect non-convex patterns and linearly non-separable clusters [Ng et al., 2002]. The general spectral clustering method needs to construct an adjacency matrix and calculate the Eigen decomposition value of the corresponding Laplacian matrix [Chung and Graham, 1997].

∗Corresponding author

In the past few years, many algorithms and frameworks have been proposed in order to improve the data clustering performance. For K-means methods, Liu [Liu et al., 2017] proposed a compressed K-means method for fast large-scale clustering. Besides, Shen [Shen et al., 2017] tried a sparse embedded algorithm to accelerate K-means clustering. For NMF-based algorithms, Cai [Cai et al., 2011b] kept the local geometry of the data in low dimensional space by proposing Graph Regularized NMF. Wang [Wang et al., 2018] incorporated the ordinal relations and proposed a novel ranking preserving NMF approach. To enhance the robust of matrix factorization, Ke [Ke and Kanade, 2005] replaced l2-norm with l1-norm to improve the robustness, but the l1-norm still cannot guarantee the feature rotation invariance. Jin [Huang et al., 2013] utilized l2,1-norm to factorize the data matrix. And then, Zhang [Zhang et al., 2017a] adopted l2,1-norm to learn the low-dimensional representation of the original data and showed its effectiveness. Lu [Lu et al., 2017] attempted to establish a connection between Linear Discriminant Analysis (LDA) and NMF in a supervised or semi-supervised way, but it can not be applied to data clustering.

In particular, the manifold regularization term has been combined with clustering models to improve the performance. Cai [Cai et al., 2011a] proposed a graph model capturing the local manifold geometry to address the underlying concepts which are consistent with the intrinsic manifold. To learn a better affinity matrix, Zhang [Zhang et al., 2017b] adopted the adaptive manifold regularization to matrix factorization. Wang [Wang et al., 2017] tried to learn a mapping in a relative low-dimensional with a more discriminative ability by a Grassmann manifold. What's more, Ma [Ma et al., 2018] extended the standard concept factorization model with an adaptive manifold regularizer which can represent of the raw data itself and a graph manifold regularizer which can reveal the local structure information of original data. And Zhang [Zhang et al., 2017a] incorporated the manifold regularization terms on both the low-dimensional feature representation and the cluster labels and get better local geometrical information. All of those models show that the manifold regularization term has an extraordinarily good clustering ability and can be expanded to other frameworks.

However, there are still some drawbacks in those existing methods. First, the framework based on NMF only factorizes matrices on high dimensional data space while ignores the intrinsic data information in the low-dimensional subspace. Thus, ordinary NMF needs more constraints to capture complicated structures. Second, local relationships among all data points got in the optimization process, such as pseudoinformation, cannot be used in most spectral clustering and manifold regularization methods. The geometry-metric structure of data distribution lacks an effective way to capture. Third, matrix factorization squares the residue error of each data point with a l2-norm objective function, so the clustering results are easily affected by the outliers.

In this paper, a novel unsupervised matrix factorization method, Pseudo Supervised Matrix Factorization (PSMF), is proposed for data clustering. The major contributions of this paper are summarized as follows:

1. A unified framework to cluster in the discriminant subspace is proposed, which combines LDA with NMF. So the intrinsic data structure in the low-dimensional subspace can be found and utilized.
2. A pseudo supervised manifold regularization term is proposed in order to find the subspace that discriminates different classes. This regularization term utilizes the pseudo-information to instruct itself and to refine the clustering results.
3. A novel Augmented Lagrangian Method (ALM) based optimization algorithm is designed to effectively and can efficiently seek the optimal solution of the problem.

The reminder of this paper is organized as follows: the derivation of our model is described in detail in Section 2. And then, the optimization algorithm is proposed in Section 3. After that, the experimental results on several benchmark datasets and further study are presented in Section 4, followed by the conclusion in Section 5.

2 ProposedModel

2.1 Non-negative Matrix Factorization

Given the data matrix X = [x1, x2, ..., xn], xi ∈ Rd×1 , among which d and n are the dimensionality and sample number, respectively. NMF approximates X with the product of two non-negative matrices:

where F ∈ Rd×c is the cluster centroid and G ∈ Rn×c is the cluster indicator matrix. Note that the product of δF and GT /δ results in the residue error when the scalar δ > 0. To get the unique solution of problem (1), Huang [Huang et al., 2013] imposed the orthogonal constraint on G. Problem (1) can be reformulated as:

where Ic ∈ Rc×c is the identity matrix. Then the optimal solution maintains its uniqueness.

2.2 Linear Discriminant Analysis

LDA tries to learn a linear projection matrix W ∈ Rd×m to project the d-dimensional data into the m-dimensional representation. Given a label matrix L = [l1, l2, ..., ln]T ∈ {0, 1}n×c , Yang [Yang et al., 2011] defined total-class scatter St, between-class scatter Sb and within-class scatter Sw:

where µ is the mean of all data points, µi is the mean of points in the i-th class, ni is the number of points in the i-th class. S = L(LT L)-1/2 is the scaled label matrix and H ∈ Rn×n is In - 1n 1n1Tn . (In is the n-dimensional identity matrix and 1n ∈ Rn×1 is a column vector with all its elements as 1)

One of the aims of LDA is to find the optimal W to make points from the same class closer to each other. Thus, the objective function can be written as:

where Im ∈ Rm×m is an identity matrix, and Tr(·) is trace operator. According to Eq.(3), problem (4) is equivalent to:

2.3 Matrix Factorization in Discriminative Subspace

For most of the applications nowadays, high dimensional data is difficult to capture its intrinsic structure but to count on its low-dimensional subspace. In order to find it, a connection between matrix factorization and LDA is proposed. Considering the framework of NMF, when G is fixed, the optimal F is computed as XG. Replace F with XG, then problem (2) becomes:

which is similar to problem (5). The difference between problem (5) and problem (6) is that X and G are replaced with W T XH and S, respectively. What's more, the optimal G for problem (2) is the cluster indicator matrix in the matrix framework while S in problem (5) is the scaled label matrix. In this way, G and S have similar meanings in practice. Therefore, the unsupervised matrix factorization framework is connected to LDA:

where G can be considered as the cluster indicator matrix or the pseudo-information matrix. What's more, W is the projection matrix that can find the discriminant subspace in the matrix factorization framework.

According to [Huang et al., 2013], we replace the Frobenius norm with the l2,1-norm to improve the robustness. In this way, we get a robust unsupervised framework:

2.4 Pseudo SupervisedManifold Regularization

Considering the original data matrix X ∈ Rd×n and its projected low-dimensional data matrix Y ∈ Rm×n , W ∈ Rd×m is the matrix which projects the d-dimensional data into the m-dimensional representation, i.e. Y =W T X.

LDA tries to find the subspace that discriminates different classes by minimizing the trace of the within-class scatter matrix SW while maximizing the trace of the between-class scatter matrix SB. According to [Zhang et al., 2009], to minimize SW , we need to solve the following problem:

where C is the number of classes; Ni is the number of samples in the i-th class; ⃗y(j) i is the j-th sample in the i-th class and ⃗ym i is the centroid of the i-th class. Problem (9) can be reduced into:

where

Yi = [⃗yi, ⃗yi1 , ..., ⃗yiNi-1 ] and ⃗eNi-1 = [1, ..., 1]T ∈ RNi-1 . LW = ∑Ni=1 SiLWi ST

i is the alignment matrix which can be obtained by an iterative procedure [Zhang and Zha, 2004]:

for i = 1, ..., N with the initialization LW = 0. With LWi and Eq.(11) [Zhang et al., 2009], we get:

According to the definition of LDA, we replace Y with the form of W T X. So problem (12) is equivalent to this problem:

Problem (13) is a manifold regularizer which is developed from LDA, so it needs data structure information to instruct itself.

In addition, we assume that any two points in high-density region of the low-dimensional manifold should share the same cluster in perspective of manifold learning. The lowdimensional data graph S ∈ Rn×n is built by W T X when W is fixed. Naturally, we minimize the following problem to get its geometry structure from low-dimensional subspace:

where Lm is the Laplacian matrix of the data graph S.

2.5 Objective Function

By combing the problem (8), (13) and (14) together, we finally get the objective function of the proposed Pseudo Supervised Matrix Factorization (PSMF) method:

where λ1 and λ2 are two non-negative manifold regularization parameters. Every time when W and G are updated, we reconstruct Lm and LW which are a Laplacian matrix and an alignment matrix for further updates.

3 Optimization

With the constraint W T StW = Im, problem (15) is difficult to solve. So we first disturb the diagonal elements of St by a scalar ϵ > 0 which is small enough, and then St is positive definite. We can decompose the constraint with Cholesky decomposition in the form of St = RT R. Thus, denoting RW as P , and denoting (R-1 )T XH and (R-1 )T X as A and B, respectively. Problem (15) is simplified into:

The above problem is not convex with three variables, so we propose to solve it with Augmented Lagrangian Multiplier (ALM) [Nie et al., 2015]. Since problem (16) depend on P and G, we introduce three auxiliary variables E = P T A - FGT , Z1 = G and Z2 = P . Then problem (16) is equivalent to the following ALM problem:

where µ ∈ R1×1 is the ALM parameter, and Λ1, Λ2 and Λ3 are ALM multipliers. Then we optimize each variable iteratively.

3.1 Update E

When fixing all the variables except E, we have:

where M = P T - FGT + Λ1 µ . According to [Huang et al., 2013], the optimal E is computed as:

3.2 Update Z1

When updating Z1, problem (17) becomes:

which can be further reduced into a closed-from problem:

where T = G+ Λ2 µ - λ1 µ LmG. Therefore, the optimal Z1 is:

3.3 Update Z2

To update Z2, problem (17) is reduced to:

By expanding the objective function and removing the irrelevant terms, we get:

where R = - λ2 µ BLWBT P + P + Λ3 µ . It is equivalent to:

And [Huang et al., 2013] proved that the optimal solution of the above problem is:

where U3 ∈ Rd×m and V3 ∈ Rm×m are the left and right singular vectors of the compact singular value decomposition of R.

3.4 Update G

To update G, problem (17) is reduced to:

Similar to Z2, the optimal solution of the above problem is:

where U1 ∈ Rn×c and V1 ∈ Rc×c are the left and right singular vectors of the compact singular value decomposition of K and K = (- λ1 µ LmZ1) + (P T A - E + Λ1 µ )T F +(Z1 - Λ2 µ ).

3.5 Update P

To update P , problem (17) is reduced to:

Similar to Z2, the optimal solution of the above problem is:

where U2 ∈ Rd×m and V2 ∈ Rm×m are the left and right singular vectors of the compact singular value decomposition of D and D = - λ2 µ BLWBT Z2 + A(FGT + E - Λ1 µ )T + (Z2 - Λ3 µ ).

3.6 Update F

Optimizing problem (17) with regard to F yields the following sub-problem:

Because GT G = Ic, the above problem is reformulated as:

Finally, the optimal F can be computed as:

3.7 Update µ,Λ1,Λ2 and Λ3

The ALM parameters are updated as follows:

where the parameter ρ controls the convergence speed.

With the updating rules discussed in this section, the optimization algorithm of problem (17) is summarized in Algorithm 1.

Algorithm 1 Algorithm to solve problem (17)

Input:

Original data matrix X ∈ Rd×n ;

Number of clusters k;

Regularization parameters λ and µ;

Output:

Cluster indicator G.

• 1: Initialize G ∈ Rn×c andW ∈ Rd×m;
• 2: Initialize F = W T XG and P = RW and compute H, A, B, St, Sb, Sw;
• 3: while not converge do
• 4: Update E by Eq.(19);

5: Update Z1 by Eq.(22);

• 6: Update Z2 by Eq.(26);
• 7: Update G by Eq.(28);
• 8: Update P by Eq.(30);
• 9: Update F by Eq.(33);
• 10: Calculate ALM parameters by Eq.(34);
• 11: end while

4 Experiments

In this section, the effectiveness of the proposed PSMF is demonstrated by nine real-world datasets. Parameter sensitivity of it is also discussed here.

4.1 Performance on Benchmark Datasets

In this part, we evaluate the performance of several methods on benchmark datasets to show the effectiveness of our algorithm on data clustering.

Datasets

There are in total nine datasets used in experiments, including one object image dataset, i.e. COIL20 [Cai et al., 2011b], two face image datasets, i.e. YALE [He et al., 2005] and UMIST [Wechsler et al., 2012], and six datasets from the UCI Machine Learning Repository, i.e. Dermatology, Movement, Scale, Iris, Automobile, and Lung-discrete. Table 1 summarizes the characteristics of the datasets used in our experiments.

Evaluation Metrics

Following [Liu et al., 2018], we adopt two widely used evaluation metrics to quantitatively measure clustering performance of our algorithm.

Clustering Accuracy (Acc) [Cai et al., 2005] discovers the one-to-one relationship between clusters and classes and measures the extent to which each cluster contains data points from the corresponding class. It is defined as follows:

where ri denotes the cluster label and li denotes the true class label, n is the total number of samples, δ(x, y) is the delta function that equals 1 if x = y and equals 0 otherwise, and map(ri) is the permutation mapping function that maps each ri to the equivalent label from the data set.

Normalized Mutual Information (NMI) [Est´evez et al., 2009] is used for determining the quality of clusters. Given a clustering result, it is estimated by:

where ni denotes the number of data contained in cluster Ci, ˆnj is the number of data belonging to class Lj, and ni,j denotes the number of data that is in the intersection between cluster Ci and class Lj .

Compared Algorithms

Six state-of-the-art clustering methods are taken for comparison, including K-means, Normalized Cut (NCut) [Shi and Malik, 2000], NMF [Lee and Seung, 1999], Graph Regularized NMF (GNMF) [Cai et al., 2011b], Robust Manifold NMF (RMNMF) [Huang et al., 2013] and Robust manifold Matrix Factorization (RMMF) [Zhang et al., 2017a]. For all the clustering methods, the number of clusters is known as input.

Initialization and Parameters Setting

For K-means, we use a faster Matlab method [Cai, 2011]. For NCut, the data graph is based on the Euclidean distances between two data points. For GNMF and RMNMF, the graph is constructed by finding the five nearest neighbors, in which every edge is weighted from 0 to 1. For NMF, GNMF will degrade to the ordinary NMF when the regularization parameter is 0. For RMMF and our method, the cluster indicator is initialized with the method in [Nie et al., 2014]. All the initial values of µ and Λ are set empirically since they have little influence on the clustering results. K-means, NCut, NMF, GNMF, and RMNMF are sensitive to the initialization, so we run 10 times and calculate the average results.

To compare these methods fairly, we run them with some selected parameter combinations and choose the best results for comparison. For GNMF and RMNMF, we set the regularization parameters α and µ by searching the grid of {10-5 , 10-4 , ..., 104 , 105 }. For RMMF, we set the number of iterations as 200 and choose the regularization parameters α and β by searching the grid of {10-5 , 10-4 , ..., 104 , 105 }. Finally, for the proposed PSMF method, we set the number of iterations as 50 and choose the regularization parameters λ1 and λ2 by searching the grid of {10-5 , 10-4 , ..., 104 , 105 }. Note that there is no parameter selection required for Kmeans and NMF since the number of clusters is given as input.

Clustering Results

Table 2 and 3 present the ACC and NMI comparison results of all included clustering algorithms on nine datasets. From two tables, we can easily observe that our method always achieves the highest ACC and NMI. In addition, we can find the following detailed points from the results:

1. Compared with other algorithms employing the idea of manifold regularization, such as GNMF, RMNMF, and RMMF, our algorithm shows a better performance both in the measurement of ACC and NMI, which suggests that our method can preserve the local geometrical structure embedded in the high dimensional space well. In other words, PSMF performs better on discovering the intrinsic geometrical and discriminative data structure for the clustering task.
2. The clustering performance of PSMF in datasets with high dimensionalities, such as COIL20, YALE, and Lung-discrete, is extremely better than other algorithms both in the measurement of ACC and NMI, which indicates that our method can effectively figure out the clustering problem of high dimensional data compared with other methods. What's more, PSMF can project the high-dimensional data into a low-dimensional subspace while maintains a good performance on clustering. It shows that PSMF can find a good intrinsic data representation in the low-dimensional subspace.
3. For datasets with categorical attributes, such as Dermatology, COIL20, and UMIST, the clustering performance of PSMF is also better than other algorithms in both two measurements. The reason for it may be that our method can extract comprehensive information of the original data, no matter what type of data is (numerical, categorical or mixed).

4.2 Parameter Sensitivity

In the proposed PSMF method, there are two manifold regularization parameters, i.e. λ1 and λ2, which determine the weight of two manifold regularizations. In order to study the influence of λ1 and λ2 on the clustering performance, we tune λ1 and λ2 in the same range of {10-5 , 10-4 , ..., 104 , 105 } and show the clustering accuracy of PSMF by a 3D visualizable way in Figure 1. According to Figure 1, it is easy to learn that the regularization parameters have much effect on clustering accuracy. Therefore, one of the important future work is to develop a regularization parameter setting rule to get the optimal parameter combination.

5 Conclusion

In this paper, we propose a novel Pseudo Supervised Matrix Factorization (PSMF) in discriminative subspace for data clustering. Different from other manifold regularized clustering methods, our method utilizes the pseudo-information to optimize the objective function iteratively. In addition, the intrinsic geometry structure can be captured in the discriminative subspace by PSMF. Furthermore, the l2,1-norm is introduced to enhance the robustness, so that our method is not sensitive to the data outliers. The proposed PSMF can be satisfactorily optimized by the suggested ALM-based method. Extensive experiments on multiple benchmark datasets illustrate that the proposed model outperforms other state-of-theart clustering algorithms. In future research, we may develop a regularization parameter setting rule for an optimal parameter combination.

Acknowledgements

This work was supported in part by the National Natural Science Foundation of China under Grants 61771349, 61822113 and 61772455, the National Key R&D Program of China under Grant 2018YFA0605501, the Natural Science Foundation of Hubei Province under Grants 2018CFA050 and 2018CFB432 and Yunnan Natural Science Funds under Grant 2018FY001(-013).