Linear Manifold Regularization with Adaptive Graph for Semi-supervised Dimensionality Reduction

1 Introduction

Dimensionality reduction is a significant topic in machine learning and other related fields. It is reasonable to presume that the naturally generated high dimensional data have a much more compact description, i.e., the high dimensional data probably lie on or close to a smooth low dimensional manifold [Roweis and Saul, 2000]. The goal of dimensionality reduction is to remove the noise and redundant information, and at the same time to preserve the desired intrinsic information of the input data. Collecting the labeled data is usually costly, while the unlabeled data are abundant and can be easily obtained. Therefore, semi-supervised dimensionality reduction has attracted great interest in recent years.

If we do not have more information than similarities between data points, a nice way to represent the data is in the form of a graph [Zhang et al., 2014; Liu et al., 2010], which aims to capture the intrinsic geometric structure of data manifold. There have been many graph-based methods for dimensionality reduction. To provide a unified perspective of various algorithms, [Yan et al., 2007] proposed a general framework known as graph embedding, in which the algorithms such as LLE [Roweis and Saul, 2000], LE [Belkin and Niyogi, 2001] and LPP [He et al., 2005] share the common formulation with different graph design. To better cope with the data sampled from nonlinear manifold, [Nie et al., 2010] proposed the flexible manifold embedding (FME) framework for semi-supervised and unsupervised dimensionality reduction. There are many other semisupervised graph-based methods that were developed with different prior assumptions [He et al., 2008; Gao et al., 2015; Chatpatanasiri and Kijsirikul, 2010] or by label propagation [Nie et al., 2009]. By adding a graph regularization term, some supervised dimensionality reduction methods can also be extended to the semi-supervised case [Cai et al., 2007; Song et al., 2008; Huang et al., 2012].

All the graph-based methods mentioned above need to construct a graph beforehand. Therefore, graph construction is a crucial step for these methods, since their performance highly relies on how well the graph models the intrinsic structure of data manifold. In general, one can construct an adjacency graph by k-nearest neighbor or ε-ball neighborhood criteria. The edge weights are then assigned by Gaussian kernel or local linear reconstruction. However, due to the noise and redundant information, such a pre-defined graph has no clear structure and may not be appropriate for the subsequent dimensionality reduction task.

To overcome the drawbacks, it is natural for us to consider how to learn an adaptive graph which is the optimal one for dimensionality reduction. The adaptive graph should better be sparse and have clear structure that the number of connected components in the graph is exactly the number of data clusters/classes. Such a structured graph would be benifical to many tasks since it contains more accurate information of the data. Motivated by these ideas, we propose a novel approach called linear manifold regularization with adaptive graph (LMRAG) for semi-supervised dimensionality reduction. It is worthwhile to highlight the main contributions of the paper as follows:

1. LMRAG performs dimensionality reduction and graph construction simultaneously, by incorporating the adaptive neighbor learning into the objective function of linear Laplacian regularized least squares (LapRLS/L). Both the optimal graph and the projection matrix can then be obtained.
2. To learn an adaptive graph that has clear structure, a structure constraint is imposed to the graph Laplacian. To the best of our knowledge, it is the first time to introduce an adaptive and structured graph for semi-supervised

dimensionality reduction.

1. A simple yet effective algorithm is developed for our new model. Extensive experiments on several widely used datasets demonstrate the effectiveness of the proposed method.

2 Background

We first introduce some notations used throughout the paper. For a matrix W ∈ Rm×n , the (i, j)-th entry and the i-th column are denoted by wij and wi, respectively. The trace and Frobunius norm of W are denoted by T r(W) and ‖W‖F , respectively. The p-norm of vector v is denoted by ‖v‖p, and Ik ∈ Rk×k is an identity matrix. 1 ∈ Rn×1 is a vector with all entries being 1.

