Cascaded Low Rank and Sparse Representation on Grassmann Manifolds

1 Introduction

Unsupervised clustering is a fundamental topic in computer vision and artificial intelligence areas, which aims to group data with different patterns into different clusters by exploring the intrinsic membership of data. Recently, benefiting from the efficiency of analyzing complex relationship and structures hidden in data, spectral clustering method is considered having better prospects. Generally, in spectral clustering method, an affinity matrix is firstly learned from the given data, and then the affinity matrix is fed to a standard clustering algorithm such as NCut or K-means to obtain the final clustering results. The data representation and its affinity matrix have heavy impacts on the clustering results, so numerous efforts have been made to exploit better affinity matrices to reveal the intrinsic structure embedded in data. Sparse Subspace Clustering (SSC) [Elhamifar and Vidal, 2011] and Low Rank Representation (LRR) [Liu et al., 2013] are two representative methods. They respectively use sparse and low rank constraints on the affinity matrix of self-expression model.

∗Corresponding author: Yongli Hu (huyongli@bjut.edu.cn)

In the view of subspace, SSC is regarded to reveal the "local" structure of the data, as the sparsity of a datum representation is resulted by the hypothesis that only the data representation coefficients corresponding to the data samples belonging to the same subspace are nonzero. LRR tries to formulate the "global" structure of data space via low rank constraint which is a holistic constraint to reflect the relationship among data. In the subspace clustering, the ideal affinity matrix need have diagonal block shape and each block can represent a subspace [Feng et al., 2014; Wu et al., 2016]. According to the conclusion from [Luxburg, 2007], the number of blocks is determined by the rank of the Laplacian matrix which is constructed by the affinity matrix; thus, LRR essentially describes the relationship among subspaces. Due to the above advantages, SSC and LRR and their varieties achieve state-of-the-art clustering performance.

With the success of SSC and LRR mentioned above, it is intuitive to combine the global and local structure constraints together to obtain further improvements. Several such combined models have been proposed in recent researches, such as Low Rank Sparse Subspace Clustering [Zhuang et al., 2012; Patel et al., 2015] and Laplacian regularized LRR [Yin et al., 2016]. These combined models which simultaneously enforce sparse and low rank constraints on the same affinity matrix have achieved some gains, but this combining strategy is not considered the best way to reveal the global and local structures of data. It is just a trade-off between the two penalties used in the current methods in general. However, the sparse and low rank properties of a matrix usually cannot be satisfied simultaneously, for example, the full rank matrix of identity matrix or diagonal matrix. From this observation, to fully reveal the global and local structures of the data, we propose a cascaded version of low rank and sparse subspace clustering methods, namely Cascaded Low Rank and Sparse Representation (CLRSR). In the proposed method, the low rank and sparse properties of data are successively mined by a global to local manner in two cascaded steps: First, the global structure of data is represented by LRR; Then the sparse representation of the former learned low rank representation is exploited to capture the local structure of data.

Another observation is that most of the current subspace clustering methods take the data as vectors with Euclidean distance as the metric. However, for many imageset, video data or other multi-dimension data, the simple vectorized data form and its metric are not applicable since it has been proven that many high-dimension data are embedded in lowdimension manifolds with nonlinear metric [Wang et al., 2008]. So it is rational to model such multi-dimension data as manifolds. For this purpose, we adopt the Grassmann manifolds to represent multi-dimension data for clustering, because Grassmann manifolds is a good tool for representing imageset or video data [Turaga et al., 2008] and has a nice property that it can be embedded into the linear space of symmetric matrices. To this end, we further conquer the challenging problem of constructing our CLRSR on the manifold space, which results in a novel data clustering on Grassmann manifolds, called G-CLRSR. We illustrate the main idea of our proposed method in Figure 1 to compare with the conventional sparse and low rank subspace clustering methods.

The contributions of this work are summarized as follows:

