Geometric Understanding for Unsupervised Subspace Learning

1 Introduction

Subspace learning is a kind of highly effective approaches to dimensionality reduction, which is widely applied in multiview data analysis [Ding and Fu, 2018], image classification [Fang et al., 2018], feature representation [Li et al., 2015], etc. In recent years, it becomes one of the mutual topics in computer vision, signal processing and statistical learning.

The goal of subspace learning is to map the highdimensional data to a low-dimensional space for representing more robust features while retaining as much information as possible. The basic idea is directly mining a latent lowdimensional manifold from the data, such that the local topological and/or geometric structures are preserved as much as possible. From the availability of training samples, subspace learning can be traditionally divided into three categories: supervised, unsupervised, and semi-supervised cases. Due to the fundamentality and simplification, in recent years, unsupervised subspace learning attracts more and more attention [Vaswani et al., 2018].

∗Contact Author

As a classical unsupervised approach, principal component analysis (PCA) finds the principal components, which are some orthogonal basis vectors most of data lies on, and all principal components span a subspace that well represents the data [Turk and Pentland, 1991]. Although PCA well concerns the globally linear structure of data, it does not catch its latently geometric structure. On the other hand, it is sensitive to the noise and outliers.

Considering the local geometric structure of data, Tenenbaum et al. [2000], Saul and Roweis [2003], and Belkin and Niyogi [2003] propose Isometric Mapping (ISOMAP), Locally Linear Embedding (LLE), and Laplacian Eigenmap (LE), respectively, by introducing the manifold assumption for data distribution. All these three approaches can effectively discover the locally geometric structure of the data by preserving data similarity before and after mapping, but they only provide the embedding results of training samples. Therefore, they cannot be extended to the classification problems. Moreover, they have high computational cost. For this, He et al. [2005b] establish the Locality Preserving Projection (LPP) method by approximating LE method, which well balances the manifold and the local structure of the data. Also considering the local structure, they propose a linear method, named Neighborhood Preserving Embedding (NPE), and apply it in classification problem [He et al., 2005a]. Later, Qiao et al. [2013] develop an explicit nonlinear mapping for manifold learning. More recently, Zhu et al. [2018] design an unsupervised spectral feature selection method by the selfexpressiveness of the features for preserving the local structure of features, and a low-rank constraint on the weight matrix to preserve the global structure.

