Compressed Spectral Regression for Efficient Nonlinear Dimensionality Reduction

Introduction

Dimensionality reduction is one of the most useful tools for data analysis in data mining, pattern recognition and many other research fields. Among various dimensionality reduction approaches, spectral dimensionality reduction methods [Belkin and Niyogi, 2001; He and Niyogi, 2003; Yan et al., 2007; Cai, 2009] have received considerable interests in recent years. These methods use information contained in the eigenvectors of a data affinity (i.e., item-item similarity) matrix to reveal the low dimensional structure of the high dimensional data.

The most popular spectral dimensionality reduction algorithms include Locally Linear Embedding [Roweis and Saul, 2000], ISOMAP [Tenenbaum et al., 2000], and Laplacian Eigenmap [Belkin and Niyogi, 2001]. These algorithms only provide the embedding results of training samples. There are many extensions [He and Niyogi, 2003; Yan et al., 2007] which try to solve the out-of-sample problem (i.e., find the embedding of the test samples [Bengio et al., 2003]) by seeking a projective function in a reproducing kernel Hilbert space. However, a disadvantage of these extensions is that their computations usually involve eigen-decomposition of dense matrices which is expensive in both time and memory [Cai, 2009].

Recently, Cai [2009] proposed an efficient spectral dimensionality reduction framework called spectral regression (SR). SR casts the problem of learning a projective function into a regression framework. By avoiding the eigendecomposition of dense matrices, SR provides an unified efficient solution for many supervised and unsupervised spectral dimensionality reduction algorithms [Cai, 2009]. However, SR in unsupervised setting still needs to construct a data affinity matrix and compute the top eigenvectors of the corresponding Laplacian matrix [Chung, 1997]. For a data set consisting of n data points, these steps have a time complexity of O(n2 ). Moreover, when the data are highly nonlinear distributed, SR has to be performed in a reproducing kernel Hilbert space to learn an effective non-linear mapping which leads to kernel SR. As the dense kernel matrix is involved, the time complexity of kernel SR becomes O(n3 ). Such a high computational complexity is an unbearable burden for large-scale applications.

Inspired by the recent progresses on scalable semisupervised learning [Liu et al., 2010] and large scale spectral clustering [Chen and Cai, 2011], we propose an efficient nonlinear dimensionality reduction algorithm termed Compressed Spectral Regression (CSR) in this paper. Specifically, CSR generates l (≪ n) representative points as the landmarks and represent the original data points as the sparse linear combinations of these landmarks. With this landmarkbased sparse representation, the spectral embedding of the data as well as the nonlinear projective function can be efficiently computed. The proposed algorithm scales linearly with the data size. Extensive experiments demonstrate the effectiveness and efficiency of proposed approach.

Background

Suppose we have n data points {xi}ni=1 ⊂ Rm, X = [x1, · · · , xn]. In the past decades, many spectral dimensionality reduction algorithms, either supervised or unsupervised, have been proposed to find a low dimensional representation of xi. Despite the different motivations of these algorithms, they can be nicely interpreted in a general graph embedding framework [Yan et al., 2007; Cai et al., 2007].

These methods first construct an undirected graph G = (V, E) represented by its adjacency matrix W = (wij)ni,j=1 , where wij ≥ 0 denotes the similarity (affinity) between xi and xj . The G and W can be defined to characterize certain statistical or geometric properties of the data set [Yan et al., 2007].

Let y = [y1, y2, · · · , yn]T be the one dimensional embedding. The optimal y∗ is given by solving the following optimization problem:

where the degree matrix D is a diagonal matrix whose entries are column (or row, since W is symmetric) sums of W, Dii = ∑j wij and L = D - W, which is called graph Laplacian [Chung, 1997]. All the three most popular spectral dimensionality reduction algorithms, LLE [Roweis and Saul, 2000], ISOAMP [Tenenbaum et al., 2000] and Laplacian Eigenmap [Belkin and Niyogi, 2001], can be interpreted in this framework with different choices of W [Yan et al., 2007].

The graph embedding approach in Eq. (1) only provides the mappings for the training data. For classification purpose, a mapping for all the samples, including the new test samples, is required. If we choose a linear projective function, i.e., yi = f(xi) = aT xi, then we have y = XT a. Eq. (1) can be rewritten as:

This approach is called Linear extension of Graph Embedding (LGE) [Yan et al., 2007]. With different choices of W, the LGE framework will lead to many popular linear dimensionality reduction algorithms, e.g., Linear Discriminant Analysis, Locality Preserving Projection [He and Niyogi, 2003] and Neighborhood Preserving Embedding [He et al., 2005].

Solving the optimization problem in Eq. (2) involves eigendecomposition of dense matrices XWXT and XDXT , which is time consuming for large scale high dimensional data. To tackle this issue, Cai [2009] proposed a spectral regression (SR) approach. SR learns the projective function in two steps. In the first step, SR solves the optimization problem (1) to get y∗ = [y∗ 1 , y∗ 2 , · · · , y∗n ]T . Since W is usually sparse and D is a diagonal matrix, the optimization problem (1) can be efficiently solved. In the second step, SR solves a regression problem to compute the projective function a∗ :

which can also be efficiently solved via some iterative algorithms (e.g., LSQR [Paige and Saunders, 1982]). It is important to note that the first step becomes trivial in supervised settings and SR (called Spectral Regression Discriminant Analysis in supervised setting) provides an efficient solution for large scale discriminant analysis [Cai et al., 2008]. Please see [Cai, 2009] for more details.

However in unsupervised settings, similar to all the other spectral dimensionality reduction methods, SR needs to construct a data affinity matrix which leads to O(n2 ) computational complexity. Moreover, when the data are highly nonlinear distributed, SR has to be performed in a reproducing kernel Hilbert space (leads to kernel SR) to learn the effective non-linear mapping.

In kernel SR, the linear regression step is replaced by the kernel regression [Cai, 2009]:

where F is the RKHS associated with Mercer kernel K and ‖ ‖K is the corresponding norm. By representer theorem [Wahba, 1990], the solution of the optimization problem (4) can be written as

where K(x, xi) is the kernel function of the corresponding Mercer kernel K. Finally, the optimal α = [α1, α2, · · · , αn]T can be computed as α∗ = (K+δI)-1y, where K is the n×n kernel matrix.

It is clear that the computational complexity of kernel SR is O(n3 ). The high computational costs of the graph construction as well as the kernel approach restrict the applicability of SR for large scale nonlinear problems.

Compressed Spectral Regression

In this section, we introduce our Compressed Spectral Regression (CSR) method for large scale nonlinear dimensionality reduction. The basic idea of our approach is compressing the data using the sparse coding technique [Olshausen and Field, 1996]. With the sparse representation, we can construct a special graph with which the spectral embedding can be efficiently computed. Moreover, the sparse representation capture the nonlinear structure of the data and the learned projective function will be nonlinear with respect to the original features.

Data Compression via Landmark-based Sparse Coding

Sparse coding [Olshausen and Field, 1996; Lee et al., 2006] is a matrix factorization technique which tries to "compress" the data by finding a set of basis vectors and the representation with respect to the basis for each data point. Let X = [x1, · · · , xn] ∈ Rm×n be the data matrix, sparse coding can be mathematically defined as finding two matrices U ∈ Rm×l and Z ∈ Rl×n by solving the optimization problem as follows:

where f is a function measuring the sparsity of each column of Z (e.g., l1 norm) and α is a parameter controlling the sparsity penalty. Each column of U can be regarded as a basis vector which captures the higher-level features in the data and each column of Z is the l-dimensional representation of the original inputs with respect to the new basis. Since each basis vector (column vector of U) can be regarded as a concept, a sparse matrix Z indicates that each data point is a combination of several selected concepts.

Comparing to some dense matrix factorization techniques (e.g., singular value decomposition and nonnegative matrix factorization), sparse coding has several advantages for data representation. The most significant one is that it yields sparse representations such that each data point is represented as a linear combination of a small number of basis vectors. Thus, the data points can be interpreted in a more elegant way.