• Proposing a novel cascaded low rank and sparse representation model to fully capture the global and local structures underling the observed data;
• Extending the CLRSR model onto Grassmann manifolds which favors the representation of multi-dimension data; and
• Giving a practical solution to the optimization problem on manifold space for the proposed G-CLRSR model.

2 Relevant Concepts

Throughout the paper, all matrices/vectors are written as bold uppercase/lowercase, i.e., X/x. A subscript is used to label the sequence of matrices/vectors, i.e., Xi/xi. Transpose matrix XT and inverse matrix X-1 are also defined. The calligraphic letter X and E represent tensor data and the italic letters (V,N, v) denote scalar values or integers. Other special symbols will be explained when they are used.

2.1 LowRank and Sparse Representation (LRSR)

To explore the hidden global structure, one often enforces the low rank constraint on the self-expression of data. The local neighbor relationship is also critical to subspace clustering; therefore, one joins the neighbor relationship and the global structure to improve the clustering performance [Zhuang et al., 2012; Patel et al., 2015] as,

where Y = [y1, y2, · · · , yN] ∈ Rd×N denotes a set of samples collected from multiple independent subspaces. rank(·) computes the rank of matrix and ℓ0-norm ‖ · ‖0 counts the number of non-zero elements of matrix. The reconstructed error is measured by Euclidean distance.

Different from above LRSR model, our proposed method exploits the sparse representation of the former learned global structure to capture the neighborhood structure. Additionally, to handle multi-dimension data, we further extend it from Euclidean space onto manifold space. Here the Grassmann manifold is adopted.

2.2 GrassmannManifolds

Definition 1 (Grassmann Manifolds) [Absil et al., 2008] The Grassmann manifold, denoted by G(p, n), consists of all the p-dimensional subspaces in n-dimensional Euclidean space Rn (0 ≤ p ≤ n).

To intuitively represent Grassmann manifolds, the popular way is to construct Grassmann manifolds by

From this, we can see that a Grassmannian point actually is an equivalent class (a subset) in which any two thin-tall orthogonal can be converted to each other applying a p×p orthogonal matrix.

Definition 2 (Embedding Distance) [Harandi et al., 2013] Grassmann manifold can be embedded into the symmetric matrices space by

Given two Grassmannian points X1,X2 ∈ G(p, n), the embedding induces the following distance on Grassmann manifolds, defined as,

There exist different types of geometries on Grassmann manifold, please refer to [Absil et al., 2008].

3 Cascaded LowRank and Sparse Representation on Grassmann Manifolds

In this Section, we will elaborate a cascaded version of low rank and sparse representation model for clustering.

3.1 CLRSR on GrassmannManifolds

To explore the global structure, we enforce the low rank constraint on the self-expression coefficients matrix of data. Then to fulfill the local properties, the low rank representation of data is cascaded with the sparse representation. So we have

where λ1, λ2, and λ3 are three major balancing parameters, and the coefficient matrices Z,C ∈ RN×N . This model is named as CLRSR model. As a common precessing, we replace the low rank and sparse items by the nuclear norm and ℓ1-norm as the effective approximations,

To handle multi-dimension data, we consider the generalization of model (6) onto Grassmann manifolds. Given a set of Grassmannian points X = {X1,X2, · · · ,XN} where Xi ∈ G(p, d), the parallelized model is given by

where we use abstract symbols ⊖, ⊎Nj=1 and ⊙ to simulate the "linear" operations to be defined on manifolds, i.e., addition, subtraction and scalar-multiplication. ‖Xi⊖(⊎Nj=1 zij⊙ Xj)‖G measures the error between the Grassmannian point Xi and its "reconstruction" ⊎Nj=1 zij ⊙Xj . According to the property of Grassmann manifolds in (3), it is easy to use the embedding distance to replace the manifold distance as,