On the other hand, to improve the robustness of PCA, in recent years, sparse representation (SR) [Wright et al., 2009] and low-rank representation (LRR) [Liu et al., 2013; Li and Fu, 2016] based methods are established. The SRbased methods can obtain a good robustness for the noisy data by seeking a best sparse representation for each group of samples, but they omit the global property of the data. For example, with the assumption of one subspace distribution, Cand`es et al. [2011] propose the Robust PCA (RPCA) method, which decomposes the data into low-rank background and sparse noise parts, and hence greatly promotes the robustness of data recovery. But in amount of real circumstances, data are not located in one subspace. Therefore, subspace clustering problem is proposed to describe the multi-subspace case. Elhamifar and Vidal [2013] discuss this task. Therein, Liu et al. [2010] adopt the nuclear norm for the wight matrix, which makes the acquirement of global representation easier. Further, to fit for data-missing case, Liu and Yan [2011] develop the Latent LRR (LatLRR) method. Later, Li and Fu [2016] develop a supervised regularization-based robust subspace (SRRS) method based on the label information and low-rank representation. More recently, Ding and Fu [2018] improve it to multi-view data analysis by collective low-rank subspace learning.

Although these methods improve more or less the robustness of subspace learning, they do not fully consider the geometric structure of the set of subspaces, which may further improve the performance of subspace learning. For example, Hamm and Lee [2008] propose a unifying view on the subspace learning method by formulating the problems as an optimization problem on the Grassmannian manifold, and performs feature extraction and classification in the same space. Li et al. [2008] develop an incremental subspace learning method by using the Log-Euclidean metric. Later, by using the non-Euclidean framework, Huang et al. [2014] consider the classification problem. More recently, Hauberg et al. [2014; 2016] find that the Karcher average of all one-dimensional subspaces spanned by normally distributed data coincides with the first principal component. Further, Chakraborty et al. [2017] extend the Karcher average computation for all one-dimensional subspaces to K-dimensional case. All these works validate the efficiency of Grassmannian manifold based methods, but there are still two shortages as follows. 1) All these methods need to have the dimension of subspaces as a prior knowledge, or need to enumerate all possible dimensions for subspace learning. 2) Grassmannian manifold is a graded manifold, hence it is difficult to represent the subspaces of different dimensions in the same framework.

Therefore, in this paper we will transfer the subspace learning on the Grassmannian manifold to the minimization problem on the rotation group, by representing all subspaces with different dimensions as an orbit of rotation group action on some standard subspaces. These standard subspaces are spanned by the orthonormal basis. Further, to avoid enumerating all possible dimensions for subspace learning, we introduce a low-rank regularizer which makes the dimension of subspace as low as possible.

The rest of this paper is organized as follows. In Section 2, we develop a new model for subspace learning from the rotation group action viewpoint. Then, an alternately iterative strategy and a structure-preserving algorithm are designed in Section 3. In Section 4, we demonstrate the effectiveness of our proposed method, and compare it with six state-of-the-art methods on two different kinds of real datasets (COIL-100, and MNIST). Finally, this paper is concluded in Section 5.

2 Geometric Model for Subspace Learning

We first consider the geometric structure of the rotation group SO(d). By the definition,

SO(d) = {R ∈ Rd×d |RRT = RT R = Id, det(R) = 1}, where Id is the identity matrix of the order d, and SO(d) is an m := d(d-1) 2 dimensional matrix Lie group with its Lie algebra so(d) defined by

which is the set of all skew-symmetric matrices. The Lie algebra so(d) is the linearization of Lie group SO(d) at the identity.

The goal of subspace learning is to find a latent lowdimensional manifold from the data for representing more robust features, while retaining as much information as possible. Therefore, it is natural to reduce the subspace learning to an optimization problem on the set of all subspaces. On the other hand, from the geometric viewpoint, all subspaces with dimension r of Euclidean spaces Rd form the Grassmannian manifold

where S is a subspace with the dimension r of Rd , and r is ergodic from 0 to d. Note that the trivial subspaces (r = 0 or r = d) should be eliminated in subspace learning, and hence we always assume that 1 ≤ r ≤ d - 1.

As mentioned above, the Grassmannian manifold Gr(d, r) is a graded manifold with respect to the dimension r, therefore, the first issue is how to uniformly parameterize this manifold for different r. After that, then, we can address the subspace learning on such manifold. Below, we will give a uniform representation for all subspaces first from the viewpoint of the rotation group action, and then we model the subspace learning as an optimization model on SO(d) ⋉ N.

2.1 Representation of Subspaces

Figure 1 gives an illustration for subspace learning. Given the set of centered samples {xi}ni=1 , we wish to find the best subspace (the hyperplane crossing the origin) such that it best describes the distribution of all samples.

To uniformly represent the subspaces, we adopt the approach of the rotation group SO(d) action. It is from the fact that there exists a rotation R ∈ SO(d), by which any subspace S with dimension r can be rotated from the standard subspace S0 spanned by r orthonormal basis. In detail, it is described as follows.

Let ei = (0, · · · , 1, · · · , 0)T be the ith orthonormal basis, Br = span{e1, e2, · · · , er} be the standard r-dimensional subspace, and Pr = [e1, · · · , er]T be the projection from Rd to Br. Then, for any r-dimensional subspace S ∈ Gr(d, r), we can find a rotation R ∈ SO(d), such that

That is,

Then, the unit normal vector Nr of the subspace S is the rotated unit normal vector of the standard subspace S0 by R. That is,

2.2 Model for Subspace Learning

By the representation of rotation group action, we reformulate the subspace learning as follows. As mentioned above, the goal of subspace learning is finding the best subspace such that it best describes the distribution of all samples. Therefore, first we should make the subspace best fitting the samples. That is, we should assume that the total distance of all samples from the subspace is as small as possible, where the distance between samples and the subspace can be formulated by

where R is the best rotation, Br is the standard subspace and N0r is the unit normal vector of Br.

Then, when the dimension of subspace is given by r, one part of model is to minimize the totally squared distance of all samples from the subspace with respect to the rotation R. That is,

On the other hand, due to the dimensional reduction, we want the dimension of the subspace is sufficiently small. Therefore, we introduce the regularization term to (1), and the model is reformulated by

where λ is a balanced parameter. Therefore, the subspace learning problem is transferred to the optimization problem on the SO(d) ⋉ N. It is remarkable that the regularity term is a low-rank constraint in previous literature because of rank(Pr) = r. On the other hand, from the Eq. (2), the first term ensures that the samples are as close as possible to the subspace. We know that this error decreases with the increment of r. Therefore, the model is well balanced the error and the dimension.

3 Algorithm

In this section, we will provide the solver for the model (2).

3.1 Iterative Strategy

It is seen that model (2) is an optimization problem with two kinds of independent variables. Therefore, we use an alternately iterative strategy. That is, we decompose the model (2) into two following subproblems.

(S1) For current rk , we update the next rotation Rk+1 by

(S2) For updated rotation Rk+1 , we update the next dimension rk+1 by

Then, by alternately updating Eqs. (3) and (4), we obtain the best rotation R∗ and the best dimension r∗ of the subspace, and moreover the best subspace S∗ = R∗Pr∗ .

Eq. (4) can be solved by traditional methods, since it is a finite and discrete minimization problem. Then, we will focus on the solution of the subproblem (S1) by the geometric structure of the rotation group SO(d).

3.2 Structure-Preserving Algorithm

To solve the subproblem (S1), we first consider the exponential map exp from so(d) to SO(d), by which we define the intrinsically iterative format as

where {Ej}1≤j≤m is the basis of Lie algebra so(d), and akj is the coefficient at step k. By this format, we insure that R at each step is an exact rotation, and hence it is a structurepreserving iterative method.

Then, we apply this method to solve the subspace (S1). Substituting the iterative format to Eq. (3), we have

Therefore, the subproblem (S1) is transferred to the conventional optimization problem. Then, to further simplify the route of solution, we adopt the quadratic approximation for the objective function by using the linearization for the exponential map [Helgason, 1978; Warner, 2013].

when aj is sufficiently small.

Let G(a) be this approximated function, then we have

where U(a) = m∑ j=1 ajEj , Y k = N0rkN T 0rk , and X is the data

matrix [x1, x2, · · · , xn].

Further, let Z k = RkT XXT Rk , then

where the last equation holds because

where M k is the matrix with element M k ij , and bk is the vector with element bki .

By the KKT condition, the optimal coefficients ak should have the equation ∇G(a) = 0 held. That is,

where [M k ]† is the Moore-Penrose pseudo inverse of M k . By this approximation, we deduce a closed form for the almost optimal coefficients. It highly improves the computational efficiency.

Therefore, we obtain the algorithm for subspace learning and summarize it as follows.

4 Experiment Results

To demonstrate the effectiveness of the proposed subspace learning algorithm, in this section, experiments of classification by using the k-NN classifier on learned subspace (reduced feature space) are conducted on two real datasets, i.e., COIL-100 object dataset [Nene et al., 1996], and MNIST digit dataset [Lecun et al., 1998]. Here, we more focus on the learned subspace, and hence classification results are Algorithm 1 Intrinsic algorithm for subspace learning

Require: Centered dataset X = {xi}ni=1 .

Ensure: The best subspace for data representation.

• 1: Initialize: The orthonormal basis {ei}, the basis of Lie algebra so(d), maximal iteration T, the balanced parameter λ, r0 = d and R0 = Id.
• 2: for each k ∈ [1, T] do
• 3: For current rk , update the next rotation Rk+1 by (5), where the coefficient vector ak is updated by (11).
• 4: For updated rotation Rk+1 , update the next dimension rk+1 by solving the finite discrete minimization problem (4).
• 5: end for
• 6: Output: The best subspace S∗ = R∗Br∗.

used to estimate the effects of the selected features. Therefore, for the classification task, we select the k-NN classifier which is more dependent with features. Further, we compare our proposed method with six state-of-the-art subspace learning methods, which include PCA (as a baseline experiment) [Turk and Pentland, 1991], LDA [Belhumeur et al., 1997], LPP [He et al., 2005b], NPE [He et al., 2005a], LSDA [Cai et al., 2007], and SRRS [Li and Fu, 2016]. All programs are written in Matlab 2013a and run by PC with Intel(R) Core(TM) i7-7500U CPU and 32 GB RAM.

4.1 COIL-100 Object Dataset

This dataset contains 100 objects. All images of one object are taken 5 degrees apart because the object is imaged on a rotated turntable, and hence each object has 72 images. The size of each image is 32 × 32 pixels with 256 grey levels per pixel. Thus, the original dimension of each image is 1024.

Then, we first randomly select 10 images of each object as the set of training samples, and the rest as the set of testing samples. Further, to make the results more repeatable, we do this operator 20 times randomly. Therefore, the results are obtained in the average sense. By executing six state-of-theart subspace learning methods and our proposed method, the numerical results are shown in Table 1, and Figures 2 and 3.

In Table 1, the first column is the data description, where the first letter P and B before the % sign represent pixel corruption and block corruption, respectively. The number in the parentheses means the dimension of the learned subspace. It is seen that the error rate of our proposed method is lower than other six state-of-the-art methods, while the dimension of the learned subspace is low. Specially, the performance is significantly improved in two cases of corruption, and the increment of performance increases with the growth of percentage of corruption. Therefore, our proposed method can obtain better feature representation as well as robustness. This is also reflected in Figure 2. Also, we provide the visualization results with all methods by t-SNE [Maaten and Hinton, 2008] in Figure 3, by which justifies that our proposed method obtains the best result.

4.2 MNIST Digit Dataset

The MNIST dataset of handwritten digits has a training set of 60,000 examples, and a test set of 10,000 examples. The digits have been size-normalized and centered in a fixed-size image. The size of each image is 28 × 28 pixels. Thus, the original dimension of each image is 784. We add 10% pixel corruption, which is the same as the case of COIL-100 dataset.

Then, we first randomly select 1000 images with 10 numbers as our dataset, such that there are 100 images in each group. Then, to test the error rate of classification in the different number of training samples, we randomly select 10, 20, and 30 images of each class as the set of training samples, respectively. Further, to make the results more repeatable, we do this operator 20 times randomly. Therefore, the results are obtained in the average sense. By executing six state-ofthe-art subspace learning methods and our proposed method, the statistical results are shown in Figure 4. From Figure 4, we see that the error rate of classification by our proposed method are almost lowest in the different number of training samples.

4.3 Discussion

For the computational time, we compare our proposed method with SRRS on COIL-100 object dataset and MNIST digit dataset without corruption. Further, to make the results more repeatable, we conduct this operator 20 times randomly. Therefore, the training time results in the average sense are shown in Table 2, which shows that our proposed method obtains the best classification result with less training time than SRRS.

5 Conclusion

In this paper, we propose an unsupervised subspace learning from a geometric viewpoint. First, we represent the set of all subspaces by the orbit of the rotation group action on the standard subspace. Then, the subspace learning is transferred to a minimization problem on the rotation group. It provides a well geometric explanation of subspace learning. Further, we construct an intrinsic algorithm by applying such geometric structure. Finally, we compare the proposed approach with six state-of-the-art methods on two different kinds of real datasets. The experimental results validate that our proposed method outperforms all compared methods.

Acknowledgments

The research is supported by the National Natural Science Foundation of China under Grant Nos. 11771276, 61573274, 61731009, and the capacity construction project of local universities in Shanghai (18010500600).