Out-of-samplealgorithmofLaplacianEigenmapsAppliedtoDimensionalityReduction.doc

out-of-sample algorithm of laplacian eigenmaps applied to dimensionality reductionjia peng, yin junsong, huang xinsheng, hu dewendepartment of automatic control, college of mechatronics and automationnational university of defense technologyhunan, china (410073)e-mail：abstractthe traditional nonlinear manifold learning methods have achieved great success in dimensionality reduction. however, when new samples are observed, the batch methods fail tolearn them incrementally. this paper presents out-of-sample extension for laplacian eigenmaps, which computes the low-dimensional representation of data set by optimally preserving localneighborhood information in a certain sense. two different incremental algorithms, thedifferential method and sub-manifold analysis method, are proposed. the algorithms are easy to be implemented and the computation procedure is simple. simulation results testify the efficiency and accuracy of the proposed algorithm.keywords: laplacian eigenmaps, manifold learning, incremental learning, dimensionality reduction1 introduction in pattern recognition and machine learning, dimensionality reduction aims to transform the high-dimensional dataset into a low-dimensional space, while retaining most of the underlying structure in the data. it can be used to solve the curse of dimensionality 1, and to visualize the data set. besides some linear algorithms, such as principal component analysis (pca) 2 and linear discriminant analysis (lda) 3, manifold learning, which serves as a nonlinear method, has attracted more and more attention recently.several representative manifold learning techniques have been proposed. isometric feature mapping (isomap) 4 estimates the geodesic distances on the manifold and uses them for projection. local linear embedding (lle) 5 projects points to a low-dimensional space that preserves local geometric properties. laplacian eigenmaps 6 uses the weighted distance between two points as the loss function to get the dimension reduction results. local tangent space alignment (ltsa) 7 constructs a local tangent space for every sample point and obtains the global low-dimensional embedding results through affine transformation of the local tangent spaces.all of the above algorithms have been widely applied. however, as the batch methods, they are computationally expensive when dealing with large scale problems where both the dimension and the number of training examples are large. to overcome the problem, many researchers have been working on incremental algorithms. literature 8 describes an incremental version of isomap: the geodesic distances are updated, and an incremental eigen-decomposition problem is solved. 9 assumes theeigenvalues of the cost matrix remain the same and tackles the incremental learning problem of lle by- 12 -solving a d dminimization problem, where d is the dimensionality of the low- dimensional space.bengio et al cast mds, isomap, lle and laplacian eigenmaps in a common framework, in which these algorithms are seen as learning eigenfunctions of a kernel, and try to generalize the dimensionality reduction results for the novel data points 10. an incremental manifold learning algorithm via tsa is proposed in 11. the ltsa algorithm is modified, and an incremental eigen-decomposition problem with increasing matrix size is solved by subspace iteration with ritz acceleration.the problem of an incremental algorithm can be stated as follows. assume that the low-dimensional coordinate yi of xi for the first n training samples are given. when a novel sample xn+1 is observed,incremental learning should figure out how to project xn+1 to the low-dimensional space. there are different ideas about incremental methods:(1) increase the dimensionality of eigenspace. some dimensionality reduction algorithms first detect the eigenspace of the data set, and then project the samples into the low-dimensional eigenspace to get the results. when a new sample comes, it is first projected to the current eigenspace and reconstructed. if the reconstruction error is unacceptable by some criterion, the eigenspace dimensionality should be increased. finally, all the low-dimensional points are re-projected to the new eigenspace 12.(2) some methods try to preserve the adjacent information within the data set. they first update adjacent matrix which represents local structure information of the data set. then based on some object function subjected to some certain constraints, the algorithms solve an eigenvalue problem to get the low-dimensional embedding results. when a novel sample xn+1 is observed, the incremental learning procedure can be summarized as: a) update the neighborhood information; b) solve an incremental eigenvalue problem to find the low-dimensional coordinate of xn+1 9, 11.this paper proposes the out-of-sample extension for laplacian eigenmaps, which is traditionallyperformed in a batch mode, thus requiring all training samples to be given in advance. two different incremental algorithms, the differential method and sub-manifold analysis method, are proposed. the algorithms are easy to be implemented and the computation procedure is simple. simulation results testify the efficiency and accuracy of the proposed algorithms.the rest of this paper is organized as follows. section 2 reviews the related work including laplacian eigenmaps algorithm. section 3 presents two different incremental algorithms for laplacian eigenmaps, the differential method and sub-manifold analysis method. section 4 shows some experimental results of the two algorithms mentioned above. section 5 provides some conclusions and discussions.2 laplacian eigenmaps laplacian eigenmaps builds a graph incorporating neighborhood information of the data set 6. using the notion of the laplacian of the graph, the algorithm computes a low-dimensional representation of the data set by optimally preserving local neighborhood information in a certain sense, as shown infigure 1.xiyix jfigure 1: the rationale sketch map of laplacian eigenmaps. when xiandy jxj are far from each other in thehigh-dimensional space,yi andymj are far in the low-dimensional space as well.letx = x1 , x2 ,k, xn be a data set, wherexi r. construct a weighted graph with n nodes, onefor each point, and a set of edges connecting neighboring points. the embedding map is now provided by computing the eigenvectors of the graph laplacian. the algorithm procedure is formally stated below.1.constructing the adjacency graph. let g denote a graph with n nodes. we put an edge betweennodes i and j if xi and xj are close. there are two variations:2(a) -neighborhoods. node i and j are connected by an edge ifnorm is the usual euclidean norm in rn.xi x j , r , where the(b)k nearest neighbors. node i and j are connected by an edge if i is among k nearest neighborsof j or j is among k nearest neighbors of i, k n .2.choosing the weights. here are two variations for weighting the edges. w is a sparse symmetric mx m matrix with wij having the weight of the edge joining vertices i and j, and 0 if there is no such edge.(a)heat kernel. if nodes i and j are connected, putwij = exi x jt2, t r .the heat kernel reflects exactly the distance information between nodes i and j.(b)simple-minded. wij = 1 if and only if vertices i and j are connected by an edge. this kind of weights only describe the simple connection information between nodes i and j.3.construct the object function. consider the problem of mapping the weighted graph g to a lowdimensional space so that connected points stay as close together as possible. lety = ( y1 , y2 ,k, yn )be such a map. a reasonable criterion for choosing a good map is to minimizei jijthe following objective function:i , j ( y y )2 w .there is a heavy penalty if neighboring points xi and xj are mapped far apart. therefore, minimizing the objective function is an attempt to ensure that if xi and xj are close, then yi and yj are close as well.it turns out that for any y, we have1( y y )2w= tr ( yt ly )2i , j i j ijwhere d is diagonal weight matrix, whose entries are column (or row, since w is symmetric) sumsof w,dii= wjij. l = d wis the laplacian matrix.the minimization problem can be reduced to findarg min tr ( yt ly )ys.t. yt dy = 1since the bigger the dii is, the more important yi is, there is a constraint asyt dy = 1.(1)4.eigenmaps. computer eigenvalues and eigenvectors for the generalized eigenvector problem,ly = dy(2)let the column vectorsy0 ,k, yd 1be the solutions of equation (1), ordered according to theireigenvalues,0 = k . we leave out the eigenvectory corresponding to eigenvalue01d 1 00 and use the next d eigenvectors for embedding in d-dimensional euclidean space:xi ( y1 (i),k, yd (i)3 out-of-sample extension for laplacian eigenmaps 3.1 differential method when a novel sample xn+1 is observed, its k nearest neighborhoods can be detected. letx n = xn (1) ,k, xn ( k ) denote this set of points as seen in figure 2, andyn = yn (1) ,k, yn ( k ) thecorresponding low-dimensional coordinates. the out-of-sample procedure is:xn +1xn ( k ) xn ( 2) xn (1) figure 2: a novel sample xn+1 and its k nearest neighborhoodsxn (1) ,k xn ( k ) .1) updating the adjacent matrix w. calculate the connection weights between xn+1 and the n formernodes, and then w is extended to the scale of ( n + 1) ( n + 1) : 2 exp xn +1xi , node i is among k nearest w( n +1)i = wi( n +1) = neighbors of node n+1, i = 1,k, n + 10, others2) calculating the low-dimensional coordinate yn+1 for xn+1. vector yn+1 should satisfy:2 yn +1 yiiw( n +1)it(3)= ( yn +1 yi ) ( yn +1 yi )iw( n +1)idifferentiate (3) by d yn+1, and let the first order derivative equal 0, we get:2 ( yn +1 yi )w( n +1)i = 0i w( n +1)1 i.e. ( y y ,k, y y )m= 0(4)n +11n +1n w ( n +1)n w( n +1)1 ( y1 ,k, yn ) m=obviously,yn +1w ( n +1)n nw( n +1)ii =1(5)3.2 sub-manifold analysis method when a novel sample xn+1 is observed, it is connected directly with its k nearest neighborhoods.x n = xn (1) ,k, xn ( k ) denotes k nearest neighborhoods of xn+1. undoubtedly,xn (1) ,k xn ( k ) , xn +1can be seen as a sub-manifold. letx s = xn (1) ,k xn ( k ) , xn +1 , and the out-of-sample procedure is:1) using laplacian eigenmaps on the sub-manifold. construct the sub-adjacent matrix ws by using the heat kernel: ew (i, j) = 2xi x jt ,t r, xi is within x j s k nearest neighborhoods or x j isxi , x j x ss within xi s k nearest neighborhoods on the sub-manifold0, otherwiseobviously, each point in xs connects with all the other points. calculate the (k+1) by (k+1) matrix dsand sub-laplacian matrix ls:ds (i, i) = ws ( j, i) ,jls = ds ws2) computer eigenvalues and eigenvectors for the generalized eigenvector problem:ls v = ds v(6)let the column vectorsv0 ,k, vd be the solutions of equation (6), ordered according to theireigenvalues,0 = k . leave out the eigenvectorvcorresponding to eigenvalue 0 and01d0use the next d eigenvectors for embedding in d-dimensional euclidean space, and then we get thelow-dimensional coordinates forxn (1) ,k xn ( k ) , xn +1on the sub-manifold:xi (v1 (i),k, vd (i) ,xi x s3) calculating the global coordinate yn+1 for xn+1. from vn+1 to yn+1, the calculating procedure can be seen as a coordinate transformation problem. during the transformation, we preserve the relationshipsbetween xn+1 and its k nearest neighborhoodsxn (1) ,k xn ( k ) . by using laplacian eigenmaps on thesub-manifold, the algorithm has detected the intrinsic structure information between xn+1 andkxn (1) ,k xn ( k ) . so compute constrained weights matrixc = c1 ,k, ck rbetween xn+1 andxn (1) ,k xn ( k ) , with the reconstruction errors minimized:k 2min vk +1 ci vi(7)ci i =1s.t. ci = 1ithus, the optimal weights ci subject to the constraint are founded by solving a least squares problem. then the global coordinate for xn+1 is computed with the weight vector c.kxn+1 yn+1 = ci yn (i ) i =1whereyn (1) ,k yn ( k )are the low-dimensional coordinates for the k nearest neighborhoodsxn (1) ,k xn ( k )of xn+1.4 simulations 4.1 incremental learning in dimensionality reduction this subsection provides incremental learning results in the domain of dimensionality reduction on different datasets.4. 1. 1 resul t s on sw i ss r o l l an d s - cur v e d a t a se tthis section presents the results of laplacian eigenmaps in learning some nonlinear manifolds: swiss roll and s-curve data. the weights matrix w is defined asewij = xi x jt2, t r, node i is among k nearest neighbors of node j ornode j is among k nearest neighbors of node i0, othersi, j = 1, 2,k, nthe number of the nearest neighbors, k, is set to 8. six hundred points are sampled randomly on the manifolds. the first 300 points serve as the training samples, and the rest ones are input in turn totestify the efficiency of the proposed incremental algorithms.figure 3: incremental dimensionality reduction learning on swiss roll data set by using the differential method. (a) the training samples; (b) the entire input samples; (c) the unfolding results of the training samples; (d) incrementallearning results by differential method.figure 4: incremental dimensionality reduction learning on s-curve data set by using the sub-manifold analysis method. (a) the training samples; (b) the entire input samples; (c) the unfolding results of the training samples; (d) incremental learning results by sub-manifold analysis method.in figures 3 and 4, the simulation results on swiss roll and s-curve data set are given, including the training samples and the entire input samples distribution graph, the unfolding results of the training samples and the incremental learning results by the proposed algorithms. as shown in the above figures,the differential and sub-manifold analysis methods both can project the novel samples into the low-dimensional embedding space accurately. when the training samples cover the whole input space well, namely the whole manifold distribution is known, the adjacent information of a novel sample is detected easily, and the low-dimensional embedding result can be computed reliably.4. 1. 2 resul t s on frey f a ce d a tase tfreyface dataset can be available at /roweis/data.html. it contains 1096 face images of one person. each image is 20 28 grayscale. some typical images are shown in figure 5.500 face images serves as the training samples, which are processed by the batch laplacian eigenmaps.600 other face images are then input in sequence. the number of nearest neighbors, k, is set to be 8, and the adjacent matrix w is constructed as follows:1, node i is among k nearest neighbors of node j orwij = node j is among k nearest neighbors of node i0, othersi, j = 1, 2,k, nfigure 5. typical samples from the freyface dataset.the simulation results are shown in figure 6. laplacian eigenmaps dimensionality reduction algorithm detects the intrinsic manifold structure of the freyface dataset. the samples are unfolded in the 2d according to different expressions and gestures. when the 600 face images are observed in sequence, the sub-manifold incremental learning method embeds them to the proper position in the low-dimensional manifold, as can be seen in figure 6 (b). the computation speed is large enough for on-line learning. in the next subsection, we will give the computational complexity analysis of the incremental learning algorithms.(a)(b)figure 6: incremental dimensionality reduction learning on freyface data set by using the sub-manifold analysis method. (a) the training samples; (b) incremental learning results by sub-manifold analysis method.4.2 computational complexity analysis except updating the adjacent matrix w, the core computation process of the differential method consists only of some easy additions and multiplications. the sub-manifold method only needs to solvea ( k + 1) ( k + 1)eigenvector problem and a simple least squares problem, where k is the number ofthe nearest neighbors. the time complexity of solving the ( k + 1) ( k + 1)eigenvector problem iso ( k + 1)3 ) . we use the fmincon function in matlab software to solve the least squares problem,and the computation time is short eno