where X = {X1XT1 ,X2XT2 , · · · ,XNXTN } is a 3-order tensor and ×3 is the mode-3 tensor-vector multiplication [Kolda and Bader, 2009]. This error measurement not only avoids using Log map operator but also has simple computation with F-norm. Finally, the CLRSR model on Grassmann manifolds can be defined as follows,

where the reconstructed error E is also a 3-order tensor. This model is named as G-CLRSR.

3.2 Optimization

Following the notation used in [Wang et al., 2014], we firstly simplify the tensorial error term E in formula (9) to avoid a complex computation between tensor and matrix,

where E(:, :, i) = XiXTi - ∑Nj=1 zij(XjXTj ) is the ith front slice of E. For convenience, we denote ∆ij = tr[(XTj Xi)(XTi Xj)] and the N × N symmetric matrix ∆ = [∆ij]Ni=1,j=1 . After a sequence of simple mathematical manipulations, we can show that problem (9) is equivalent to

We use Alternating Direction Method [Lin et al., 2011] to solve such a convex problem. For this purpose, we introduce two auxiliary variables J = Z and Q = C to separate variables Z and Q, and reformulate the optimization problem as

which can be solved by minimizing its ALM formula

where A and B are the Lagrange multipliers, and µ is the penalty parameter. The above ALM is appealing on if we can find closed form solutions to the following subproblems (13), (14), (15) and (16).

Updating J. Fixing other variables, we update J by,

A closed-form solution is suggested as the following theorem [Cai et al., 2010],

Lemma 1 Given that UDVT = SVD(Z + Aµ ) as defined above, the solution is given by

where D is the diagonal matrix of singular values. The singular value thresholding operator is defined as Sτ(x) = sign(x) · max(|x| - τ, 0) where τ is a positive number.

Updating Z. To update Z, we solve the following problem by fixing other variables,

This is a quadratic optimization problem with respect to Z. We set its derivation as zero and get the closed solution as,

Updating Q. By fixing other variables, we update Q through,

We know that the problem in (15) has a closed-form solution, given by [Beck and Teboulle, 2009],

Lemma 2 Let matrix M = C+ Bµ and a positive number τ , a closed-form global optimal solution of above problem is

where the singular value thresholding operator is also defined as Sτ = sign(x) · max(|x| - τ, 0).

Updating C. The final self-expression C is updated by solving the following problem,

Similar to update the subproblem Z, we also set the derivation of (16) to zero and get the closed-form solution,

Updating A and B. After we solving the above subproblems on J,Z,Q and C, respectively, we can easily update these two linear Lagrangian multipliers by

Adapting penalty parameter µ. For the penalty parameter µ > 0, we could update it by

where µmax is the pre-defined upper bound of µk .

Once obtaining the coefficient matrix C, any clustering algorithm, i.e., NCut used in this paper, can be performed on affinity matrix W = (|C| + |C|T )/2 to receive final clustering results. We summarize pseudo code in Algorithm 1.

Algorithm 1 The G-CLRSR Optimizing Problem.

Input: Grassmannian points {Xi}Ni=1 , Xi ∈ G(p, d), the balancing penalty parameters λ1, λ2 and λ3. Output: The coefficient matrix C.

• 1: Initialize: J = Z = N = Q = 0, A = B = 0, ρ = 1.1, ϵ = 10-8 , µ0 = 10-6 and µmax = 106 .
• 2: ##Calculating ∆
• 3: for i=1:N do
• 4: for j=1:N do
• 5: ∆ij ← tr[(XjXTj )(XTi Xi)]
• 6: end for
• 7: end for
• 8: ##Optimizing variables
• 9: while not converged do
• 10: Updating variable J,Z,Q,Cvia (13-16), respectively.
• 11: Updating multipliers A and B via (17).
• 12: Updating parameter µ by µ =min(ρµ, µmax).
• 13: Checking the convergence conditions:

• 14: end while

3.3 Complexity and Convergence Analysis

