Continuous nonlinear dimensionality reduction by kernel Eigenmaps

1 Background: Graph embeddings

Consider a connected graph with weighted undirected edges specified by edge matrix W. Let W;1 = W1; be the posi­ tive edge weight between connected vertices i and j, zero otherwise. Let D = diag(Wl) be a diagonal matrix where D;; = 2, 1 W;1, the cumulative edge weights into vertex i. The following points are well known or easily derived in spectral graph theory [Fiedler. 1975; Chung, 1997]:

1. The generalized eigenvalue decomposition (EVD)

has real eigenvectors V = [vi, ···, vN] and eigenvalues A= diag([AI 2: A2 2 00 · 2 AN]).

1. Premultiplying equation (1) by n-I makes the general­ ized eigenproblem into a stochastic eigenproblem

where n - I W is a stochastic transition matrix having nonnegative rows that sum to one. The largest eigen­ value of equation ( 1) is therefore stochastic (AI = 1) and its paired eigenvector is uniform (vi = 1/ VN).

1. Expanding and collecting terms in Wu reveals the geo­ metric meaning of the eigenvalues:

The d eigenvectors paired to eigenvalues A2 through Ad+I therefore give an embedding of the vertices in JRd with minimal distortion vis-a-vis the weights, in the sense that a larger W; 1 stipulates a shorter embedding distance. Formally, the embedding

minimizes the distortion d- trace(YWYT)

for any integer d E [ 1, N]. The norm constraint YDY T = I sets the scale of the embedding and causes vertices of high cumulative weight to be embedded nearer to the origin.

1. Y can be rigidly rotated in JRd without changing its dis­ tortion. The distortion measure is also invariant to rigid translations, but the eigenproblem is not, thus there is an unwanted degree of freedom (DOF) in the solution. Due to stochasticity, this DOF is isolated in the uniform eigenvector VI, which is suppressed from the embedding without error (because 1- AI = 0). Adding tvJ toY rigidly translates the embedding bytE JRd .
2. Premultiplying by y T and rearranging equates equa­ tion 1 to the EVD of the graph Laplacian D- W:

1. Premultiplying by n - I/ 2 connects equation to the (symmetric) EVD of the normalized Laplacian:

In summary: Equation 1 gives an optimal embedding of a graph in JRd via eigenvectors v2 · · ·Vd+I; eigenvalue AI is stochastic and the corresponding eigenvector VI is uniform; this is an important property of the EVD solution because it isolates the problem's unwanted translational degree of free­ dom in a single eigenvector, leaving the remaining eigenvec­ tors unpolluted.

Many embedding algorithms can be derived from this anal­ ysis, including the Fiedler vector [Fiedler, 1975], locally lin­ ear embeddings (LLE) [Roweis and Saul, 2000], and Lapla­ cian eigenmaps [Belkin and Niyogi, 2002]. For example, di­ rect solution of equation I gives the Laplacian eigenmap; as a historical note, the symmetrized formulation was proposed by Fiedler in the 1970s and has been used for heuristic graph layout since the 1980s [Mohar, 1991].

2 Transformational embeddings

Now consider a more general problem: We are given some information about the vertices in a matrix Z ~ [z1, · · · ,zN] E !Rd xN, whose columns are generated by applying a vector­ valued function z( ·) __, z E JRd to each vertex of the graph. We seek a linear operator which transforms Z to the optimal graph embedding: G(Z) __, Y. We will call this the "transfor­ mational embedding," to distinguish it from the "direct em­ bedding" discussed above.

A natural candidate for the algebraic statement ofthe trans­ formational embedding problem is the generalized EVD

because setting Y = y Tz makes this equivalent to the orig­ inal direct embedding problem. Again, there is an equiva­ lent symmetric eigenproblem: Make Cholesky 1 decomposi­ tion R T R +- ZDZ T into upper-triangular R E !Rdxd and let

Then

with

