Robust Subspace Segmentation by Simultaneously Learning Data Representations and Their Affinity Matrix

1 Introduction

In scientific data analysis applications, we often have to face a set of data X = [X1, ...,Xn] ∈ Rd×n derived from a union of c subspaces {Ss}cs=1 , where d is the feature dimension and n is the number of data vectors. To characterize the given data as different groups such that the data in the same group are highly similar to each other (ideally drawn from one subspace), the subspace segmentation recently has been focus of considerable research in machine learning, computer vision and pattern recognition [Hong et al., 2006; Rao et al., 2010; Ho et al., 2003; Liu et al., 2013; Nie et al., 2009].

1.1 Notation Summary

Lowercase letters (u, v, ...) mean scalars and bold lowercase letters (u, v, ...) vectors. uj represents the j th entry in u. Bold uppercase letters (U,V , ...) stand for matrices. U T and U -1 are the transpose and inverse of U, respectively. Ui stands for the ith column of U, while Uij the j th element in the ith column of U. |Uij| is the absolute value of Uij. ‖U‖0, ‖U‖1, ‖U‖F and ‖U‖∗ denote the ℓ0 norm (number of nonzero entries), ℓ1 norm (∑i,j |Ui,j|), ℓ2 or Frobenius norm (√∑i,j U 2i,j ), and nuclear norm (sum of all the singular values) of U, respectively. ‖U‖2,0 and ‖U‖2,1 stand for the ℓ2,0 norm (∑i ‖√∑j U 2ij ‖0) and ℓ2,1 norm (∑i √∑j U 2ij ), respectively. 〈U,V 〉 is the inner product of two matrices with identical size, which is equal to the trace of U T V , i.e. tr(U T V ). U⊙V presents the Hadamard product of two matrices with identical size. Moreover, 0, 1 and I denote the vectors of all zeros, all ones and identity matrix with compatible sizes, respectively.

1.2 RelatedWork

Recently, many subspace segmentation methods have been proposed. From the perspective of their mechanisms of representing the subspaces, existing approaches can be roughly divided into four categories: iterative [Bradley andMangasarian, 2000], statistical [Ma et al., 2007], algebraic [Ma et al., 2008; Vidal et al., 2005] and spectral clustering methods [Lu et al., 2012; Elhamifar and Vidal, 2009]. An elaborate review of these methods can be found in [Vidal, 2010]. Our method belongs to the spectral clustering based one, therefore we review the related work along this direction in the following.

The key of spectral clustering based approaches is to construct a "good" affinity matrix A ∈ Rn×n , in which each element Aij reflects the similarity between data points Xi and Xj . Ideally, the affinity should be 1 if they are from the same cluster, 0 otherwise. Directly computing distances on the raw data (e.g. k-NN using cosine or heat kernel distances) is possibly the most intuitive way to conduct the data affinity matrix. Nie et al. [Nie et al., 2014] develop a more sophisticated method to learn the similarity matrix by adaptively assigning neighbors for each data point based on the local connectivity. But, the affinity matrix constructed on the raw data is unable to well reveal the global subspace structure of data.

