Subspace Clustering via New Low-Rank Model with Discrete Group Structure Constraint

1 Introduction

Subspace clustering is an important problem in machine learning, computer vision, such as motion segmentation [Costeira et al., 1997], image compression [Hong et al., 2005] and human face image clustering [Ho et al., 2003]. Subspace clustering assumes that the data points are drawn from multiple low dimensional subspaces. And the goal of subspace clustering is to find such multiple low dimensional subspaces such that all data points are segmented correctly and each group of data points fits into one of the low-dimensional subspaces. Due to the assumption that all data points lie in the multiple subspaces, traditional clustering methods such as Kmeans, spectral clustering are often difficult to cluster the data points with relatively high accuracy. The subspace clustering is suggested to solve such problem and has gained promising result in the previous research.

There has been a lot of work presenting many different subspace clustering methods. Generally speaking, these methods can be grouped into the following categories: factorization based method such as [Gruber and Weiss, 2004], algebraic based method such as [Vidal et al., 2005], statistic based method such as [Ho et al., 2003], spectral clustering based method such as [Elhamifar and Vidal, 2013]. Recently, sparse learning based methods such as Sparse Subspace Clustering (SSC) [Elhamifar and Vidal, 2013], Subspace Segmentation via Quadratic Programming (SSQP) [Wang et al., 2011] and low-rank based methods such as Robust Subspace Segmentation by Low-Rank Representation (Robust LRR) [Liu et al., 2013], Latent Low-Rank Representation for Subspace Segmentation and Feature Extraction [Liu and Yan, 2011] are all proved to get good performance. For the sparse coding model, it aims to represent one data point with a sparse linear or affine combination of other data points to achieve the subspaces efficiently and naturally. The low-rank based methods utilize low-rank representation opposed to sparse representation to recover the low-dimensional subspaces. The optimal subspace clustering assumptions include different data groups are assigned to different subspaces, there exist enough data points to shape each subspace, and no further subspaces can be derived from the current subspaces. In the other work [Yu and Schuurmans, 2011], the fast closed-form solutions were proposed to solve subspace clustering problem.

⇤To whom all correspondence should be addressed. This work was partially supported by US NSF-IIS 1117965, NSF-IIS 1302675, NSF-IIS 1344152, NSF-DBI 1356628, NIH R01 AG049371.

We propose a novel low-rank segmentation model to solve the above assumption based subspace clustering problem. Instead of learning the affinity matrix done by previous subspace clustering methods including SSC [Elhamifar and Vidal, 2013] and LRR [Liu et al., 2013], we propose to learn the group indicator directly such that the low-rank structure for each subspace is clear. We utilize the Schatten p-norm instead of trace norm to relax the rank constraint to better approximate the low-rank constraint, tactically avoid the integer programming problem imposed by the group indicator constraint, and extend our new model to discuss the general case in which the subspaces don't pass the original point. An efficient method is suggested with convergence prove to solve the proposed subspace segmentation model. Both synthetic and real world data sets demonstrate our proposed new subspace clustering model's effectiveness.

2 New Subspace ClusteringModel

To begin with, we denote the groups for which data points in X 2 Rd⇥n are clustered as ¶X1↪X2↪ ▷▷▷↪Xk♦. Here we first discuss the condition when each subspace passes through the origin. We will talk about the case when each subspace might not pass through the origin later.

In order to get the optimal subspace clustering result, we need seek the low-rank subspaces to accommodate different groups of data points. To achieve the optimal low-rank solutions, we can minimize the rank of all subspace simultaneously as:

However, the above objective often has trivial solution. If all data points are allocated into a single subspace, the objective value in Eq. (1) is min(d↪ n). If the ground-truth are k groups, denote ri = rank(Xi), the objective value Pki=1 ri is often larger than min(d↪ n). For example, there are 5 subspaces of dimension 3 in R10 , and 20 points per subspace are sampled. As a result, each data matrix Xi is 10 by 20 and of rank 3. The objective value is 15 (sum of ground-truth ranks). But, if all data points are allocated into a single subspace, the objective value is 10, which is even smaller than that of the ground-truth. Thus, the objective function in Eq. (1) is not suitable for the subspace clustering task.

Because in subspace clustering, ri << min(d↪ n), the value of Pki=1 r2 i is smaller than (min(d↪ n))2 in general cases. For example, in the above example, the sum of the square ranks Pki=1 r2 i = 45. If all data points are in a single subspace, (min(d↪ n))2 = 102 = 100 > 45. Thus, we can avoid the trivial solution. Based on this consideration, we propose to solve the following objective function:

We denote a cluster indicator matrix by Idx 2 Rk⇥n , Idx(i↪ j) = 1 if the j-th data point is grouped into the ith cluster or subspace, and Idx(i↪ j) = 0 otherwise. The k diagonal matrices G1↪G2↪ ▷▷▷↪Gk are defined, where the diagonal elements of Gi(1  i  k) are formed by the i-th row of matrix Idx. Identity matrix is denoted by I. Note that rank(Xi) = rank(XGi), the problem (2) is equivalent to the following problem:

where the summation constraint ensures each data point is assigned to one and only one group.

The objective function is non-convex, so does the discrete constraint. It is well known that trace norm (denote as k·k⇤ ) is the convex envelope of the rank function, and is thought to be the best convex approximation [Recht et al., 2010]. The discrete constraint for the group indication matrix Gi♣ki=1 ✓ ¶0↪ 1♦n⇥n can be relaxed to be Gi♣ki=1 ✓ [0↪ 1]n⇥n which is explained as probabilities for each data point belonging to certain group. As a result, based on the convex relaxation of both objective function and constraint, the problem (3) now becomes:

But such convex relaxation format may not be the best approximation to the original non-convex problem (3), because it includes two relaxation including trace norm and continuous constraint. So we want to use other relaxation which can approximate the original problem (3) better. The Schatten pnorm is a better relaxation for the rank objective function than the trace norm when p < 1 based on the definition [Nie et al., 2012; 2015]. As a result, in this paper we relax the rank based objective function to the following one with binary constraint for each Gi♣ki=1 :

where kFkpSp = Tr((FF T ) p2 ) is power p of the Schatten p-norm of matrix F .

In the next section, an iteration based re-weighted method [Nie et al., 2010; 2012; 2014] will be applied to solve it with convergence analysis.

3 Optimization Algorithm

According to the definition of Lagrange function, the Lagrangian function of problem (5) is

where g(⇤↪ Gi♣ki=1 ) encodes the constraints in problem (5).

Taking derivative w.r.t. Gi♣ki=1 and setting to zero, it becomes

where

Note that Di depends on Gi, we use an iteration based algorithm to get the solution that satisfies Eq. (6). First, according to Eq. (7), Di♣ki=1 can be calculated based on the current solution Gi♣ki=1 . If Di♣ki=1 is given, the optimal solution Gi♣ki=1 to the following problem will satisfy Eq. (6):

Then update the current solution Gi♣ki=1 by the optimal solution to the problem (8). It can be proved that this iterative procedure converges to a local optimum of the problem (5). Here is the intuition on the changes from Eq. (6) to Eq. (8). Eq. (6) which has the Lagrange multipliers is one of the KKT conditions w.r.t. the derivative of primal variables Gi, and after that we set Di fixed. The same KKT condition w.r.t. the derivative of primal variables Gi is used to Eq. (8). The objective derivative w.r.t. Gi in Eq. (8) is the same as the first part of Eq. (6). Because our constraint does not change, it has the same form as the second part for multipliers in Eq. (6). The reason why we introduce the multipliers in Eq. (6) is to show the connection between Eq. (6) and Eq. (8). Because the constraints are not changed and we want to use iteration based re-weighted method to solve this optimization problem instead of KKT condition involved method, the multipliers disappear.

Problem (8) can be rewritten as

where Ai = XT DiX. the problem (9) can be further changed to the following problem because of the discrete constraint Gi♣ki=1 ✓ ¶0↪ 1♦n⇥n :

Note that Gi♣ki=1 are n by n diagonal matrices, the problem (9) can be rewritten as the following problem:

where gci is the c-th diagonal element of matrix Gi and aci is the c-th diagonal element of matrix Ai.

The optimal solution to the problem (11) can be easily obtained as follows

The algorithm to solve problem (5) is summarized in Algorithm 1. As can be seen, the algorithm is very simple and concise.

Algorithm 1 Algorithm to solve problem (5).

Initialize Gi♣ki=1 such that the constraints in the problem (5) are satisfied

while not converge do

1. Calculate Di♣ki=1 = p kXGikpSp (XG2i XT ) p-2 2
2. Update Gi♣ki=1 , where the c-th diagonal element gci of Gi is updated by Eq. (12)
3. Calculate Ai♣ki=1 = XT DiX

end while

3.1 Convergence Analysis

We have the following theorem on the convergence of the proposed algorithm.

Theorem 1 The Algorithm 1 will finally converge to a local optimal solution of the problem (5) when p = 1.

Proof: Assume the updated Gi in step 3 of Alg. 1 is ˜Gi, Since ˜Gi is the optimal solution to the problem (8), we have