The data matrix is denoted by X ∈ Rd×n (n = l+u), where the first l samples {xi}li=1 are labeled and the last u samples {xi}ni=l+1 are unlabeled. c is the number of data classes. The label matrix Y ∈ Rc×n is defined as yji = 1 if xi has label j ∈ {1, 2, . . . , c} and yji = 0, otherwise. Let G = {X, S} be an undirected and weighted graph, in which X is viewed as the vertex set and S ∈ Rn×n is the similarity matrix. The entry sij measures the similarity between xi and xj. The graph Laplacian is then defined as L = D - S, where the diagonal matrix D has entry dii = ∑nj=1 sij (i = 1, . . . , n).

2.1 Linear Manifold Regularization

Manifold regularization [Belkin et al., 2006; Sindhwani et al., 2005a; 2005b] is a widely used geometric framework that brings together three distinct concepts from the theory of regularization in reproducing kernel Hilbert spaces (RKHS), manifold learning and spectral methods. It has successfully extended linear regression and support vector machine (SVM), respectively, to the semi-supervised learning methods Laplacian regularized least squares (LapRLS) and Laplacian SVM. We take LapRLS/L as an example to briefly introduce the linear manifold regularization. The formulation of LapRLS/L is as follows:

where W ∈ Rd×c is the projection matrix and b ∈ Rc×1 is the bias term. The third term is the label fitness term. γA, γI are two regularization parameters that control the RKHS norm and the intrinsic norm, respectively.

2.2 Adaptive Neighbor Learning

We consider the probabilistic neighbors to learn the similarity matrix. The probability of two data points to be neighbor can be regarded as their similarity [Nie et al., 2014]. It is natural to presume that a smaller distance should be assigned a larger probability, and vice versa. For simplicity, we adopt the Euclidean distance. Therefore, we can adaptively determine the probabilities by solving the following problem:

where γ > 0 is the regularization parameter, and si ∈ Rn×1 is a vector with the j-th entry as sij . The regularization term s2ij is used to avoid the trivial solution that the nearest neighbor has the probability of 1 while the others are all 0. We do not consider xi itself is the neighbor of xi. For the diagonal entries of S, we simply have {sii = 0}ni=1 . In Eq.(2), we can measure the distance in the projected space by replacing xi, xj with W T xi, W T xj , respectively. Moreover, we enforce S to be symmetric by (ST + S)/2.

3 The ProposedMethod

3.1 Formulation

It is not difficult to verify that

Based on Eq.(1) to Eq.(3), by incorporating the adaptive neighbor learning into the objective function of LapRLS/L, the proposed LMRAG is formulated as follows:

where α, β and γ are three trade-off parameters. The fourth term is the label fitness term. U is a diagonal matrix with the first l and the last u diagonal entries being 1 and 0, respectively.

By solving Eq.(4), we can learn an adaptive graph while in most cases all the data points are in one connected component. According to [Mohar et al., 1991], the multiplicity c of eigenvalue 0 of the graph Laplacian matrix is equal to the number of connected components in the graph. Therefore, to make the adaptive graph structured, we can add a structure constraint by restricting the rank of L to be (n - c).

However, it is challenge to directly solve the problem of Eq.(4) with the rank constraint. Suppose σi(L) is the i-th smallest eigenvalue of L, we can transform the rank constraint to the sum of the first c smallest eigenvalues. Note that σi(L) ≥ 0, since L is positive semi-definite. The objective function of LMRAG then becomes:

As we can see, for a large enough λ, solving Eq.(5) will make ∑ci=1 σi(L) get infinitely close to zero, then the rank constraint is approximately satisfied. Such a relaxation is beneficial to the subsequent optimization, while Eq.(4) with the rank constraint is hard to tackle.

Further, we have the following equation:

where the optimal F is formed by eigenvectors of L corresponding to the first c smallest eigenvalues as rows. fi ∈ Rc×1 can be seen as a kind of embedding of xi. The term on the right hand side of Eq.(6) is actually the objective function of spectral clustering [Von Luxburg, 2007]. Therefore, our final objective function is formulated as follows:

3.2 Optimization

We divide the problem in Eq.(7) into three subproblems, and propose an alternative and iterative algorithm to optimize them. The whole procedure is summarized in Algorithm 1.

Step 1: Update F with W, b and S fixed. The problem in Eq.(7) becomes:

The optimal solution F is formed by eigenvectors of L corresponding to the first c smallest eigenvalues.