This gives an embedding [v2 , v3 , · · ·JTZ, and a computational advantage: IfZ E JRdxN is a short matrix (d « N), the origi­ nal N x N eigenproblem can be reduced to a very small d x d problem, and the matrix multiplications also scale as O(d 2 N) rather than O(N 3 ), due to the sparsity ofW and D.

2.1 Correcting problematic eigenstructure

It is generally the case that yT rf. range(ZT)-there is no lin­ ear combination of the rows ofZ giving Y, so the desired lin­ ear mapping G(Z) __, Y does not exist. Equations 8-11 give the optimal least-squares approximation G(Z) = yTz ~ Y. This approximation can have a serious flaw: If 1 rf. range(ZT) then the first eigenvector v1 is not uniform; it cannot be dis­ carded as the unwanted translational DOF. Worse, all the other eigenvectors will be variously contaminated by the un­ wanted DOF, resulting in an embedding polluted with arti­ facts. For this reason, we call direct solution of equation 8 a raw approximation.

Our options for remedy are limited to those that modify the row-space of Z to reintroduce the uniform eigenvector. For reasons that will become obvious below, we will restrict ourselves to operations that can be applied to any column of Z without knowing any other column.

The simplest such o~eration is to append a uniform row to Z, so that zi __, [zJ, 1] . This makes the relation between Z

1 Any gram-like factorization will work. For example, given EVD Af!A T <- ZDZ T, R =fl 1 1 2 AT. The Cholesky is especially attrac­ tive for its numerical stability, sparsity, and easy invertibility.

and Y affine and guarantees that v J Z is uniform, but it can also force the eigenvectors to model additional variance that is not part of the problem.

Workin~ backward from _the desiderat~ th~t the leading col­ umn ofV Z should be umform, let K = drag(vJ z) - 1 such that ZK is a modified representation of the vertices with val­ ues of z( ·) reweighted on a per-vertex basis: zi __, zd ( v J zi) . Clearly (VTZK)T has a uniform first column, since each row is divided by its first element.

It follows immediately that the related eigenproblem

is stochastic, and Y" ~ [v~, v~ , · .. jT ZK is an embedding with the unwanted translational degree of freedom totally re­ moved. Note that the raw and stochastic approximations are orthogonal (under metric D): Y"DY"T is a diagonal matrix; the other methods are not.

It should be noted that-when scaled to have equal norm trace(YDY T )-none of these approximations has uniformly superior distortion scores; but in Monte Carlo trials with ran­ dom graphs, we find a clear ordering from lowest to highest distortion: reweighted, affine, stochastic, raw (see figure 1 ).

The raw approximation is suboptimal because information about the d-dimensional embedding is spread over d + 1 eigenvectors, no subset of which is optimal. The stochas­ tic approximation is also suboptimal-it optimizes a differ­ ent measure implied by equation 12. In practice, when com­ puting embeddings of graphs whose embedding structure is known a priori, we find that the reweighted and stochastic approximations give results that are clearly very similar, and superior to the other approximations.

The need for any such correction stems from the fact that-the literatures of spectral graph theory and NLDR notwithstanding-equation 1 is not a completely correct statement of the embedding problem. We will show in a forthcoming paper that, as a statement of the embedding problem, equation 1 is both algebraically underconstrained and numerically ill-conditioned. In particular, point #2 is not strictly true: The stochastic eigenvalue is not always paired to a uniform eigenvector. This leads to patholo­ gies that can ruin the embedding, whether obtained from the basic or derived formulations. NLDR algorithms that can be derived from equation 1 (e.g., [Roweis and Saul, 2000; Belkin and Niyogi, 2002; Teh and Roweis, 2003]) do notre­ mediate the problem.

A forthcoming paper makes a full analysis of these is­ sues, identifies the correct problem statements for both equa­ tions 1 & 8, and gives closed-form optimal solutions to both problems. The approximation methods discussed in this sec­ tion are still useful in that they are faster and give reasonably high-quality embeddings. For the NLDR method and datasets considered below, the result of the reweighted approxima­ tion is almost numerically indistinguishable from the optimal embedding, and requires substantially less calculation. The reweighting method can also be justified as a Pade approxi­ mation of the optimal solution.

3 Nonlinear dimensionality reduction

Let X~ [x1, · · ·, XN] E JR 0 be a set of points sampled from a low-dimensional manifold embedded in a high-dimensional ambient space. A reduced-dimension embedding Y ~ [y~, · · ·YN] E JRd with d < D::; N is a set of low-dimensional points with the same local neighborhood structure. We de­ sire instead a mapping G : JR0 -7 JRd, which will general­ ize the correspondence to the whole continuum, with rea­ sonable interpolation and extrapolation to be expected in the neighborhood of the data. Spectral methods for NLDR typ­ ically require solution of many and/or very large eigenvalue or generalized eigenvalue problems [Kruskal and Wish, 1978; Kambhatla and Leen, 1997; Tenenbaum et at., 2000; Roweis and Saul, 2000; Belkin and Niyogi, 2002], and with the ex­ ception of [Teh and Roweis, 2003; Brand, 2003], offer em­ beddings of points rather than mappings between spaces.

Here we show how a leverage the transformational embed­ ding of section 2 into a continuous NLDR algorithm, specifi­ cally a kernel-based mixture of affine maps from the ambient space to the target space. To do so, we must show how the edge weight matrix W and vertex matrix Z are specified. Let Wij ~ f(xi,Xj) iff Xi and Xj satisfy some locality criterion, e.g., llxi -xjll <£,otherwise Wij = 0. As stated above, an embedding Y of X should satisfy

where larger Wij penalize large distances between Yi and Yj·

How should Wij be computed? f is a measure of similar­ ity: The graph-theoretic literature usually takes !(·,-) = 1, while NLDR methods typically take f(xi,Xj) oc exp( -llxi­ Xjll 2 /2cr 2 ) to be a Gaussian kernel, on analogy to heat dif­ fusion models [Belkin and Niyogi, 2002]. The uninforma­ tive setting Wij = 1 is only usable when there is a very large number of points (and edges), so that connectivity informa­ tion alone suffices to determine metric properties of the em­ bedding. The Gaussian setting has a complementary weak­ ness: It can be very sensitive to small variations in distance to neighbors (that may be introduced by the curvature of the data manifold or measurement noise in the ambient space).

f should be monotonically decreasing, relatively insensitive to noise (df should be small), and it should lead to exact re­ constructions of data sampled from manifolds that are already flat. Straightforward calculus shows that equation 13 has the desired minimum when f(xi,Xj) oc llxi -xjll- 1 , or more gen­ erally, the multiplicative inverse of whatever distance mea­ sure is appropriate in the ambient space 2 . (By contrast, the LLE weightings are not correlated with distances.) To make the problem scale invariant, we scale W such that its largest nonzero off-diagonal value is 1 (consequently df::; 1 every­ where f is computed).

Let us now situate some Gaussian kernels Pk(x) ~ !N:(xl,uk.Lk) on the manifold. In this paper, we will take a random subset of data points as kernel centers, and set all Lk = cr 2 1; these kernels are radial basis functions. Let vector

be the kth local homogeneous coordinate of Xi scaled by the posterior of the kth kernel. Ki is an optional local dimensionality-reducing linear projection. Let representation vector

be the vertical concatenation of all such local coordinate vec­ tors. Collect all such column vectors into a basis matrix Z ~ [z(xi), · · · ,z(xN )].

To summarize thus far, our goal now is to find a linear transform Yi ~ G(z(xi)) of the basis (kernel-weighted coor­ dinates) that is maximally consistent with our local distance constraints, specifically

This is isomorphic to the graph embeddings of section 2; the methods developed there apply directly toW and Z. The con­ tinuous mapping from ambient to embedding space immedi­ atly follows from the continuity and smoothness of z( ·):

where the EVD determines the transformation G ~ [v2, · · ·, vd+dT of the continuous kernel representation de­ fined over the entire ambient space:

2 Proof: Consider three points {x1 = 0, x2 =A, x3 = I} on a Io manifold. What similarity measure Wab = f(llxa -xbll) causes the distortion (y2- y!) 2W12 + (y2- y3?Wz3 to have a global minimum at y2 = A? Without loss of generality, we fix the global location and scale of the embedding by fixing the endpoints: {YI = 0, Y3 = I}. Solving for the unique zero of the distortion's first derivative, we obtain the optimum at yz = Wz3/(W12 + Wz3). Since this is a har­ monic relation, the unique continuous satisficing measure is Wab = (llxa -xbl[)- 1 . This sets W12 = I/A and W23 = I/(I-A); some simple algebra confirms that indeed yz = A at the optimum. The induction to multiple points in multiple dimensions is direct.

As a matter of numerical prudence, we recommend using the reweighted approximation:

At first blush, it would seem that reweighting should not be necessary: By construction, 1 E range(ZT), thus vi z(x)­ and the denominator-should be uniform at the datapoints. However, as mentioned above, even when the algebra predicts this structure, numerical eigensolvers may not find it.

To obtain an approximate inverse mapping, we map the means and covariances of each kernel Pk( ·) into the target space to obtain kernels p~(y) ~ 9£(yl,u',:E') there. Then, breaking G = [G,, · · · ,GK] into blocks corresponding to each kernel, take the Moore-Penrose pseudo-inverse of each, and set G+ ~ [G{, · · · ,GtJ. If using the reweighted map, the ap­ proximate inverse map is

4 Illustrative example

We will use a variant of the "swiss roll", a standard test man­ ifold in the NLDR community, to illustrate the arguments and meth­ ods developed in this paper. We sampled a twisted version of the manifold regularly on a 30 x 30 grid and added a small amount of Gaussian isotropic noise. Fig­ ure 2 shows the ideal JR 2 param­ eterization and two views of the ambient JR 3 embedding. Points are shown joined into a line to aid Figure 2: The swiss roll. visual interpretation of the embeddings. All experiments in this

section use a W matrix that was generated using the 12 near­ est neighbors to each point and the inverse distance function.

The Laplacian eigenmap embedding (figure 3) shows the embedding specified by the W matrix. Note that it exhibits some fold­ ing at the comers and top and bottom edges, due partly to problems with the uniform eigenvector and exacerbated by the fact that spectral embeddings tend to compress near the boundaries. The Laplacian of a large 900 x 900 eigen- embedding.

problem, and offers no mapping off the points. Kernel eigen­ maps will be approximations to this embedding.

We now show some kernel eigenmaps computed using the transformational embedding of section 2. All embedding methods are given the same inputs.

Figure 4 shows a raw kernel eigenmap embedding computed using a basis (Z matrix) created from 64 Gaussian unit-cr kernels placed on random points. This required solving a much more manageable 256 x 256 eigenproblem. 100 trials were performed with different sets of ran­ domly placed kernels. In all stochastic maps gave the best reconstructions, while

the raw and affine maps exhibited substantial folding at the edges and comers of the embedding.

Figure 5 shows a reweighted kernel eigen­ map computed from the same W and Z as figures 3 & 4. The result is smoother and actually exhibits less folding than the original Laplacian eigenrnap.

The problem can be reg­ ularized by putting positive mass on the diagonal of W (e.g., W -t W +I), thereby making the recovered ker­ nel eigenrnap more isomet­ ric (bottom figure 5). This regularization is appropri­ ate when it is believed that all neighborhoods are roughly the same size.

The recently proposed Locality Preserving Projec­ tion (LPP) [He and Niyogi, 2002], is essentially the raw approximation (direct solu­ tion of equation 8) with

Wij = e-llx;- xjll 2 /t and Z = X, thereby giving a linear pro­ jection from the ambient space to the target space that best preserves local relationships.

LPP is admirably simple, but it can be shown that the affine approximation from section 2 will always have Jess distortion. LPP can also suffer from Joss of the uni­ shows embeddings of the

swiss roll produced by LPP and by an affine modification of it that is equivalent to our method with a trivial single uniformdensity kernel. Upgrading LPP to an affine projection cap­ tures more of the data's structure. Even so, there is no affine "view" of this manifold that avoids folding.

5 Visualizing word usages

In statistical analyses of natural language, similar usage pat­ terns for two words are taken to indicate that they have sim­ ilar meanings or strongly related meanings. Latent semantic analysis (LSA) is a linear dimensionality reduction of a term­ document co-occurence matrix. The principal components of this matrix give an embedding in which similarly used words are similarly located. Literally, co-location is a proxy for col­ location (the propensity of words to be used together) and synonymy. We may expect that the kernel eigenmap offers a more powerful nonlinear analysis:

The NIPS 12 corpus 3 features a matrix counting occur­ rences of 13000+ words in 1700+ documents. We mod­ eled the first 1 000 words and the last 200 documents in the matrix. This roughly corresponds to one year's papers, a reasonable "snapshot" of the ever-changing terminology of the field. We stemmed the words and combined counts for the same roots, then determined distance between two word roots as the cosines of the angles between their log-domain­ transformed occurrence vectors (x; 1 -> log 2 ( 1 + XiJ) ). The W matrix was generated by adding an edge from each word to its 30 closest neighbors in cosine-space. The representation Z was made using 4 random words as kernel centers. Fig­ ure 7 discusses the resulting 2D embedding, in which techni­ cal terms are clearly grouped by field and many of the more common English words are tightly clustered by common se­ mantics. The first two LSA dimensions (also shown in fig­ ure 7) of the same data are reveal significantly less semantic structure.

6 Discussion

The kernel eigenmap generates continuous nonlinear map­ ping functions for dimensionality reduction and manifold re­ construction. Suitable choices of kernels can reproduce the behavior of several other NLDR methods. One could put a ker­ nel at every local group of points, perform local dimensional­ ity reduction (e.g., a PCA) at each kernel, and thereby obtain from equations 8 and 17 an NLDR algorithm much like chart­ ing [Brand, 2003] or automatic alignment [Teh and Roweis, 2003]. Or, as in the demonstrations above, the kernel eigen­ map can simultaneously determine the local dimensionality reductions and their global merger.

The kernel eigenmap typically substitutes a small dense EVD for the the large sparse EVD of graph embedding prob­ lems. In the sparse case, a specialized power method can compute the desired eigenvectors in significantly less than the O(N 3 ) time required for a full EVD. In the kernel setting, similar efficiencies apply because both Wand Z are typically sparse, allowing fast construction of the reduced EVD prob­ lem ZWZ T; this too is amenable to fast power methods. Of course, the most important efficiency of the kernel method is its ability to embed new points quickly via the function G( x )-there is no need to compute a new global embedding or revise the EVD.

3 Courtesy S. Roweis, available from the U. Toronto website.

The reweighting scheme, although theoretically mooted by our subsequent discovery of a better problem formulation and closed-form solution, is still practically viable as a fast ap­ proximation for large problems, and as a post-conditioning step for unavoidable numerical error of any NLDR algorithm based on eigenvalue decompositions.

In this paper we have used random kernels. There are nu­ merous avenues to discovering stronger methods by investi­ gating placement and tuning of the kernels, stability of the embedding and its topological structure, and sample com­ plexity. In short, all the issues that proved fertile ground for research in classification and regression can be studied anew in the context of estimating the geometry and topology of manifolds.