which can be written as

According to step 1 in Alg. 1, Eq. (13) can be written as

According to the Cauchy-Schwarz inequality, it can be proved that, when p = 1 we have

Therefore, based on Eqs. (14-15), we have

⇤

Eq. (16) indicates that the objective function in the problem (5) will monotonically decrease during the iteration until the algorithm converges. Since the objective in problem (5) has clearly lower bound 0, the algorithm will converge. When it converges, we can see that Eq. (6) will always be satisfied, so the algorithm will converge to a local optimum solution to the problem (5) when p = 1.

In practice, we observed the algorithm is also converged for 0 < p < 1. If the objective in problem (5) is changed to Pi (kXGikpSp )m(m > 2) (the Di in step 1 is changed to mp 2 (kXGikpSp )m-1 (XG2i XT ) p-2 2 accordingly), and/or the constraint is changed to the one as in Eq.(4) (the Eqs. (11-12) are changed accordingly), the convergence is also observed 1 .

3.2 Complexity Analysis

We suppose d < n in the following analysis. In step 1, we need to compute Di♣ki=1 . Computing kXGikpSp need SVD of XGi, which takes O(nd2 ). Suppose the SVD of XGi is XGi = U⌃V T , then (XG2i XT ) p-2 2 = U⌃p-2U T , which

1 Gi♣ki=1 = 1n is a trivial optimal solution to problem (4) because of the following facts: problem (4) is convex, Alg. 1 will converge to a local minimum, and the solution will be unchanged with this trivial solution as initialization in Alg. 1.

takes O(d3 ). So computing Di takes O(d2n) and computing Di♣ki=1 takes O(d2nk). Step 2 only needs to compute the diagonal elements of Ai, so step 2 takes O(d2nk). Computing Gi♣ki=1 in step 3 takes O(nk). In summary, the computational complexity of Algorithm 1 is O(d2nkt), where t is the iteration number. In the experiments we found the algorithm converges very fast, and always converges in 5-10 iterations. According to the above analysis, the time consumption of Algorithm 1 is linear w.r.t. the number of data. Therefore, the algorithm can easily handle large scale data if the dimensionality of data is not too high.

4 Subspace Clustering without Passing through Original Point

In the previous section, we assume that each subspace passes through the origin. If a subspace does not pass through the origin, the rank of the group will reduce one as long as the subspace is shifted to pass through the origin, as can be seen in Fig. 1. Therefore, the shifts can also be learned such that the ranks of the groups are minimized. In this case, instead of minimizing Eq. (2), we minimize the following objective function:

where 1 denotes a column vector with all ones. Similar to the previous section, minimizing Eq. (17) is equivalent to solving the following problem:

Similarly, the following problem for the subspace clustering problem in this case is proposed to solve:

This problem can also be efficiently solved by iteration based re-weighted algorithm. We define Xui = X - ui1T . Similarly to Eq. (8), we need to solve the following problem:

where

When fix ui♣ki=1 , the Gi♣ki=1 can be solved as in Eqs. (11-12) with Ai = XT ui DiXui .

When fix Gi♣ki=1 , the problem (20) becomes

The optimal solution to the problem (22) can be easily derived as

The algorithm to solve problem (19) is summarized in Algorithm 2.

Algorithm 2 Algorithm to solve problem (19).

Initialize Gi♣ki=1 such that the constraints in the problem (19) are satisfied. Initialize ui♣ki=1 = 0.

while not converge do

1. Calculate Xui ♣ki=1 = X - ui1T
2. Calculate Ai♣ki=1 = XT ui DiXui
3. Calculate Di♣k1 = p kXui GikpSp (Xui G2i XT ui ) p-2 2
4. Update Gi♣ki=1 , where the c-th diagonal element gci of Gi is updated by Eq. (12)
5. Update ui♣ki=1 = 1 1T G2i 1 XG2i 1

end while

5 Experiments

5.1 Experimental Results on Synthetic Data

First, we test our algorithm on three challenge Synthetic datasets, each of which comprises several groups distributed on different lower dimensional subspaces but with 5% level noises to deviate from the subspaces. These subspace intersect each other, so it is a very difficult subspace clustering problem.

We run our algorithm with ten different initializations and select the results with the best objective values. The results are shown in Fig. 2, which we can see that our algorithm can correctly find the low-dimensional structure and find the correct groups on these challenge datasets. The results indicate that the proposed method can effectively solve the subspace clustering problem and has the power to find the lowdimensional structure information hidden in data, which is also the goal in Generalized Principal Component Analysis (GPCA) [Vidal et al., 2005].