Step 2: Update W, b with F and S fixed. The problem in Eq.(7) becomes:

To obtain the optimal solution, by setting the derivatives of the objective function with respect to W and b equal to zero, respectively, we have:

where Hcu = U - 1l U11T U is the centering matrix for the labeled data.

Step 3: Update S with W, b and F fixed. The problem in Eq.(7) becomes:

We have 2Tr(FLF T ) = ∑ni,j=1 ‖fi - fj‖22 sij which is similar to Eq.(3). Note that it is independent to conduct adaptive neighbor learning for each data point. Thus we can solve the following problem for the i-th sample:

Algorithm 1 The Proposed Method LMRAG

Input: Data matrix X ∈ Rd×n , where {xi}li=1 are labeled and {xi}ni=l+1 are unlabeled, label matrix Y ∈ Rc×n , trade-off parameters α,β, and the neighbor number k.

• 1: Initialize S, λ, γ according to the initialization section.
• 2: while not converge do
• 3: Update F, which is formed by eigenvectors of L corresponding to the first c smallest eigenvalues.
• 4: Update W,b by Eq.(10).
• 5: Update S by solving Eq.(13) for each sample.
• 6: end while

Output: Projection matrixW.

Denote dij = ‖W T xi -W T xj‖22 + λ‖fi - fj‖22 , and denote di ∈ Rn×1 as a constant vector with the j-th entry as dij, Eq.(12) can be rewritten as follows:

The problem in Eq.(13) naturally has a sparse solution and can be solved by an efficient iterative algorithm [Huang et al., 2015]. We can also just update the k nearest similarities for each sample to ensure a sparse solution.

3.3 Initialization

We can learn an initial graph by solving the problem of Eq.(2), and the algorithm proposed in [Huang et al., 2015] can be adopted again. Alternatively, based on k-nearest neighbor (KNN) assumption, we apply another strategy to tackle the problem, and at the same time to determine the parameter γ. The Lagrangian function of Eq.(2) for the i-th sample can be written as follows:

where zij = ‖xi - xj‖22 , η and ξ ∈ Rn×1 are the Lagrangian multipliers. zi ∈ Rn×1 is a constant vector with the j-th entry as zij, and the overall γ can be set to the average of {γi}ni=1 . Based on the KKT condition, the optimal si has sij = (- zij 2γi + η)+, where (z)+ = max(z, 0).

We consider that each sample has k nearest neighbors, i.e., si has k nonzero entries. Let us rank zi in ascending order, we have

In above derivations, we get a value range for γi and we set it to the maximum. Consequently, the initial S can be computed by

As for parameter λ, in practice, we can use a dynamic strategy to determine it and to accelerate the iterative process. Specifically, we can initialize λ = γ. Denote the number of connected components in S as ncc, then in each iteration, we double λ if c > ncc, halve λ if c < ncc, and we stop the iteration, otherwise.

3.4 Computational Complexity

The complexity of step 1 is O(n3 ). Considering that L is sparse, the ARPACK eigensolver [Lehoucq et al., 1998] can be adopted to reduce the cost to (O(p3 )+[O(np)+O(nk)]× O(p - c)) × T , where p is a value several times larger than c, and T is the number of times of restarted Arnoldi.

Step 2 mainly takes O(n2d + nd2 + d3 ). We can use the Nystr¨om method [Fowlkes et al., 2004] to reduce the cost of the inverse operation performed on a symmetric matrix. Woodbury formula can also be used when d > n.

Comparing to step 1 and step 2, the complexity of step 3 can be ignored since the algorithm proposed in [Huang et al., 2015] is based on Newton method that has quadratic convergence rate, and in practice we can just update the local similarities.

4 Discussions

There have been extensive study on the problem of dimensionality reduction. Besides the traditional KNN graph, there also exist many graph construction methods [Liu et al., 2010; Zhang et al., 2014]. However, most graph-based methods conduct graph construction and dimensionality reduction in two separate steps, and a very limited number of works have devoted to learning an optimized graph for dimensionality reduction.