As shown in Algorithm 1, the major time-consuming components lie in steps 5, 10, 11, and 13. (i) The complexity of step 5, the computation of the symmetric matrix ∆, is O(p2d + p3 ) as its main complexity comes from 4 matrices multiplication; (ii) Step 10 takes a SVD operation for updating J, leading to a complexity of O(N 3 ); and (iii) Both steps 11 and 13 involve matrix multiplication and inversion operations, thus the complexity is O(N 3 ). Overall, the total complexity is O(N 2 (p2d + p3 ) +N 3 ) for each iteration. Under the condition p ≪ d, the total complexity is basically O(N 2 p2d + kN 3 ) where k denotes the number of iterations.

To experimentally show the convergence behavior of the proposed method, we provide an intuitive curve illustration of the convergence with respect to the iteration number on all datasets, shown in Figure 2. Our proposed algorithm converges in about 230 steps, which reflects that our optimization method has a good convergence property.

4 Experiments

4.1 Datasets

Extended Yale B dataset is collected from 38 subjects and each subject has 64 front face images. We randomly choose 8 images from the same subject to construct an imageset.

CMUPIE 1 dataset contains totally 1632 front face images evenly from 68 subjects. Every 4 images from the same subject are chosen to construct an imageset.

UCF sport dataset 2 includes a total of 150 video clips with 13 actions, which has a wide range of scenes and viewpoints. Each clip has 22 to 144 frames.

SKIG gesture dataset 3 contains 1080 RGB-D video clips with 10 gesture types of 6 subjects. All the gestures are performed under different backgrounds and illumination. Each sequence contains 63 to 605 frames.

Highway traffic dataset 4 contains 253 video clips labeled with three levels. Each clip has 42 to 52 frames.

Each video clip is regarded as an imageset to construct a point on Grassmann manifolds.

4.2 Comparison Algorithms and Settings

For the comparison algorithms, we have the following (i) vector-form input methods, LRR [Liu et al., 2013], SSC [Vidal, 2011], LS3C [Patel et al., 2013] and CLRSR, (ii) Grassmann-form input methods, SCGSM [Turaga et al., 2011], SLRR [Yin et al., 2015], FGLRR [Wang et al., 2014], LapFGLRR [Wang et al., 2018] and a baseline G-LRSR.

To execute LRR, SSC, LS3C and LRSR on imageset and video data, we have to vectorize each imageset as inputs for them. However, in most experiments, we cannot simply take these long vectors into the algorithms due to the high dimension for a larger imageset. So, we reduce the dimension of these vectors to a low one which equals to the dimension of PCA components retaining 95% of its variance energy.

SLRR extends LRR onto Stiefel manifolds which has the same form with Grassmann manifolds. SCGSM and FGLRR can be considered as the Grassmann manifolds version of Kmeans and LRR models. In LapFGLRR, a local constraint among data, Laplacian regularization, is introduced into the FGLRR model. In G-LRSR, the LRSR model is embedded into Grassmann manifolds as a important baseline.

For a fair comparison, we execute the code of compared methods provided by authors. The parameters λ1, λ2 and λ3 in CLRSR and G-CLRSR, one parameter λ1 in LRR, SSC, FGLRR and SLRR, two parameters λ1 and λ2 in LapFGLRR, G-LRSR and LS3C, need to be tuned, respectively. We tune their values within the set of {0.1, 0.2, · · · , 1, 2, · · · 10}, and we report the best value in each experiment.

To make a comprehensive evaluation, we employ five clustering validation measurements, including ACC (accuracy), NMI (normalized mutual information), RI (rand index), Purity and FM (F-measure). For all the metrics, higher values indicate better performance.

1 http://www.cs.cmu.edu/afs/cs/project/PIE/MultiPie/Multi-

Pie/Home.html

2 http://crcv.ucf.edu/data/

3 http://lshao.staff.shef.ac.uk/data/SheffieldKinectGesture.htm