Alternatively, inspired by the success of compressed sensing [Cand`es et al., 2006; Donoho, 2006], a large body of research on exploiting the relationship of data representations R ∈ Rn×n has been carried out [Elhamifar and Vidal, 2009; Lu et al., 2012; Liu et al., 2013; Lu et al., 2013; Saha et al., 2013; Feng et al., 2014], the formulation of which can be generally written as follows:

where E ∈ Rd×n denotes the residual, Ψ(R) stands for the regularizer on R, and λ is the coefficient controlling the importance of the regularizer. Θ(E) is the model of E, which can be with different forms depending on the characteristic of data. For instance ‖E‖1 is optimal for Laplacian residuals while ‖E‖2F for Gaussian. The choice of the model is important for the task of subspace segmentation, but not the main focus of this work. Hence, we simply adopt the ℓ2 for the rest of the paper as [Elhamifar and Vidal, 2009; Lu et al., 2012; 2013; Nie et al., 2014].

The main difference among the methods mentioned above lies in the regularization on R. Specifically, Sparse Subspace Clustering (SSC: Ψ(R) ≡ ‖R‖0) [Elhamifar and Vidal, 2009] introduces compressed sensing techniques to subspace segmentation. Because the ℓ0 norm is non-convex and NP-hard to approximate, replacing it with its tightest convex surrogate ℓ1 norm makes the problem tractable and gives the optimal solution to the original problem under some condition. The main drawback of SSC is that it processes data individually and thus lacks optimality due to the existence of inherent joint structure between the representations of data points. Instead of a sparse representation, Liu et al. [Liu et al., 2013] propose the Low Rank Representation (LRR: Ψ(R) ≡ rank(R)) to jointly find a low rank representation by minimizing the rank of R. As the rank(·) is also intractable to directly optimize, its convex replacement, the nuclear norm, is employed. Least Squares Regression (LSR: Ψ(R) ≡ ‖R‖2F ) [Lu et al., 2012] is a much more efficient solver for subspace segmentation than LRR with a similar grouping effect. To further refine the representation, Feng et al. [Feng et al., 2014] impose a block diagonal prior on the representation, which shows a reasonable improvement on the segmentation results. To simultaneously take into account the grouping effect and sparsity of representation, Grouping Sparse Coding (GSC: Ψ(R) ≡ ‖R‖2,0) [Saha et al., 2013] was developed. For efficiently solving the GSC problem, the ℓ2,0 norm needs to be convex relaxed to ℓ2,1 . Besides, Correlation Adaptive Subspace Segmentation (CASS: Ψ(R) ≡ ∑ni=1 ‖XDiag(Ri)‖∗) is designed [Lu et al., 2013] to better explore the subspace structure, which can be viewed as an adaptive interpolation between SSC and LSR. However, the computational load of CASS is relatively heavy as it involves a series of SVD operations for dealing with every single data.

Traditionally, after solving the problem (1), the representation is utilized to define the affinity matrix of an undirected graph in the way of |Rij|+|Rji| 2 for the data vectors Xi and Xj, then spectral clustering algorithms such as Normalized Cuts [Shi and Malik, 2000] are employed to segment the data into c clusters. Although this kind of affinity is somehow valid, the meaning of which is already not the same as the original definition. In this paper, we propose a method to construct a meaningful affinity matrix under the assumption that the data points should have a larger probability to be in the same cluster if their representations have a smaller distance.

1.3 Contribution

The contribution of this paper can be summarized in three aspects. 1) We propose a novel subspace segmentation method to jointly learn the representations of data and their affinity matrix in a unified optimization framework. 2) We design a new Augmented Lagrange Multiplier based algorithm to efficiently and effectively seek the solution of the associated optimization problem. 3) To demonstrate the efficacy and the superior performance of the proposed algorithm over the state-of-the-art alternatives, extensive experiments on both synthetic data and several datasets are conducted.

2 ProblemFormulation

Given a set of clean data points sufficiently sampled from c independent subspaces {S c s=1}, X = [X1, ...,Xn]. By introducing a hypothesized permutation matrix Γ that arranges the data to the true segmentation of data, we have X∗ = XΓ = [X1 , ...,Xc ] ∈ Rd×n , where Xs denotes a collection of ns data points from the sth subspace Ss with n = ∑cs=1 ns. The data can be self-represented by a linear combination of the items in X as X = XR. To avoid the trivial solution, we should impose some constraint on R. Considering the simplicity and the effectiveness leads us to choose the ℓ2 norm ‖R‖2F to do the job, any other choices can, of course, be selected. Recall our assumption that the data points should have a larger probability to be in the same cluster if their representations have a smaller distance, thus we naturally propose the following constraint:

where A ∈ Rn×n is the desired affinity matrix, Aij reflects the probability of the data points Xi and Xj from the same cluster based on the distance between their representations Ri and Rj . The constraints ATi 1 = 1 and Ai ⪰ 0 are to guarantee the probability property of Ai. However, simply solving the problem (2) results in that only the nearest representation (or equally the nearest data) is assigned as the neighbor of Ri (or Xi) with probability 1 and all the others with probabilities 0. Similar to R, we again enforce minimizing ‖A‖2F to prevent from the trivial solution.

Putting the concerns together with slight algebraic transformation gives the following formulation:

where LA is the Laplacian matrix of A, which is constructed in the way of DA - A. The degree matrix DA is defined as a diagonal matrix where the ith diagonal element is ∑j Aij . In addition, λ1, λ2 and λ3 are three non-negative weights balancing the corresponding terms.

In real world applications, the noise free and independent subspaces assumption may not be satisfied. It is desirable to extend the problem (3) to be robust to noises. With the introduction of the noise term ‖E‖2F , the problem can be finally formulized as follows:

In the next section, we will propose a novel algorithm to effectively and efficiently solve the problem (4).

3 Optimization

As we have seen in (4), it has combined all aforementioned priors and constraints for learning the representations of data and finding the meaningful affinity matrix with respect to the data representations in a unified optimization framework. Although the objective (4) is not jointly convex in A and R, but convex with respect to each of them when the other is fixed. The Augmented Lagrange Multiplier (ALM) with Alternating Direction Minimizing (ADM) strategy [Lin et al., 2011] has proven to be an efficient and effective solver of problems like (4). To apply ALM-ADM on our problem, we need to make our objective function separable. Thus we introduce one auxiliary variable Q to replace R in the trace term of the objective function (4). Accordingly, Q = R acts as the additional constraint. Note that the probability properties of every Ai are enforced as hard constraints. The augmented Lagrangian function of (4) L{∀i ATi 1=1;Ai⪰0} is

with the definition Φ(Z, C) ≡ µ2 ‖C‖2F +〈Z, C〉, where µ is a positive penalty scalar and, Z1 and Z2 are the Lagrangian multipliers. Besides the Lagrangian multipliers, there are four variables, including E, R, Q and A, to solve. The solver iteratively updates one variable at a time by fixing the others. The solutions of the subproblems are as follows:

E-subproblem: For computing E(t+1) , we take derivative of L with respect to E with the unrelated terms fixed and set it to zero, then obtain E(t+1) =:

where {µ(t) } is a monotonically increasing sequence. R-subproblem: It is obvious that all the terms related with R are quadratic, thus dropping the constant terms and taking derivative of L with respect to R gives R(t+1) =:

Q-subproblem: In a similar way to updating E and R, the closed form solution of this subproblem can be easily calculated by Q(t+1) =:

A-subproblem: The update of A(t+1) can be done via solving the following optimization problem:

Algorithm 1: Proposed Robust Subspace Segmentation

Input: Data matrix X ∈ RD×n , cluster number c, nearest neighbor number k, λ1 ≥ 0, λ2 ≥ 0 and

while not converged do

Segment the data into c groups by Normalized cuts.

Output: Final Data Segmentation

As can be seen from (9), it can be separated into a set of smaller independent problems, i.e.:

where dQ(t+1) i ∈ Rn×1 is a vector, the j th element of which is dQ(t+1) ij = λ2‖Q(t+1) i -Q(t+1) j ‖2F 4λ3 . For each Ai, the closed form solution is:

where the operator (u)+ turns negative elements in u to 0 while keeps the rest. Please notice that the parameter k ∈ {1, ..., n} is introduced to control the number of nearest neighbors of Qi (or Xi) that could have chance to connect to Qi (or Xi). In addition, the elements of ˜ dQ(t+1) i are those of dQ(t+1) i but with the ascending order. For clarity and completeness, the detailed proof of the closed form solution of (10) can be found in the appendix. As the graph constructed according to A obtained by (10) is generally an unbalanced digraph, we employ A+AT 2 to achieve the balance.

Multipliers: Besides, there are still two multipliers to update, which are simply done through:

The procedure of solving the problem (4) terminates when ‖X - XRt+1 - Et+1‖F ≤ δ‖X‖F with δ = 10-7 or the maximal number of iterations is reached. After obtaining the affinity matrix A and the representation matrix R, the spectral clustering techniques, Normalized Cuts for all the involved methods in this paper, can be executed to segment the input data into c clusters. The entire algorithm of subspace segmentation is summarized in Algorithm 1.

It is worth noting that although there is no established theory of global convergence in literature for ADM algorithms applied to non-convex problem as the one solved in this work, it is guaranteed that the proposed algorithm converges to at least a stationary point (first order optimality condition). In addition, empirical evidence on both synthesized and real data presented in the next section suggests that the proposed algorithm have very strong and stable convergence behavior.

4 Experimental Verification

In this section, we conduct experiments on synthetic data to reveal the efficacy of our proposed method, and on real data to demonstrate the superior performance of our method over the state-of-the-art alternatives including k-NN using heat kernel distance, CAN and PCAN [Nie et al., 2014], SSC [Elhamifar and Vidal, 2009], LRR 1 [Liu et al., 2013], LSR 2 [Lu et al., 2012] and CASS [Lu et al., 2013], the codes for which are downloaded from the authors' webpages. To obtain the best possible performance of the compared methods for different cases, we tune their corresponding parameters. Specifically, for k-NN, CAN and PCAN, the free parameter k is tuned from 1 to 10. For SSC, the space of the regularizer weight on R is α ∈ {2, 4, ..., 20}, λ ∈ {0.1, 0.2, ..., 5.0} for LSR, λ ∈ {0.0001, 0.001, 0.01, 0.1, 1.0, 2.0, 3.0} for LRR, while λ ∈ {0.0001, 0.001, 0.01, 0.1, 1.0} for CASS. To simplify our parameters, we let λ1 = λ2 = λ3 = ˆλ ∈ {0.1, 0.2, ..., 1.0}, although the simplification may very likely exclude the best performance for our method. By doing so, the parameter space is significantly shrunken, but contained by the original space. In other words, if the proposed algorithm with the shrunken parameter space outperforms the state-of-the-arts, the original parameter space can also achieve the same or possibly better performance. Due to space limit, we do not give the influence analysis of each parameter individually. Please note that our model also involves the parameter k, which will be fixed in the experiments according to a k effect testing (discussed in Sec. 4.2). Normalized Cuts [Shi and Malik, 2000] is employed to segment the input data into clusters for all the competitors, average segmentation accuracies 3 over 10 independent trials are finally reported.

1 As the authors proposed two versions of LRR with different models of E, i.e. the ℓ1 and ℓ2,1 , we test both of them and denote them as LRR1 and LRR21, respectively. For more details, please refer to [Liu et al., 2013].

2 LSR has two implementations, which are denoted as LSR1 and LSR2, respectively. Please refer to [Lu et al., 2012].

4.1 Synthetic Data

This part attempts to verify the robustness of our method to different levels of noise. We generate 5 independent subspaces {Ss}5s=1 of 4 dimensions, whose ambient dimension is 250. There are 100 unit data points randomly sampled from each subspace, a part of which are chosen to be corrupted with different levels of white Gaussian noise N(0, 1). For this experiment, the proportion of polluted data is fixed to 50%, the nearest neighbor number and ˆλ are set to 10 and 0.1 respectively, and the noise level varies from 0 (no noise) to 0.5. We evaluate the clustering performance of executing the Normalized Cuts on the affinity only (A, red curve), the presentation only (R, blue curve) and their Hadamard product (A⊙R, green curve). The combination of R and A is motivated by that the probability of two data points drawn from different subspaces simultaneously having high responses in R and A should be low. As can be seen from the first picture in Fig. 1, all the three schemes achieve very high accuracies when data are clean. But, as the noise level increases, the clustering performance of R only sharply drops to about 0.6. While the accuracies of A only (for most cases) and A⊙R are much higher. Please notice that, the red curve is always superior to the green until the noise level is up to 0.2. Afterwards, the green precedes the red. That is to say, using the affinity matrix only can provide a promising result on slightly polluted data, while further introducing the representation matrix is able to significantly boost the robustness to heavily corrupted data. The second picture in Fig. 1 displays the affinity matrix obtained by our method corresponding to the case of zero noise, which shows the perfect block sparsity. The rest two pictures in Fig. 1 give the affinity matrix and, the combination of affinity and representation with respect to the case with 0.25 noise, the block sparsities of which are not perfect but very well preserved.

3 The metric, segmentation accuracy, is calculated by finding the best matching between cluster labels and ground truth labels.

4.2 Extended Yale B

We compare the proposed method with other state-of-the-art methods for face clustering on the Extended Yale B dataset [Lee et al., 2005]. The dataset contains face images of 38 subjects. For each subject, there are about 64 frontal face images taken under different illuminations. More than half of the data vectors in this dataset have been corrupted by "shadows", which makes the task difficult. In this experiment, we resize the images into 32×32 and use the raw pixel values to form data vectors of 1024 dimensions.

k Effect. We use first 10 subjects to test the parameter effect of k, say the number of nearest neighbors. To eliminate the effect from other parameters, we empirically set ˆλ to 0.1. In addition, the three kinds of matrix including R only, A only and A⊙R are again employed to see the difference. As displayed in the left graph of Fig. 2 , it is easy to see that using A only and A⊙R give much more promising results when k ranges from 2 to 30 than using R only. Similar to the conclusion drawn from Sec. 4.1, A⊙R shows a better robustness than A only in this experiment. Based on this testing, we will fix k = 3 for our method for the rest experiments.

Convergence Speed. Without loss of generality, the convergence speed of Algorithm 1 by setting ˆλ = 0.1 on 10 subjects is given in the right picture of Fig. 2, in which the stop criterion sharply drops to the level of 10-6 with about 10 iterations and to 10-7 using 27 iterations. Our algorithm takes 4s to finish the computation on our PC, which is slower than LSR that spends 0.06s, but much more efficient than SSC (39s), LRR (60s) and CASS (34, 560s). This indicates that our proposed algorithm can converge sufficiently fast. Moreover, all the experiments conducted in this paper by our algorithm are converged with about 25 - 40 iterations.

Performance Comparison. Table 1 provides the quantitative comparison among the competitors on the Extended Yale B dataset. We evaluate the performance of the competitors on the tasks with different numbers of subject including 5, 10, 20, 30 and 38. The bold and italic numbers in each row represent the best and the second best results, respectively, for the corresponding task. Our parameter ˆλ is determined according to the highest accuracy of the case with 5 subjects, i.e. 0.1, and fixed for all the cases in this experiment. While, for different cases, each of the alternatives reports its best accuracies individually obtained by tuning its parameters in the corresponding parameter space. Please notice that the results of CASS with respect to the 5- and 10-subject cases are the best results reported by the authors of CASS [Lu et al., 2013]. For the rest cases, we do not provide the results as CASS takes too much time to handle even the 20-subject task (and not reported in [Lu et al., 2013]). As can be observed from Table 1, both our methods using A only and A ⊙ R with the uniform setting greatly outperform the others with tuned parameters for all the involved cases. We can also see that in this experiment A⊙R is slightly behind A only. The reason may be that although the data in this dataset are corrupted by different illuminations, they are well aligned and thus largely preserve the subspace structure. In [Feng et al., 2014], the authors state that their proposed scheme can improve the performance of LRR and LSR by 3%-6% on the (only reported) cases with 5 and 10 subjects 4 , even though, our method still significantly outperforms [Feng et al., 2014].

4.3 USPS

Further, we compare the performance of SSC, LRR, LSR, CASS and our method on the USPS dataset 5 , which consists of 10 classes corresponding to 10 handwritten digits, 0 ∼ 9. We use the first 100 examples with the size 16 × 16 of each subject for this experiment. The examples of each class are with many variations on appearance, and of different classes may share some features (e.g. digits 3 and 8), which violates the assumption of independent subspaces and thus increases the difficulty of clustering.

4 Since the code of [Feng et al., 2014] is not available when this paper is prepared, we do not explicitly compare with it.

5 www.csie.ntu.edu.tw/˜cjlin/libsvmtools/datasets/multiclass.html

Performance Comparison. As shown in Table 2, the best possible clustering accuracies of SSC, LRR, LSR and CASS are very close to each other, which fall into the range [72.20, 74.40]. These results are reasonably good although the USPS is more challenging than the Extended Yale B, as the subject number of the USPS is only 10 and the amount (sampling) of each subject of USPS is more than that of the Extended Yale B. As for our method, the performance of the scheme using A achieves 80.79, while that of the scheme using A ⊙ R obtains 82.58, which significantly improve the clustering accuracy on USPS compared with the others. It is worth mentioning that, in this experiment, the A⊙R scheme emerges its advantage over the A only.

4.4 UMIST

Moreover, we attempt to test the abilities of different approaches on a more challenging dataset UMIST [Graham and Allinson, 1998]. The UMIST collects 575 images from 20 subjects, which are resized into 56 × 46 for this experiment. Each subject has about 28 images with different poses, which significantly breaks the assumed subspace structure. Therefore, the performance of SSC, LRR, LSR, CASS and our method on subspace segmentation may degenerate or even fail, because they are primarily designed under the assumption of strong subspace structure.

Performance Comparison. Table 3 shows the performance comparison on the UMIST, from which we can see that the segmentation accuracies of LRR, LSR and CASS are around 0.51. This verifies the fact that the UMIST is a very challenging subspace segmentation dataset. Surprisingly, SSC and our method achieve reasonably high accuracies, i.e. 67.84, 68.09 and 70.12, respectively. This is mainly because SSC processes data individually instead of enforcing the subspace structure, which is more suitable for this dataset than the strategies of LRR and CASS. While our method introduces k nearest neighbor concept (with relatively small k) in the affinity matrix, which may connect the faces of the same subject with slight pose changes.

5 Conclusion

Subspace segmentation is an important yet challenging problem in many research fields, such as machine learning, computer vision and pattern recognition. Differently to previous spectral clustering based work that computes the affinity based either directly on the distance between data or on the similarity of data representations, this paper has shown a novel method that simultaneously learns data representations and their affinity matrix. We have formulated the problem into a unified optimization framework and designed an efficient Augmented Lagrangian Multiplier based algorithm to seek the solution. The experimental results, compared to the state-of-the-art alternatives, have demonstrated the clear advantages of the proposed method.

Appendix

Given a problem with the following shape:

where a ∈ Rn×1 is the target, and d ∈ Rn×1 is a known (distance) vector. To be more general, we can further appoint the number of nonzero elements in a as k ∈ {1, ..., n}. The closed form solution of the problem (12) is as follows:

where the elements of ˜ d ∈ Rn×1 are those of d but with the ascending order.

Proof. The Lagrangian function of the problem (12) is as:

where α ⪰ 0 and ω ⪰ 0 are the Lagrangian multipliers. It is easy to verify the optimal solution of a can be obtained through solving the following equation system, with the help of the KKT condition:

The third equation holds when the condition that if aj > 0 then ωj = 0 is satisfied, thus we have

And there are k positive elements in a ⪰ 0. Namely:

According to Eq. (16) together with aT 1 = 1, we have

Substituting α in Eq. (16) with 1+∑kj=1 ˜ dj k in Eq. (18) recognizes the form stated in Eq. (13). Please notice that only those data points with representation distances to the target smaller than 1+∑kj=1 ˜ dj k are its neighbors.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant 61402467 and in part by the Excellent Young Talent Program through the Institute of Information Engineering, Chinese Academy of Sciences.