Graph optimized locality preserving projection (GoLPP) [Zhang et al., 2010] is the first attempt to perform graph optimization during a specific dimensionality reduction task according to the authors. The idea of GoLPP is to regularize the objective function of asymmetrical LPP [He et al., 2005] by an entropy term. However, GoLPP suffers the nonuniqueness of the solutions, since GoLPP is formulated in the trace ratio form while solved in the ratio trace form [Wang et al., 2007; Jia et al., 2009]. Due to the entropy regularizer, the graph learned by GoLPP is dense even though a sparse initial graph is given.

To address the problems of GoLPP, graph optimization for dimensionality reduction with sparsity constraints (GODRSC) was then proposed [Zhang et al., 2012], based on the orthogonalization of sparsity preserving projections (SPP) [Qiao et al., 2010]. GODRSC obtains the sparsity of graph by replacing the entropy regularizer in GoLPP with an ℓ1-norm minimization, and avoids nonunique solution by directly solving the trace ratio formulation.

GoLPP and GODRSC are both proposed for unsupervised dimensionality reduction. Therefore, LMRAG is of great value as an effective extension of the existing graph optimized dimensionality reduction methods in semi-supervised case. In fact, LMRAG can be easily extended to the unsupervised case, by removing the label fitness term in the formulation and adding an orthogonal constraint to the projection matrix to avoid trivial solution.

Recently, [Meng et al., 2015] proposed the adaptive semisupervised dimensionality reduction (ASSDR), trying to optimize the graph by a heuristic iteration scheme. Two matrices of size n × n need to be stored in each iteration, which is quite memory consuming. ASSDR may rely on kmeans in the second step while kmeans itself is sensitive to initialization. Compared to ASSDR, LMRAG has several advantages: (1) LMRAG has a specific objective function, while ASSDR does not have one since it is based on a heuristic scheme. (2) LMRAG adaptively learns the graph, while ASSDR still uses the pre-defined way in each iteration. (3) The adaptive graph learned by LMRAG is sparse and structured, and the initial graph computed by Eq.(15) is also scale invariant.

5 Experiments

5.1 Datasets

We use several widely used benchmark datasets JAFFE 1 , CMU PIE [Sim et al., 2003], UMIST 2 , YALE, YALE-B 3 , Corel [Chen et al., 2011] and COIL-20 4 to evaluate the proposed LMRAG in our experiments. We provide a brief description of these datasets below.

JAFFE contains 213 images of 7 facial expressions posed by 10 Japanese female models. We used the frontal pose subset (C27) of CMU PIE, in which the images were acquired under variable illuminations and with different expressions. The images in UMIST cover a wide range of poses from profile to frontal views. YALE contains 15 individuals and each one has 11 grayscale images under variable illuminations. YALE-B is an extended version of YALE. Corel has 2074 images, which are represented by color, texture, and shape. COIL-20 is an object dataset that the images were captured from varying angles.

These datasets were first scaled to [0,1] by feature. We cropped UMIST to the size of 56×46. Except for Corel, PCA is then conducted on them with 98% information reserved. The detailed statistics can be seen in Table 1.

5.2 Comparison Algorithms

We compare LMRAG with five existing methods: semisupervised discriminant analysis (SDA) [Cai et al., 2007], trace ratio based flexible SDA (TR-FSDA) [Huang et al., 2012], stable semi-supervised discriminant learning (SSDL)

1 http://www.kasrl.org/jaffe.html

2 http://www.cs.nyu.edu/ roweis/data.html

3 http://www.cad.zju.edu.cn/home/dengcai/Data/data.html

4 http://www.cs.columbia.edu/CAVE/software/softlib/coil20.php

[Gao et al., 2015], FME [Nie et al., 2010] and LapRLS/L [Sindhwani et al., 2005b].

SDA is a representative method that imposes the graph Laplacian regularization into the objective function of linear discriminant analysis (LDA) [Belhumeur et al., 1997]. Based on SDA and FME, TR-FSDA was proposed as the first semisupervised dimensionality reduction method using trace ratio criterion. SSDL considers both the similarity and diversity of data to design the graph, which is then incorporated into the objective function of LDA. We also evaluate the projection ability of LapRLS/L to verify the effectiveness of incorporating LapRLS/L and the adaptive neighbor learning.