However, solving the optimization problem (5) is very time consuming. Most of the existing approaches [Lee et al., 2006] compute U and Z iteratively. Apparently, these methods are not suitable for compressing the data efficiently. Following [Liu et al., 2010; Chen and Cai, 2011], we use a much more simplified coding strategy by treating the basis vectors as the landmark points of the data set. A good set of landmarks should be able to cover the data set well. Thus, we use kmeans clustering 1 to generate the landmarks (taking the cluster centers as the landmarks), which gives us the basis matrix U.

To compute the sparse representation matrix Z, we again use a simple yet efficient way. Recall that sparse coding essentially finds two matrices U ∈ Rm×l and Z ∈ Rl×n , whose product can approximate X ∈ Rm×n ,

For any data point xi, its approximation ˆxi can be written as

where uj is j-th column vector of the basis matrix U and zji is ji-th element of Z. A natural assumption here is that zji should be larger if uj is closer to xi. We can emphasize this assumption by setting the zji to zero as uj is not among the r (≤ l) nearest neighbors (for all the column vectors of U) of xi. This restriction naturally leads to a sparse representation matrix Z [Chen and Cai, 2011].

Let N〈i〉 denote the index set which consists r indexes of r nearest landmarks of xi, i.e., uj is among the r nearest landmarks of xi if j ∈ N〈i〉, we compute zji as

where K(·) is a kernel function. One can simply choose the most commonly used Gaussian kernel, K(xi, uj) = exp(-‖xi - uj‖2 /2σ2 ).

1 There is no need to wait the k-means converge and we can stop the k-means after t iterations, where t is a parameter (5 is usually enough). This will be discussed in the experiments.

Spectral Analysis on Landmark-based Graph

Now we have the compressed sparse representation Z ∈ Rl×n for the input X ∈ Rm×n . Instead of constructing a k nearest neighbor graph to model the geometric structure of the data as most of the spectral dimensionality reduction algorithms [He and Niyogi, 2003; Cai, 2009] do, we simply compute the affinity matrix as

where ˆ Z = D-1/2Z and D is a l × l diagonal matrix whose entries are the row sums of Z (dii = ∑j zij).

The advantages of using this graph instead of a k-nearest neighbor graph are as follows:

1. Constructing a k-nearest neighbor graph requires O(n2 ) while computing W in Eq. (8) only needs O(nl2 ). Considering l ≪ n, this is a significant efficiency boosting.
2. The top eigenvectors of W in Eq. (8) can be efficiently computed as we will discuss next.

Noticing that each column of Z sums up to 1, it is easy to check thatn∑

Thus the degree matrix of W in Eq. (8) is I, i.e., the affinity matrix W is normalized [Chung, 1997].

Directly computing the top eigenvectors of W ∈ Rn×n requires O(n2 ) time. By noticing the special structure of W in Eq. (8), we can use an efficient way as follows. Let the Singular Value Decomposition (SVD) of ˆ Z is as follows:

where Σ = diag(σ1, · · · , σl) and σ1 ≥ σ2 ≥ · · · ≥ σl ≥ 0 are the singular values, A ∈ Rl×l and B ∈ Rn×l are the left and right singular vector matrices. It is easy to check that the column vectors of B are the eigenvectors 2 of matrix W = ˆ Z T ˆ Z and the column vectors of A are the eigenvectors of matrix ˆ Z ˆ Z T . Since the size of matrix ˆ Z ˆ Z T is l × l, we can compute A within O(l3 ) time. B can then be computed as

The overall time is O(l3 + l2n), which is a significant reduction from O(n3 ) considering l ≪ n.

Regression with Landmark-based Sparse Representation

The eigenvectors of W in Eq. (8) (the column vectors of B) can be regarded as an approximation of y∗ in the optimization problem (1) with a k-nearest neighbor graph. We then follow the idea of spectral regression [Cai, 2009] to learn the projective functions.

2 It is not hard to verify the the first column of B (i.e., the eigenvector of W corresponding to the largest eigenvalue) is a vector with all ones (without normalization). Thus, we removed this column in real implementation.

For each column vector b in B, we compute the corresponding projective functions p∗ ∈ Rl by solving the regression problem as follows:

where bi is the i-th element in b and α is the ridge regularization parameter.

It is easy to verify that the optimal p∗ has the close form solution as follows:

If we want to reduce the original data into a c-dimensional subspace, the optimal projective function P ∈ Rl×c can be computed as:

where B(c) is the first c columns of B.

The proposed method first compresses the data using landmark-based sparse coding; then performs the spectral analysis on the landmark-based graph; finally, the regression is used to learn the projective functions. Thus, we name this approach Compressed Spectral Regression (CSR).

The final learned CSR model consists two parts: the landmark matrix U ∈ Rm×l and the projective function matrix P ∈ Rl×c . Given a data point x ∈ Rm, CSR model first learns its l-dimensional landmark-base sparse representation z using Eq. (7), then output the c-dimensional reduced representation as ˆx = P T z. The computational procedure of CSR is summarized in Algorithm 1.

It is important to note that CSR is a nonlinear dimensionality reduction method. The mapping from x to z is nonlinear and the mapping from z to ˆx is linear. Thus, the overall mapping from x to ˆx is nonlinear.

Computational Complexity Analysis

Given n data points with dimensionality m, the computational complexity of CSR is as follows:

1. O(tlnm): k-means with t iterations to select l landmarks.
2. O(lnm): learn the landmark-based sparse representation Z.
3. O(l3 + l2n): compute the eigenvectors of W in Eq. (8).
4. O(l2n + l3 + lnc): compute the projection matrix P in Eq. (12).

Considering l ≤ n, the overall complexity of CSR is O(mn) which is linear with respect to number of training samples.

All the existing spectral dimensionality reduction algorithms (e.g., LPP [He and Niyogi, 2003], SR [Cai, 2009], Laplacian Eigenmap [Belkin and Niyogi, 2001], LLE [Roweis and Saul, 2000], ISOMAP [Tenenbaum et al., 2000]) require to construct the affinity graph which is O(mn2 ). Thus, our CSR has a significant computational advantage.

Algorithm 1 Compressed Spectral Regerssion

Input:

• n data points x1, x2, . . . , xn ∈ Rm; Reduced dimensionality c; Parameters l, t, r, α.
• 1: Produce l landmark points using k-means with t iterations, leading to the basis matrix U;
• 2: Learn the landmark-based sparse representation matrix Z ∈ Rl×n according to Eq. (7);
• 3: Compute the first c+1 eigenvectors of ˆ Z ˆ Z T , denoted by A;
• 4: Compute B according to Eq. (9) and remove the first column of B;
• 5: Compute the projection matrix P according to Eq. (12).

Output:

The model: basis matrix U and the projection matrix P ; The dimension reduced data, ˆX = P T Z ∈ Rc×n .

Experiments

In this section, we conduct the clustering experiment on the MNIST (handwritten digits image) data set to demonstrate the effectiveness of the proposed Compressed Spectral Regression (CSR) approach. Each image is represented as a 784 dimensional vector. It has 60000 training images and 10000 test images.

Compared Methods

The compared methods are listed as follows:

• Principle Component Analysis (PCA) [Duda et al., 2000] and Kernel Principle Component Analysis (KPCA) [Scholkopf et al., 1998].
• Locality Preserving Projection (LPP) and Kernel Locality Preserving Projection (KLPP) [He and Niyogi, 2003].
• Spectral Regression (SR) and Kernel Spectral Regression (KSR) [Cai, 2009].
• Laplacian Eigenmap [Belkin and Niyogi, 2001].
• Compressed Spectral Regression (CSR) in this paper.

Given a sample set with k clusters, we use all the above methods to reduce the dimensionality to k. Then the k-means is applied on the reduced subspace. All the methods learn the projective function on the training data and embed the test data using the learned function. The Eigenmap can only provide the embedding of the training data, thus, we do not report its performance on test data.