We also test our algorithm on a synthetic high-dimensional dataset. In this dataset, we randomly generate 250 data points with 50 dimensions lying in 5 different subspaces. Each subspace has 50 data points, and the dimensions of the 5subspaces are 5↪ 10↪ 15↪ 20↪ 25, respectively. In the data, we also add 5% level noises to deviate from the subspaces. The nearest neighbor graph of the data based on Euclidean distance is shown in Fig. 3(a), and the result Idx0 Idx of the proposed method is shown in Fig. 3(a), where Idx is defined in Section 2. This result also verify that the proposed method is an effective subspace clustering method.

5.2 Experiments on Real World Data

Here we just evaluate Algorithm 2, and leave out the Algorithm 1, because Algorithm 2 is the subspace clustering in the case which has more wide application. We evaluate the algorithm 2's performance in handling two real-world problems: motion segmentation and human images' clustering.

Motion Segmentation in Hopkins 155 Dataset

We first test our model on the Hopkins 155 motion dataset. This dataset consists of 155 sequences, which have about 39550 points tracked from 2 or 3 motions, separately. And in these 155 datasets, there are 120 datasets having two motions and 35 datasets having three motions. Our goal is to cluster these points into groups according to their motions in each sequence. Each sequence is a dataset, so there are 155 datasets with different properties on subspace's number, data samples' number. There are about 10 sequences that are grossly corrupted, having high error levels, although most of the outliers in the datasets have been removed manually.

We combine PCA and K-means methods to initialize our Gi(1  i  k). We use the PCA to project the coordinates in each sequence into the dimensions ranging from 5 to 20. Then we use K-means method to get our initialized Gi(1  i  k). We compare the performance of our algorithm with the best state-of-the-art motion segmentation algorithms: Sparse Subspace Clustering (SSC) [Elhamifar and Vidal, 2013], Subspace Segmentation via Quadratic Programming (SSQP) [Wang et al., 2011], Robust Low-Rank Representation (Robust LRR) [Liu et al., 2013], Groupwise Constrained Reconstruction (GCR) [Li et al., 2012], Groupwise Constrained Reconstruction with Dirichlet Process (GCRDP) [Li et al., 2012]. The final result with the known best results for each algorithm are reported in Table 2.

From the Table 2 we can observe that: Our proposed method outperforms other state-of-art methods in all accuracy and error's measurement including mean accuracy, minimum accuracy and standard deviation of error. The main reason is that our method takes the Schatten p- norm into account, which is the relaxation that better approximates the original problem (3) than the trace norm approximation when p < 1. Also, our method has considered the condition when the subspace does not pass through the origin to learn the shifts such that the ranks of the groups are minimized.

Subspace Clustering on Real-World Image Datasets

In order to check our method's suitability in the real human images' clustering condition, face datasets including JAFFE

[Lyons et al., 1998], MSRA [Liu et al., 2007], XM2VTS [XM2, ] and another human palm image dataset called PALM [Yan et al., 2007] are used.

We include Robust LRR [Liu et al., 2013], Convex LRR [Liu et al., 2013], its very similar model SSQP [Wang et al., 2011], SSC [Elhamifar and Vidal, 2013], and the classic Principle Component Analysis (PCA) to be our comparatives in the human images' clustering.

Similar to the motion segmentation method, we also combine PCA and K-means methods to initialize our Gi(1  i  k) . We use the PCA to project the coordinates in each dataset into the dimensions ranging from 10 to 50. Because the dataset for MSRA, XM2VTS, and PALM are relatively large, we first use PCA to project the original data samples to a low dimension based on the PCA ratio equaling to 0▷98 to facilitate the whole progress. The best results for each methods can been seen on Table 3.

From Table 3 we can find that: our proposed method performs better than other methods in four datasets, which proves that our method is more suitable for the subspace clustering problem than other proposed subspace clustering methods. Note that in the MSRA dataset, our method's clustering accuracy is 93▷11 %, which is the best result so far to the best of our knowledge. It is because the MSRA face samples are more satisfied to our proposed subspace clustering model than other 3 human image samples.

6 Conclusions

In this paper, a new subspace clustering model based on the Schatten p-norm is proposed, and it is a better approximation to the rank minimization problem than the trace norm approximation and others. According to the subspace model we propose, an iteration based re-weighted method with both effectiveness and efficiency is suggested to solve this model. The convergence of our proposed algorithm is proved. What's more, we consider the condition when the subspace does not pass through the origin to learn the shifts such that the ranks of the groups are minimized. Both the Synthetic datasets and Real world datasets demonstrate our method's effectiveness in dealing with subspace clustering problems.