4 http://www.svcl.ucsd.edu/projects/traffic/

4.3 Data Representation

An imageset F = {f1, · · · , fP}; fi ∈ Rd , with fi being the feature descriptors of i-th frame. In this paper, we just use the image intensity as the raw feature. The Grassmannian point can be represented by an orthogonal basis. Similar to the work in [Wang et al., 2014], the feature matrix F ∈ Rd×P is decomposed by SVD as F = UΣVT ; then we pick up the first p (p ≤ P) singular-vectors of U to represent a point X = [u1, u2, ..., up] on Grassmann manifolds G(p, d).

4.4 Experimental Results

Image Set Data Clustering

Although facial recognition and clustering has been a mature application in computer vision, it often gets unfavorable experimental results since the varieties of facial expressions, poses or light conditions. Imageset data collect complementary information from multiple views of each object, which can naturally overcome the environment affection. Here, we regard each imageset as one point on Grassmann manifolds.

Table 1 reports all clustering results in terms of different validation criteria on two facial imageset datasets. From an overall perspective, manifold representation based methods obtain better performance than vector-form input clustering methods (i.e., LRR, SSC and LS3C), which demonstrates that manifold representation is a good tool for imageset data. Compared with SLRR, SCGSM and FGLRR, our proposed method simultaneously considers meaningful global and neighbor structures, which makes our proposed method more robust in clustering tasks. As for the important baselines, i.e., LapFGLRR and G-LRSR, they both involve global and local constraints on only one layer of representation, while our proposed approach can extract meaningful representation layer by layer. Through keeping the local neighbor relation on the base of global structure information, our proposed method receives the best experimental results.

Video Data Clustering

Action recognition is another open and hot issue in computer vision literature. Most action data are in the form of video clips, which is a set of sequential frame images and suitable to our manifold representation methods. So we select a human action and a gesture video datasets.

Table 2 lists the results of different clustering methods on video datasets. These datasets are much more challenging than the previous two fical imageset datasets, as the environment is more complex and varies dramatically. Obviously, for video datasets, manifold representation based methods reflect more impressive superior to vector-input based methods, which is further to verify the advantage of manifolds representation. What is more, LapFGLRR, G-LRSR and our proposed methods report the superior results in all five metrics, which also verify that fusing global and local information are important for video data on clustering tasks. It should be noted that the performance of most compared methods are not robust to all datasets. For example, LapFGLRR achieves the best performance on SKIG in terms of RI . However, the performance on the other datasets are not such promising.

In addition, we wish to evaluate our proposed method in a real and practical application; therefore, a natural scenario dataset, Highway Traffic videos, is chosen as test data. Table 3 shows the clustering results of all the algorithms, our proposed method still get satisfied performance, i.e., almost 4 percent rasing performance in ACC. For the general application of G-CLRSR, we can use it to learn the threshold of different level or part the period for the traffic jam on specific roads based on the historical traffic data.

5 Conclusion

In this paper, we proposed a novel cascaded low rank and sparse representation clustering method and extended it onto Grassmann manifolds to handle multi-dimension or imageset data. Compared with existing works which integrate low rank and sparse constraints on the same data representation simultaneously, the proposed methods execute a sparse constraint on the low-rank representation matrix in a cascaded manner which contributes to obtaining an improved affinity matrix. In addition, an efficient alternative based algorithm is proposed to solve for the optimal solution. Intensive experiments were conducted on imagesets and video datasets to demonstrate the superior performance of our methods.

Acknowledgements

The research project is supported by the Australian Research Council (ARC) through the grant DP140102270 and also partially supported by National Natural Science Foundation of China under Grant No. 61390510, 61672071, 61632006, 61772048, Beijing Natural Science Foundation No. 4172003, 4162010, 4152009, 4184082, Beijing Municipal Science and Technology Project No. Z171100000517003, Z171100000517004, Z171100004417023, Z161100001116072.