We use the same 5-nearest neighbor graph in LPP (KLPP), SR (KSR) and Eigenmap. For CSR, we empirically set the parameters l = 1000 (# landmarks), t = 5 (# iterations in k-means), r = 5 (# nearest landmarks) and α = 0.01 (regression regularization). We use the same Gaussian kernel for all the kernel methods and our CSR in Eq. (7). The kernel width parameter σ is empirically 3 set to be 10.

3 we estimate σ by randomly choose 3000 samples and let σ equal to the average of the pair-wise distances.

Table 1: Clustering results on the training set (%)

∗The kernel methods can not be applied due to the memory limitation

Evaluation Metric

The clustering result is evaluated by comparing the obtained label of each sample with the label provided by the data set. We use the normalized mutual information (NMI) metric to measure the clustering performance. The NMI ranges from kmeans 0 to 1 and a higher value indicates a better clustering result. Please see [Cai et al., 2005] for a detailed description on NMI. We also record the running time of each method. All the codes in the experiments are implemented in MATLAB R2011b and run on a Windows 7 machine with 3.40 GHz i72600K CPU, 16GB main memory.

Clustering Results

Table 1, 2 and 3 show the clustering performance on training set, the subspace learning time and the performance on test set, respectively. The evaluations were conducted with the cluster numbers ranging from two to ten. For each given cluster number k (except k = 10), 20 test runs were conducted on different randomly chosen clusters and the average performance as well as the standard deviation are reported in the tables. These results reveals a number of interesting points as follows:

• Comparing to the baseline (k-means in the original feature space), PCA and KPCA do not show any significant improvement while all the other methods do. This demonstrates the effectiveness of the spectral dimensionality reduction methods. The nonlinear spectral methods (KLPP, KSR, CSR and Eigenmap) are significantly better than the linear spectral methods (LPP and SR), especially when the cluster number is large. This probably due to that the data are non-linear distributed. However, due to the memory limitation, all the kernel methods can not be applied when k is larger then 3. The proposed CSR approach achieves the similar performance as the Eigenmap does.
• Considering the learning time, our CSRmethod is linear with respect to the number of samples thus it is very efficient. It only needs 15.7 seconds learning with 10 classes (60000 samples). While all the other spectral methods need to build the affinity graph and need more than 1000 seconds on the same data set.
• Moreover, our CSR model learns the projective function which is defined everywhere (on both training samples and test samples). While the Eigenmap can only provide the embedding results of the training data.

Parameters Selection

Our CSR model has four parameters, # of landmarks (l), # of iterations (t) in k-means to generate the landmarks, # of nearest landmarks (r) for sparse representation and the ridge regularization parameter (α). In this subsection, we will examine how the performance of CSR influenced by these parameters. All the remaining experiments are conducted on the entire MNIST data set (10 classes, 60000 training samples and 10000 test samples). The default settings for CSR are l = 1000, t = 5, r = 5 and α = 0.01. When we study the impact of one parameter, the other parameters are fixed as the default.

Figure (1) shows how the performance of CSR changes as the number of landmarks varies. As the number of landmarks increases, both the performance and the learning time of CSR increase. As we use 3000 landmarks, CSR achieves almost the same performance as Eigenmap does while CSR only needs 44.4 seconds which is much less than the learning time of Eigenmap (1141.4 seconds).

Figure (2) shows how the performance of CSR changes as the number of iterations (t) in k-means varies. t = 0 indicates we use random sampling to select the landmarks. As the number of iterations increases, it is reasonable to see that both the performance and the learning time of CSR increase. Even with random landmark selection, CSR achieves better performance than LPP and SR.

Figure (3) shows how the performance of CSR changes as the number of nearest landmarks (r) for sparse representation varies. As the number of nearest landmarks increases, the performance of CSR consistently decreases. The best performance is achieved when r = 2. Our CSR model is very stable with respect to the regularization parameter α as shown in Table (4).

All the dimensionality reduction algorithms have a critical parameter: the number of reduced dimensions. Figure (4) shows how the performances of different methods vary as the number of reduced dimensions changes. All the algorithms achieve reasonable stable performance as the number of reduced dimensions is between 8 and 20.

Conclusion

In this paper, we have presented a novel efficient nonlinear spectral dimensionality reduction algorithm, called Compressed Spectral Regression (CSR). CSR has O(n) computational complexity where n is the number of samples, which is a significant improvement over all the other spectral methods (require at least O(n2 )). Extensive experiments show the effectiveness and efficiency of our proposed approach.

Acknowledgments

This work was supported in part by the National Basic Research Program of China(973 Program) under Grant 2013CB336500 and National Nature Science Foundation of China (Grant Nos: 61222207, 91120302).