Since ASSDR [Meng et al., 2015] is based on pairwise constraints rather than directly uses the label information, to be fair, we do not consider it as a comparison method.

5.3 Experimental Setting

The parameters α and β in LMRAG, SDA, TR-FSDA and SSDL 5 , µ and γ in FME, γA and γI in LapRLS/L need to be tuned, respectively. We searched their values in the range of {10-6 , 10-4 , 10-2 , 100 , 102 , 104 , 106 }. For fair comparison, the reduced dimensionality was fixed as c in SDA, TR-FSDA and SSDL. We randomly chose 40% samples per class as the training data, and used the remaining 60% as the test data. Among the training data, we randomly selected p = {1, 2, 3} samples per class as the labeled data, and used the remaining as the unlabeled data.

To evaluate the projection ability, the nearest neighbor classifier was performed on the projected data for final classification. We uniformly set the neighbor number k to 5 and chose the band width σ of Gaussian kernel in a self-tuning way [Chen et al., 2011] while evaluating the classification performance. We report the best mean accuracy and standard deviation (std) over 20 random splits on each dataset.

5 We applied the Tikhonov regularization to handle the singular problem, thus an additional parameter β is introduced to SSDL.

5.4 Experimental Results

In Figure 1, we tested on JAFFE to give intuitive and practical illustrations of the initial graphs and the adaptive graphs learned by LMRAG with the different neighbor number k. For comparison, the traditional KNN graphs were also illustrated. As can be seen, there are many strong inter-class connections in the KNN graph when k is just a small value 5, and the situation becomes much worse as k increases to 10. We see that the initial graphs of LMRAG are sparse and the sparsity marginally changes with k increasing. With the good initialization and structure constraint, the final adaptive graphs of LMRAG are indeed sparse and structured.

Since parameter γ can be initialized adaptively, and in practice λ can be tuned by a dynamic strategy, we only studied the effect of parameters α and β to the final classification performance on three datasets Corel, CMU PIE and YALE-B. The parameter p was set to 3 during the tests. Figure 2 displays the 3D mesh plots, from which we have the following observations:

1. On each dataset, the performance on the unlabeled training data is basically consistent with the performance on the test data.
2. Comparing to the performance on Corel and YALE-B, the performance on CMU PIE is pretty robust to parameters α and β in a wide range, perhaps because there are more training data in CMU PIE to make up for the accuracy loss caused by inappropriate parameter settting.
3. As parameter α goes up in a range, which means the label fitness term plays a more and more important role in Eq.(7), the performance on all datasets tends to become better. This point is consistent with the first observation from Table 2 listed below. The accuracy may drop when α gets too large, since the model can not make the best use of the unlabeld training data.

Table 2 shows the classification performance in the projected space. Several observations can be made as follows:

1. As the number of labeled samples goes up, the performance of all the methods tends to be better, which demonstrates the usefulness of the labeled data.

1. With respect to the mean accuracy, LMRAG outperforms the other five methods in 40 out of 42 cases, showing the effectiveness of learning an adaptive graph.
2. On four face datasets, when p equals to 1, the performance of LMRAG has a great improvement over others. Specifically, for unlabeled training data, LMRAG exceeds the second best results 6.22%, 7.52%, 7.71% and 13.32% on JAFFE, CMU PIE, UMIST and YALE-B, respectively. For test data, LMRAG exceeds 10.39%, 8.12%, 6.55% and 10.89%, respectively.
3. The superior of LMRAG over LapRLS/L demonstrates the effectiveness of incorporating the adaptive neighbor learning into the objective function of LapRLS/L.

6 Conclusion

In this paper, we have proposed a novel approach denoted by LMRAG which incorporates the graph construction into the semi-supervised dimensionality reduction. The projection matrix and the optimal graph for the specific task are then both obtained. The proposed LMRAG is meaningful as an effective extension and supplement of the existing graph optimized dimensionality reduction methods. Extensive experiments have demonstrated the superiority of LMRAG, comparing to other state-of-the-art methods.

Acknowledgments

This work was supported in part by the National Science Foundation of China under Grants 61522207 and 61473231.