Chapter 5 Kernel Methods and Manifold Learning

5.1 Kernels and the Kernel Trick

The techniques considered in the previous chapter (PCA, NMF, SVD, and classical Scaling) are ill suited to identify nonlinear structure and dependence in data. If we wish to most efficiently reduce dimensions without loss of information, we will need techniques which incorporate nonlinear structure. One can expand a data matrix by including specific nonlinear relationships then apply PCA or SVD but there are numerous problems with this approach. In particular, which relationships does one choose to include? Even including simple quadratic or cubic terms (features) can result in a data matrix with a massive increase in the number of columns. Even when the original dimensionality of the data is moderate, the including of polynomial terms can quickly result in a data matrix of nonlinear features with an untenable number of columns which can make application of the linear methods we have discussed much more computationally demanding to implement.

Kernels are a important class of functions which can be used to kernelize the methods we have discussed before. In theory, these kernelized versions of the linear methods we have discussed can identify and use nonlinear structure for better dimensionality reduction while circumventing the issue of higher dimensional featurized data. This approach follows from an application of the so called ’kernel trick` which we now discuss.

Briefly, a kernel is a function \[k:\mathbb{R}^d\times \mathbb{R}^d \to \mathbb{R}\] which has an associated feature space, \(\mathcal{H}\) and (implicity defined, possibly nonlinear) feature mapping \(\varphi:\mathcal{R}^d \to \mathcal{H}\) such that inner products in the feature space, denoted \(\langle \varphi(\vec{x}), \varphi(\vec{y})\rangle_{\mathcal{H}}\) can be obtained through an evaluation of the kernel, namely \[\begin{equation} k(\vec{x},\vec{y}) = \langle \varphi(\vec{x}), \varphi(\vec{y})\rangle_{\mathcal{H}} \end{equation}\]

Any method which can be expressed involving inner products can be kernelized by replacing terms of the form \(\vec{x}^T_i\vec{x}_j\) with the quantity \(k(\vec{x}_i,\vec{x}_j)\). Thus, we are replacing inner products of our original \(d\)-dimensional data with inner products in the associated feature space \(\mathcal{H}\). Importantly, if we only need inner products, we never need to explicitly compute the feature map \(\varphi\) for any of our data! At first glance this connection may seem minor, but by using kernels we can turn many linear techniques into nonlinear methods including PCA, SVD, support vector machines, linear regression, and many others.

There are some limits though. Not every choice of \(k\) has an associated feature space. A function is only a kernel if it satisfies Mercer’s Condition.

Theorem 5.1 (Mercer's Condition) A function \[k:\mathbb{R}^d\times \mathbb{R}^d \to \mathbb{R}\] has a an associated feature space \(\mathcal{H}\) and feature mapping \(\varphi:\mathbb{R}^d \to \mathcal{H}\) such that \[k(\vec{x},\vec{y}) = \langle \varphi(\vec{x}), \varphi(\vec{y})\rangle_{\mathcal{H}}, \qquad \forall \vec{x},\vec{y}\in\mathbb{R}^d\] if and only if for any \(N \in \{1,2,\dots\}\) and \(\vec{x}_1,\dots,\vec{x}_N\in\mathbb{R}^d\) the kernel matrix \({\bf K}\in \mathbb{R}^{N}\) with entries \({\bf K}_{ij} = k(\vec{x}_i,\vec{x}_j)\) is positive semidefinite. Equivalently, it must be the case that \[\int_{\mathbb{R}^d}\int_{\mathbb{R}^d} g(\vec{x})g(\vec{y}) k(\vec{x},\vec{y}) d\vec{x}d\vec{y} \ge 0\] whenever \(\int_{\mathbb{R}^2}[g(\vec{x})]d\vec{x}<\infty.\)

We will only consider symmetric functions such that \(k(\vec{x},\vec{y}) = k(\vec{y},\vec{x})\) for all \(\vec{x},\vec{y}\in\mathbb{R}^d\). It may not be immediately obvious if a symmetric function satisfies Mercer’s condition, but there are many known examples. A few are shown in the following table.

Name	Equation	Tuning Parameters
Radial Basis Function	\(k(\vec{x},\vec{y} = \exp\left(-\sigma\|\vec{x}-\vec{y}\|^2\right)\)	Scale \(\sigma >0\)
Laplace	\(k(\vec{x},\vec{y} = \exp\left(-\sigma\|\vec{x}-\vec{y}\|\right)\)	Scale \(\sigma >0\)
Polynomial	\(k(\vec{x},\vec{y}) = (c+ \vec{x}^T\vec{y})^d\)	Offset \(c >0\), Degree \(d \in \mathbb{N}\)

The radial basis function (rbf) is the most commonly used kernel and has an associated feature space \(\mathcal{H}\) which is infinite dimensional! The associated feature map \(\varphi\) for the rbf kernel is \[\varphi(\vec{x}) = e^{-\sigma\|\vec{x}\|^2}\left(a_{\ell_0}^{(0)}, a_{1}^{(1)},\dots,a_{\ell_1}^{(1)}, a_{1}^{(2)},\dots, a_{\ell_2}^{(2)},\dots \right)\] where \(\ell_j = \binom{d+j-1}{j}\) and \(a_\ell^{(j)} = \frac{(2\sigma)^{j/2}x_1^{\eta_1}\dots x^{\eta_d}}{\sqrt{\eta_1!\dots\eta_d!}}\) when \(\eta_1+\dots+\eta_d = j.\) The preceding expression is quite cumbersome, but there is one important point to emphasize. Every possible polynomial combination of the coordinates of \(\vec{x}\) appears in some coordinate of \(\varphi(\vec{x})\) (though higher order terms are shrunk by the factorial factors in the denominator of \(a_\ell^{(j)}\)). Thus, the rbf kernel is associated with a very expressive feature space which makes it a potent but dangerous choice since risks overfitting. To explore these details more, let’s discuss one very important application of kernels in unsupervised learning.

5.2 Kernel PCA

Suppose we have a kernel \(k\) and associated feature map \(\varphi\). In kernel PCA, we want to apply to PCA to the featurized data \(\varphi(\vec{x}_1),\dots,\varphi(\vec{x}_N)\) rather than the original data. The idea is that by studying the featurized data, we can identify additional nonlinear structure in the features that provides a better lower-dimensional representation of the data. We have discussed three approaches to computing PC scores to data: (i) diagonalization of the sample covariance, (ii) applying SVD to the centered data, and (iii) using the duality of PCA and classical scaling.

For the rbf kernel and its infinite dimensional feature map, approaches (i) and (ii) are impossible. Why? The centered data matrix of features \[{\bf H}\tilde{\bf X} = {\bf H}\begin{bmatrix} \varphi(\vec{x}_1)^T \\ \vdots \\ \varphi(\vec{x}_N)^T\end{bmatrix} =\begin{bmatrix} \varphi(\vec{x}_1)^T-\bar{\varphi}^T \\ \vdots \\ \varphi(\vec{x}_N)^T-\bar{\varphi}^T\end{bmatrix} \] has a infinite number of columns so that we cannot compute its SVD. In the above expression, \(\bar{\varphi} = \frac{1}{N} \sum_{i=1}^N \varphi(\vec{x}_i)\) is the mean for the feature vectors. The associated sample covariance matrix \[{\bf \Sigma}_F = \frac{1}{N} \tilde{\bf X}{\bf H}\tilde{\bf X } = \frac{1}{N} \sum_{i=1}^N \left(\varphi(\vec{x}_i) - \bar{\varphi}\right)\left(\varphi(\vec{x}_i) - \bar{\varphi}\right)^T\] will have an infinite number of rows and columns so we cannot hope to diagonalize it either.

Fortunately, the third option, using duality of classical scaling and PC, provides a workaround. Observe that the inner product matrix of the centered feature data \({\bf H}\tilde{\bf X} ({\bf H}\tilde{\bf X})^T\) can be written in terms of the kernel since \[{\bf H}\tilde{\bf X} ({\bf H}\tilde{\bf X})^T = {\bf H} \begin{bmatrix} \varphi(\vec{x}_1)^T \\ \vdots \\ \varphi(\vec{x}_N)^T\end{bmatrix} \begin{bmatrix} \varphi(\vec{x}_1) & \dots & \varphi(\vec{x}_N)\end{bmatrix} {\bf H} = {\bf H K H}\] where \({\bf K}\) has the inner products in the feature space which we can calculate using the kernel function \[{\bf K}_{ij} = \varphi(\vec{x}_i)^T\varphi(\vec{x}_j) = k(\vec{x}_i,\vec{x}_j).\]

Since \(k\) is a symmetric kernel, it follows that \({\bf K}\) is positive semidefinite. Using this property, one can argue that \({\bf HKH}\) will also be positive semidefinite. We can use the eigendecomposition of the doubly centered kernel to compute the kernel principal component scores. Specifically, if \({\bf HKH}\) rank \(r\) with eigenvalues \(\lambda_1\ge \dots \ge \lambda_r >0\) and corresponding eigenvalues \(\vec{u}_1,\dots,\vec{u}_r \in \mathbb{R}^N\), then \({\bf HKH}\) factorizes as \[{\bf HKH} = \underbrace{\begin{bmatrix}\vec{u}_1 & \dots &\vec{u}_r\end{bmatrix} \begin{bmatrix} \lambda_1^{1/2} &0 &0 \\ 0& \ddots & 0 \\ 0 &0 & \lambda_r^{1/2} \end{bmatrix}}_{{\bf U\Lambda}^{1/2}} \left({\bf U\Lambda}^{1/2}\right)^T. \]

The rows of the matrix \({\bf U\Lambda}^{1/2}\) are almost the kernel PC scores. The only issue is an additional the identity \[{\bf HKH} = ({\bf H}\tilde{\bf X})({\bf H}\tilde{\bf X})^T\] is missing the factor of \(1/N\) appearing in the covariance calculation. Accounting for this, the first \(r\) non-zero kernel PC scores are the rows of the matrix \[\frac{1}{\sqrt{N}} {\bf U\Lambda}^{1/2}\] and the corresponding nonzero PC variances are \(\lambda_1/N,\dots,\lambda_r/N.\)

Notably, at no point do we compute the PC loadings! However, similar to standard PCA, we use the scores for dimension reduction and the PC variances for choosing a dimension. Without the loadings, we cannot recompute the original data. Below, we show an application of kernel PCA to the helix and demonstrate its ability to identify the one-dimensional structure of the helix and its sensitivity to kernel selection and tuning.

Example 5.1 (Kernel PCA applied to the Helix) First, we show the kPCA variances for three different kernels and tuning parameters. The data are regularly spaced points along the helix.

From these graphs, one would infer very different lower dimensional choices depending on the kernel and parameters. The polynomial kernel provides the most robust estimate of the one-dimensional nature of the data.

Figure 5.1: kPCA Variances for different Kernels

Below, we show the recovered one-dimensional coordinates for the polynomial kernel with offset 1 and degree 4 shown below, which is good, but do not quite reflect the equal spaced nature of the points.

As the preceding example demonstrates, kernel PCA can identify nonlinear structure, but is quite sensitive to kernel selection and tuning. More advanced implementations make use of cross-validation to aid in the selection and tuning of the kernel [22].

5.3 Brief primer on manifolds and differential geometry

In the preceding section, we focused our attention on linear manifolds and saw cases where this structure was insufficient. Using kernel PCA, we tried to find a workaround by first (implicity) mapping our data to a higher dimensional feature space then approximating results with linear subspaces (of feature space). We now turn to a family of methods that take a more direct approach, assuming the data live on or near a manifold – the manifold hypothesis introduced in Section 3.1 – rather than relying on a fixed kernel. While our data may exhibit some low-dimensional structure, there is no practical reason to expect such behavior to be inherently linear. In the remaining sections of this chapter, we will explore methods which consider nonlinear structure directly. Such methods are referred to as nonlinear dimension reduction or manifold learning. Critical to this discussion is the notion of a manifold.

Definition 5.1 (Informal Definition of a Manifold) A manifold is a (topological) space which locally resembles Euclidean space. Each point on a \(t\)-dimensional manifold has a neighborhood that can be mapped continuously to \(\mathbb{R}^t\). We call \(t\) the intrinsic dimension of the manifold.

We’ll return to mathematical details shortly, but let’s stick with intuition for now. If you were to stand at any point on a \(k\)-dimensional manifold, the portion of the manifold closest to you look just like a \(k\)-dimensional hyperplane – though you might need to be extremely near-sighted for this to be true. The canonical example is the surface of the Earth. If we take all points on the surface of the Earth, they form a sphere in \(\mathbb{R}^3\). However, if we focus on the area around any point it looks like a portion of the two-dimensional plane. In fact, we can represent any point on the surface of the Earth in terms of two numbers, latitude and longitude. In the language of manifolds, we can say the surface of the Earth is a manifold with intrinsic dimension two which has been embedded in \(\mathbb{R}^3\). Here are a few more concrete examples.

Example 5.2 (Examples of Manifolds) Many familiar geometric objects are manifolds such as lines, planes, and spheres. For example, Let \(\vec{w}_1,\dots,\vec{w}_k\in\mathbb{R}^d\) be a set of linearly independent vectors. Then \(\text{span}(\vec{w}_1,\dots,\vec{w}_k)\) is a \(k\)-dimensional manifold in \(\mathbb{R}^d.\) When \(k=1\), the span is a line; for \(k>1\) the span is a hyperplane. The sphere unit sphere \(\{\vec{x}\in\mathbb{R}^d:\,\|\vec{x}\|=1\}\) is a \(d-1\) dimensional manifold. For example, when \(d=3\) manifold resembles the surface of the earth which locally looks like a portion of the 2-dimensional plane. The Swiss roll from Example 3.1 is a two-dimensional manifold in \(\mathbb{R}^3\).

In short manifold as a nice smooth, curved surface. There are several reasons why a subset may not be a manifold. The most straightforward examples are cases where the manifold as a self intersection or a sharp point.

Example 5.3 (Non-manifold) Consider a figure eight curve. At the midpoint where the upper and lower circle meet, there is not neighborhood that resembles Euclidean space so the surface is not a manifold.

For a second example, take any two vectors \(\vec{h}_1\) and \(\vec{h}_2\) and consider there convex hull, that is \(\{a\vec{h}_1+b\vec{h}_2: a,b > 0, a+b\le 1\}\) such the orange triangle in the figure below. At the three vertices of the orange triangle, there is no local neighbor that looks flat. However, if we were to exclude the edges and vertices of the triangle, the subset would be a two-dimensional manifold; thus, \(\{a\vec{h}_1+b\vec{h}_2: a,b > 0, a+b< 1\}\) is a manifold while \(\{a\vec{h}_1+b\vec{h}_2: a,b > 0, a+b\le 1\}\) is not.

5.3.1 Charts, atlases, tangent spaces and approximating tangent planes

Hereafter, we will focus on manifolds which are a subset of \(\mathbb{R}^d\), which are often referred to as submanifolds. Differential geometry can be made far more abstract, but that is unnecessary for the discussion here. After all, we’re dealing with finite dimensional data so any nonlinear surface containing our data must be a submanifold of \(\mathbb{R}^d\)!

Now, some mathematical foundation. Let \(\mathcal{M}\) be a manifold in \(\mathbb{R}^d\) with intrinsic dimension \(t<d\). For every point \(\vec{x}\in\mathcal{M}\), there is a neighborhood \(U_x\subset \mathcal{M}\) containing \(\vec{x}\) and a function \(\phi_x: U_x \to \phi_x(U_x) \subset \mathbb{R}^t\) which is continuous, bijective, and has continuous inverse (a homeomorphism if you like greek). The pair \((U_x,\phi_x)\) is called a chart and behaves much like a map of the area around \(\vec{x}\). There are many choices for \(\phi_x\), but we can always choose one which maps \(\vec{x}\in\mathbb{R}^d\) to the origin in \(\mathbb{R}^t\). For vectors \(\vec{z}\in U_x\), we call \(\phi_x(\vec{z})\in\mathbb{R}^t\) the local coordinates of \(\vec{z}.\)

Example 5.4 (Charts and Manifold Maps) Let’s revisit the Swiss roll example (Example 3.1) where \(A = (\pi/2,9\pi/2)\times (0,15)\). We defined the map \(\Psi:A\to \mathbb{R}^3\) as follows \[\Psi(\vec{z}) = \Psi(z_1,z_2) = \begin{bmatrix} z_1\sin(z_1) \\ z_1\cos(z_1) \\ z_2 \end{bmatrix}.\] For other manifolds, we may need multiple charts to cover the manifold, but we only need one chart the swiss roll (or any manifold defined through a homeomorphic manifold map). Let the neighborhood be the entire manifold, i.e. \(U = \mathcal{M}\) Given \(\vec{x}=(x_1,x_2,x_3)^T\), let \(\phi(\vec{x}) = (\sqrt(x_1^2+x_2^2),x_3)^T.\) This chart essentially unrolls the swiss roll and turns it back into a rectangle in \(\mathbb{R}^2\) so that the local coordinates are the original coordinates!

If we take a collection of charts \(\{U_x,\phi_x\}_{x \in \mathcal{I}}\) such that \(\cup_{x\in\mathcal{I}}U_i = \mathcal{M}\), we have at atlas for the manifold. Here \(\mathcal{I}\) is a subset of \(\mathcal{M}\) which could be countable or finite. With a chart, we can consider doing differential calculus on the manifold. In particular, if \(f:\mathcal{M}\to\mathbb{R}\), then for a chart \((U_x,\phi_x)\), the function \(f\circ \phi_x^{-1}\) is a map from a subset of \(\mathbb{R}^t\) (namely \(\phi_x(U)\)) to \(\mathbb{R}\) so we might hope that we could apply the typical rules of calculus. However, if we have two charts \((U_y,\phi_y)\) and \((U_x,\phi_x)\) which overlap – \((U_x\cap U_y) \ne \emptyset\) – then derivatives \(f\circ \phi_x^{-1}\) and \(f\circ \phi_y^{-1}\) should agree on \(U_x\cap U_y\). More succinctly, the rules of calculus should remain consistent across charts!

Definition 5.2 (Differentiable Manifolds) An atlas \(\{U_x,\phi_x\}_{x\in\mathcal{M}}\) is differentiable if the transition maps \(\phi_x \circ \phi_y^{-1}: \phi_y(U_y) \to \mathbb{R}^t\) are differentiable functions. Recall \(\phi_y(U_y)\subset\mathbb{R}^t\) so differentiable in this case follows the traditional Euclidean definition from calculus.

With some additional properties and manifold with a differential atlas is a differentiable manifold allowing us to compute derivatives of function from the manifold to the reals. We’ll revisit this detail again when discussion Hessian Local Linear Embeddings. For now, let’s assume we have a differentiable manifold with intrinsic dimension \(t\). To every point on the manifold, we can attach a \(t\)-dimensional tangent space. There are several methods for defining the tangent space, but the most straightforward involves the case where we have a manifold map. We’ll restrict our attention to this case.

Definition 5.3 (Tangent Space) Let \(A\subset\mathbb{R}^t\) and let \(\Psi:A\to \mathcal{M}\subset \mathbb{R}^d\). Furthermore, suppose \(\Psi\) has coordinate functions \(\Psi_1,\dots,\Psi_d: A \to \mathbb{R}\) such that for \(\vec{y}=(y_1,\dots,y_t)^T\in A\), \(\Psi(\vec{y}) = (\Psi_1(\vec{y}),\dots,\Psi_d(\vec{y}))^T.\) The Jacobian of \(\Psi\), denoted \({\bf J}_{\Psi}\) is the \(d\times t\) dimensional matrix of partial derivative such that \[({\bf J}_{\Psi})_{ij} = \frac{\partial \Psi_i}{\partial y_j}.\] The tangent space of the manifold \(\mathcal{M}\) at the point \(\vec{p}=\Psi(\vec{y})\), denoted \(T_p(\mathcal{M})\) is the column span of \(({\bf J}_{\Psi})\) evaluated at \(\vec{y}\).

The manifold locally resembles \(\mathbb{R}^t\) so the tangent space should also be \(t\)-dimensional. As such, we’ll require \({\bf J}_{\Psi}\) to be a full rank matrix (with rank \(t\) since \(t < d\)) for every \(\vec{y}\in A\). Importantly, the tangent space is a linear subspace meaning it passes through the origin. This should not be confused with tangent plane to the manifold which we now define.

Definition 5.4 (Tangent Plane) Suppose a manifold \(\mathcal{M}\) has tangent space \(T_p(\mathcal{M})\). Then the approximating tangent plane to the manifold is the affine subspace obtained by translations, namely \(\{\vec{x}\in\mathbb{R}^d: \vec{x}-\vec{p} \in T_p(\mathcal{M})\}.\)

Let’s return to the swiss roll to make these details explicit.

Example 5.5 (Jacobians and Tangent Space for the Swiss Roll) The Jacobian of the swiss roll map is \[{\bf J}_{\Psi} = \begin{bmatrix} z_1\cos(z_1) + \sin(z_1) & 0 \\ -z_1\sin(z_1) + \cos(z_1) & 0 \\ 0 & 1 \end{bmatrix}.\] At the point \(\Psi(3\pi/2,5)=(-3\pi/2,0,5)^T\) the Jacobian is \[{\bf J}_{\Psi}\mid_{(3\pi/2,5)^T} = \begin{bmatrix} -1 & 0 \\ -3\pi/2 & 0 \\ 0 & 1 \end{bmatrix}\] so \[T_{(3\pi/2,0,5)^T}(\mathcal{M}) = \text{span}\{(1,3\pi/2,0)^T, (0,0,1)^T\}.\] We can view the associated approximating tangent plane (translucent blue) at the point \((3\pi/2,0,5)^T\) (in red).

5.3.1.1 Estimating properties from data

In the preceding subsection, we needed to manifold map to define charts, local coordinates, and tangent spaces. However, we do not have access to the manifold map only samples we assume are living near or on an \(t\)-dimensional manifold. In fact, we do not typically know \(t\) either. Fortunately, ideas we have discussed previously can allow us to estimate intrinsic dimensionality, local coordinates, and tangent spaces from data. The key idea is to zoom in on a sufficiently small neighborhood of the manifold that looks inherently Euclidean.

Suppose we have sample \(\vec{x}_1,\dots,\vec{x}_N\in\mathcal{M}\). Given a point \(\vec{x}\in\mathcal{M}\) – typically a sample in our data set – we can first find \(k\) nearest points to \(\vec{x}\) using Euclidean distance. Without loss of generality, let’s call them \(\vec{x}_1,\dots,\vec{x}_k\). If the nearest neighbors are sufficiently close to \(\vec{x}\) they should reside close to a \(t\)-dimensional hyperplane! We can then apply PCA to \(\vec{x}_1,\dots,\vec{x}_k\) or SVD to the displacements \(\vec{x}_1-\vec{x},\dots,\vec{x}_k-\vec{x}\). We expect to see a sharp drop after \(t\) eigenvalues (or singular values) allowing us to estimate the intrinsic dimension \(t\). Subsequently, the first \(t\) PCA scores (or first \(t\) columns of \({\bf US}\) in the SVD) serve as local coordinates for \(\vec{x}_1,\dots,\vec{x}_k\). Finally, the first \(t\) PC loadings (or first \(t\) right singular vectors) are an approximate basis for the tangent space.

Example 5.6 (Estimating Local Properties from Samples) Suppose we are given \(N=10^4\) samples from the Swiss roll. Let’s see how the SVD appraoch outlined above performs at estimating the tangent space at \(\Psi(3\pi/2,5) = (-3\pi/2,0,5)^T=\vec{x}.\) As an example, we’ll plot the singular values of the SVD of the data matrix with rows \(\vec{x}_1^T-\vec{x}^T,\dots, \vec{x}_k^T-\vec{x}^T\). We’ll investigate for a range of \(k\).

The sharp drop after \(t=2\) is striking for \(k >5\) nearest neighbors. Furthermore, this method does an excellent job at estimating the tangent space. Below, we show the associated approximating hyperplane which we estimate using the first two right singular vectors of the data matrix using the first 15 neighbors. The results look visually indiscernible from the case where the Jacobian was used.

SVD and PCA still have their uses thanks to the locally Euclidean nature of manifolds. Naturally, the quality of the estimations depends on having sufficiently dense sampling on the manifold. Verifying results are robust to \(k\) is a good place to start, but this approach has been formalized into a more reliable method for estimating the intrinsic dimensionality of a manifold[23].

5.4 Isometric Feature Map (ISOMAP)

5.4.1 Introduction

The first manifold learning method we are going to cover is the Isometric Feature Map (ISOMAP), originally published by Tenenbaum, de Silva, and Langford in 2000 [24]. As suggested by the name, we will see that the assumption of isometry is central to this method. ISOMAP combines the major algorithmic features of PCA and MDS — computational efficiency, global optimality, and asymptotic convergence guarantees. Thanks to these extraordinary features, ISOMAP is capable of learning a broad class of nonlinear manifolds.

5.4.2 Key Definitions

Different notions of pointwise distance

Prior to discussing the ISOMAP algorithm, let’s briefly discuss the notion of isometry through an example which motivates different notions of distance between two points.

Example 5.7 (Distance between points on a Helix) Consider the helix map \(\Psi:\mathbb{R}\to\mathbb{R}^3\) given by the formula \[\begin{equation} \Psi(t) = \begin{bmatrix} \frac{1}{\sqrt{2}}\cos(t) \\ \frac{1}{\sqrt{2}}\sin(t) \\ \frac{1}{\sqrt{2}}t \\ \end{bmatrix} \end{equation}\] Below, we show the result of applying the Helix map to each point in the interval \((0,25)\). Let’s focus on two points \(\vec{x}_1 = \Psi(2\pi)= (1/\sqrt{2},0,\sqrt{2}\pi)^T\) and \(\vec{x}_2 = \Psi(4\pi)=(1/\sqrt{2},0,2\sqrt{2}\pi)^T\) in particular which are shown as large black dots in the figure below.

There are a few different ways we could measure the distance between the two black points. The first approach would be to ignore the helix (manifold) structure viewing them as vectors in \(\mathbb{R}^3\) and directly measure their Euclidean distance which gives \[\|\vec{x}_1 - \vec{x}_2\| = \sqrt{2}\pi.\] However, we also know that these points are images of the one-dimensional coordinate \(z_1 = 2\pi\) and \(z_2 = 4\pi\) respectively. Thus, we could also consider the Euclidean distance of the lower-dimemsional coordinates which is \(|2\pi - 4\pi| = 2\pi\), which notably differs from the Euclidean distance.

A third option is to return to the three-dimensional representation but to also account for the manifold structure when considering distance. Recall Euclidean distance gives the length of the shortest, straightline path connecting the two points. Instead, let’s restrict ourselves to only those paths which stay on the helix (manifold). You may correctly conclude that the curve starting at \(\Psi(2\pi)\), rotating up the helix one rotation, and ending at \(\Psi(4\pi)\) is the shortest such path. Fortunately, computing arc-length is relatively friendly in this example since \(\Psi\) already parameterizes the path connecting these two points. The arc-length is then \[\int_{2\pi}^{4\pi} \left\|\frac{d\Psi}{dt}\right\| dt = \int_{2\pi}^{4\pi} dt = 2\pi.\] Jumping slightly ahead, we then say the manifold distance between \(\Psi(2\pi)\) and \(\Psi(4\pi)\) is \(2\pi\). Importantly, the manifold distance coincides exactly with the Euclidean distance between the lower-dimensional coordinates. In fact, for any two points, \(s\) and \(t\), on the real line their Euclidean distance, \(|s-t|\) will be the same as the manifold distance between \(\Psi(s)\) and \(\Psi(t)\). Thus, the helix map \(\Psi\) above serves as our first example of an isometric (distance preserving) map.

We may generalize this idea to any smooth manifold to define a new notion of distance. Given a manifold \(\mathcal{M}\), we define the manifold distance function \(d_{\mathcal{M}} : \mathcal{M} \times \mathcal{M} \to [0,\infty)\) as follows

Definition 5.5 (Manifold Distance Function) Given two points \(\vec{x}\) and \(\vec{y}\) on a smooth manifold, \(\mathcal{M}\), let \(\Gamma(\vec{x},\vec{y})\) be the set of all piecewise smooth curves connecting \(\vec{x}\) and \(\vec{y}\) constrained to stay on \(\mathcal{M}\). Then, we define the manifold distance to be \[\begin{equation} d_{\mathcal{M}}(\vec{x},\vec{y}) = \inf_{\gamma \in \Gamma(\vec{x},\vec{y})} L(\gamma) \tag{5.1} \end{equation}\] where \(L(\gamma)\) is the arclength of \(\gamma.\)

As we reviewed above, the helix example with the arclength formula is one example of a manifold and distance function. Additional examples of a manifold and manifold distance include,

• Euclidean space \(\mathbb{R}^d\) where standard Euclidean distance gives the manifold distance.
• The sphere in \(\mathbb{R}^3\) which is a two-dimensional manifold. Its manifold distance is also called the Great Circle Distance.

We may now define the notion of isometry which is a central assumption of ISOMAP.

Definition 5.6 (Isometry) Let \(\mathcal{M}_1\) be a manifold with distance function \(d_{\mathcal{M}_1}\) and let \(\mathcal{M}_2\) be a second manifold with distance function \(d_{\mathcal{M}_2}\). The mapping \(\Psi:\mathcal{M}_1 \mapsto \mathcal{M}_2\) is an isometry if \[d_{\mathcal{M}_1}(x,y) = d_{\mathcal{M}_2}\left(\Psi(\vec{x}),\Psi(\vec{y})\right) \qquad \text{ for all } \vec{x},\vec{y}\in \mathcal{M}_1.\]

For the purposes of ISOMAP, we will think of \(\mathcal{M}_1\) as some subset of a \(\mathbb{R}^k\) for \(k\) small where we measure distances using the Euclidean norm. Then \(\mathcal{M}_2\) will be a \(k\)-dimensional manifold in \(\mathbb{R}^d\) containing our data. Our first assumption is that the manifold mapping \(\Psi\) is an isometry. Unfortunately, in practice we do not know the manifold nor will we have a method for parameterizing curves on the manifold to compute distances.

5.4.3 Algorithm

Instead, ISOMAP makes use of a data-driven approach to estimate the manifold distance between points following a three-step procedure.

1) Construct Weighted Neighborhood Graph:

MDS uses Euclidean distance to measure pairwise distance between points \(\vec{x}_i\) and \(\vec{x}_j\) (data points in space \(\mathcal{M}_2\)), while ISOMAP uses the geodesic distance in order to reveal the underlying manifold structure. However, when the data points in the high dimensional space \(\mathcal{M}_2\) have a manifold structure, usually the Euclidean pairwise distance is quite different from their pairwise geodesic distance. Fortunately, for small distances on a smoothly embedded manifold, the geodesic path between two close-by points lies nearly flat in the ambient space. So, the length of this path will be very close to the straight line (Euclidean) distance between those points in the ambient space.

The key intuition is that as the density of data points on the manifold increases (i.e., points get closer and closer), the straight line segment in the ambient space connecting two neighboring points becomes a better and better approximation of the shortest path between those points on the manifold. In the limit of the density going to infinity, these distances converge.

Let’s elucidate this concept with two illustrative examples. Firstly, imagine a two-dimensional surface, like a Swiss Roll, situated within a three-dimensional space. For an ant journeying across the Swiss Roll, the vast size difference means its immediate surroundings appear flat. From its perspective, the distance between its consecutive steps closely mirrors the distance a human might measure (Euclidean distance) – both virtually equating to the roll’s geodesic distance. For a larger-scale analogy, think of Earth. Suppose extraterrestrial beings possessed technology allowing them to traverse straight through Earth’s crust and mantle, thus following the shortest Euclidean path. Their journey from Los Angeles to New York might save them hundreds of miles compared to humans. However, when moving between closer landmarks, such as the Science Center to the Smith Center, their advantage diminishes.

As a result, when it comes to the measurement of geodesic distance, it is reasonable to only look at those data points that are close to each other. First, calculate all the pairwise Euclidean distance \(d_{ij}=||\vec{x}_i - \vec{x}_j||_2\), then determine which points are neighbors on the manifold by connecting each point to Either (i) All points that lie within a ball of radius \(\epsilon\) of that point; OR (ii) all points which are K-nearest neighbors with it. (Two different criteria, \(K\) and \(\epsilon\) are tuning parameters)

According to this rule, a weighted neighborhood graph \(G = G(V,E)\) can be built. The set of vertices (data points in space \(\mathcal{M}_2\)): \(V = \{\vec{x}_1, \dots , \vec{x}_N\}\) are the input data points, and the set of edges \(E = \{e_{ij}\}\) indicate neighborhood relationships between the points. \(e_{ij} = d_{ij}\) if (i) \(||\vec{x}_i - \vec{x}_j||_2 \leq \epsilon\); OR (ii) \(\vec{x}_j\) is one of the K-nearest neighbors of \(\vec{x}_i\), otherwise \(e_{ij} = \infty\). Sometimes, the tuning of \(\epsilon\) (or \(K\)) is quite decisive in the output of ISOMAP, we will explain this later with a simulation example.

2) Compute graph distances

In this step, we want to estimate the unknown true geodesic distances \(\{d^{\mathcal{M}}_{ij}\}\) between all pairs of points with the help of the neighborhood graph \(G\) we have just built. We use the graph distances \(\{d^{\mathcal{G}}_{ij}\}\)— the shortest distances between all pairs of points in the graph \(G\) to estimate \(\{d^{\mathcal{M}}_{ij}\}\). For \(\vec{x}_i\) and \(\vec{x}_j\) that are not connected to each other, we try to find the shortest path that goes along the connected points on the graph. Following this particular sequence of neighbor-to-neighbor links, the sum of all the link weights along the path is defined as \(\{d^{\mathcal{G}}_{ij}\}\). In other words, we use a number of short Euclidean distances (representing the local structure of the manifold) to approximate the geodesic distance \(\{d^{\mathcal{M}}_{ij}\}\).

This path finding step is usually done by Floyd-Warshall algorithm, which iteratively tries all transit points \(k\) and find those that \(\tilde{d}_{ik} + \tilde{d}_{kj} < \tilde{d}_{ij}\), and updates \(\tilde{d}_{ij} = \tilde{d}_{ik} + \tilde{d}_{kj}\) for all possible combination of \(i,j\). The algorithm works best in dense neighboring graph scenario, with a computational complexity of \(O(n^3)\).

The theoretical guarantee of this graph distance computation method is given by Bernstein et al. [25] one year after they first proposed ISOMAP in their previous paper. They show that asymptotically (as \(n \rightarrow \infty\)), the estimate \(d^{\mathcal{G}}\) converges to \(d^{\mathcal{M}}\) as long as the data points are sampled from a probability distribution that is supported by the entire manifold, and the manifold itself is flat.

The distance matrix \(\Delta\) can be expressed as: \[\Delta_{ij} = d^{\mathcal{G}}_{ij}\]

Simulation Example

Here we provide a randomly generated Neighborhood Graph for six data points, it uses the K-nearest neighbor criteria (can easily tell this since the matrix is not symmetric, \(K=2\))

# Define the matrix
matrix <- matrix(c(
  0,   3,  4,   Inf, Inf, Inf,
  7,   0,  Inf, 2,   Inf, Inf,
  6,   Inf,0,   Inf, 7,   Inf,
  Inf, 5,  Inf, 0,   Inf, 10,
  Inf, Inf,8,   Inf, 0,   13,
  Inf, Inf,Inf, 9,   14,  0
), byrow = TRUE, nrow = 6)
print(matrix)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    0    3    4  Inf  Inf  Inf
## [2,]    7    0  Inf    2  Inf  Inf
## [3,]    6  Inf    0  Inf    7  Inf
## [4,]  Inf    5  Inf    0  Inf   10
## [5,]  Inf  Inf    8  Inf    0   13
## [6,]  Inf  Inf  Inf    9   14    0

Shown below is the implementation of Floyd-Warshall algorithm in R. As you can see from the three for loops, its computation complexity is \(O(n^3)\).

# Adjusting the matrix to set d_ij and d_ji to the smaller value
n <- dim(matrix)[1]

for (i in 1:n) {
  for (j in 1:n) {
    if (i != j && is.finite(matrix[i, j]) && is.finite(matrix[j, i])) {
      min_val <- min(matrix[i, j], matrix[j, i])
      matrix[i, j] <- min_val
      matrix[j, i] <- min_val
    }
  }
}

# Floyd-Warshall Algorithm
floyd_warshall <- function(mat) {
  n <- dim(mat)[1]
  dist <- mat
  
  for (k in 1:n) {
    for (i in 1:n) {
      for (j in 1:n) {
        dist[i, j] <- min(dist[i, j], dist[i, k] + dist[k, j])
      }
    }
  }
  
  return(dist)
}

# Get the result
result <- floyd_warshall(matrix)

# Print the result
print(result)

##      [,1] [,2] [,3] [,4] [,5] [,6]
## [1,]    0    3    4    5   11   14
## [2,]    3    0    7    2   14   11
## [3,]    4    7    0    9    7   18
## [4,]    5    2    9    0   16    9
## [5,]   11   14    7   16    0   13
## [6,]   14   11   18    9   13    0

3) Applying MDS to \(\Delta\)

As mentioned before, ISOMAP can be viewed as the application of classical MDS in a non-linear case. As a result, the reconstruction of \(\{\vec{z}_i\}\) in the \(k\) dimensional \(\mathcal{M}_1\) follows similar steps as that of classical MDS. The main goal is to preserve the geodesic distance of the manifold in \(\mathcal{M}_2\) as much as possible.

Without any additional information, there are infinite \(\{\vec{z}_i\}\) that can be viewed as the optimal solution. For some invertible function \(\Phi:\mathbb{R}^k\to\mathbb{R}^k\), a new manifold mapping \(\Psi \circ \Phi^{-1}\) can be constructed. \(\vec{x}_i = \Psi \circ \Phi^{-1} (\Phi(\vec{z}_i))\), which proves that \(\{\Phi(\vec{z}_i)\}\) is equivalent to \(\{\vec{z}_i\}\) when it comes to the reconstruction of the lower dimensional configuration.

Without loss of generality, we assume that \(\{\vec{z}_i\}\) are actually centered. So the distance matrix of \(\{\vec{z}_i\}\) can be expressed as \(B=Z^T Z\), so that \(B_{ii}=||z_i||^2_2\) and \(B_{ij}={z_i}^T z_j\).

The embedding vectors \(\{\hat{z}_{i}\}\) (estimate of points in lower dimensional feature space \(\mathcal{M}_1\)) are chosen in order to minimize the objective function:

\[\sum_{i,j} \left(\|\vec{z}_i - \vec{z}_j\|_2 - \Delta_{ij}\right)^2\]

Following the same procedure explained in classical MDS chapter, we can compute each entry of \(B\): \[B_{ij}= -\frac{1}{2} \left(\Delta^2_{ij} - \frac{1}{N} \sum^{N}_{k=1} \Delta^2_{kj} - \frac{1}{N} \sum^{N}_{k=1} \Delta^2_{ik} + \frac{1}{N^2} \sum^{N}_{k=1} \sum^{N}_{l=1} \Delta^2_{kl}\right)\]

To express it in matrix form, it is actually, \(B = - \frac{1}{2} H \Delta^{(2)} H\), where \(\Delta^{(2)}\) is the matrix of elementwise-squared geodesic distances, i.e. \(\Delta^{(2)}_{ij} = \Delta_{ij}^2\), and \(H = I_N - \frac{1}{N} \mathbb{1} \mathbb{1}^T\).

The next step is just a PCA problem. Implement eigen decomposition on matrix B, \(B=U \Lambda U^T= (\Lambda^{1/2} U)^T (\Lambda^{1/2} U)\), then arrange the singular value in descending order, find the first \(k^{\prime}\) ones. We acquire \(\Lambda_{k^{\prime}}\) and \(U_{k^{\prime}}\). \[(\hat{z}_1 | \hat{z}_2 | \dots | \hat{z}_N) = \Lambda_{k^{\prime}} U_{k^{\prime}}\]

Since we don’t know the dimension of the underlying feature space, here \(k^{\prime}\) is a tuning parameter. Usually, we use a scree plot (\(k^{\prime}\) against the sum of the omitted eigenvalues ) and find the elbow point.

5.4.4 Limitations of ISOMAP

Though ISOMAP is a powerful manifold learning method that works well under most circumstances. It still has some limitations in certain scenarios.

1. If the noises \(\{\epsilon_i\}\) is not negligible, then ISOMAP may fail to identify the manifold. Also, ISOMAP is quite sensitive to the tuning parameters. To alleviate the negative impact, it’s highly suggested to start with a relatively small \(\epsilon\) or \(K\), and increase them gradually.

2. When data points are sparse in certain areas or directions of the manifold, the integrity of the learned manifold structure can be compromised. The following example will clarify this notion:

3. One of the two major assumptions of ISOMAP is the convexity of the manifold, that is to say, if the manifold contains many holes and concave margins, then the result of ISOMAP will probably be not ideal.

5.5 Locally Linear Embeddings (LLEs)

5.5.1 Introduction

Locally linear embedding (LLE) is an unsupervised learning algorithm first introduced in 2000 by Sam T. Roweis and Lawrence K. Saul [26]. In the original four-page paper, the two authors introduced the LLE algorithm and demonstrated its effectiveness in dimensional reduction, manifold learning, and in handling real-world high-dimensional data. Unlike clustering methods for local dimensional reduction, LLE maps its inputs into a single global coordinate system of lower dimensionality, and its optimizations do not involve local minima. By exploiting the local symmetries of linear reconstructions, LLE is able to learn the global structure of nonlinear manifolds, such as those generated by images of faces or documents of text. Thanks to its great mathematical properties and relatively low computing cost (compared to other manifold learning methods, like ISOMAP), LLE quickly became attractive to researchers after its emergence due to its ability to deal with large amounts of high dimensional data and its non-iterative way of finding the embeddings [27]. Compared to ISOMAP and some other previous manifold learning methods, LLE is computationally simpler and can give useful results on a broader range of manifolds [28].

For dimensional reduction, most methods introduced before LLE need to estimate pairwise distances between even two remote data points, no matter it is the simple Euclidean distance (classical MDS) or more sophisticated manifold distance (ISOMAP). The underlying main idea of these methods is actually finding a configuration that recovers all pairwise distances of original data points as much as possible. LLE, however, is quite different from these previous methods as it focuses on preserving locally linear relationships.

5.5.2 Algorithm

LLE algorithm is actually built on very simple geometric intuitions. As explained in MDS Part, if we consider a small enough region on a manifold in \(D\) dimensional space, in most cases, it can be regarded as a \(d\) dimensional hyperplane (\(d \ll D\)). LLE also makes use of this intuition and assumes that the whole manifold consist of numerous \(d\)-dimensional patches that have been stitched together. Assuming that there exists sufficient data (data points are compact), it is reasonable to expect that each data point and its neighbors lie on or close to a locally linear patch of the manifold.

Following this idea, LLE approximates each data point by a weighted linear combination of its neighbors and proceeds to find a lower-dimensional configuration of data points so that the linear approximations of all data points are best preserved.

Specifically speaking, LLE algorithm consists of three steps. The initial step involves selecting a certain number of each data point’s nearest neighbors based on Euclidean distance. Following this, the second step calculates the optimal reconstruction weights for each point using its nearest neighbors. The final step carries out the embedding while maintaining the local geometry depicted by the reconstruction weights.

5.5.2.1 Construct Neighborhood Graph

This step is actually very similar to that of ISOMAP. The process of finding neighbors in LLE is typically conducted using grouping methods like k-nearest neighbors (KNN) or selecting neighbors within a fixed radius ball (\(\epsilon\)-neighborhoods), based on the Euclidean distance for each data point, in the provided data set. The KNN method is predominantly utilized for its straightforwardness and ease of implementation. The following explanations are based on KNN method.

Denote \(N\) data points in original \(D\) dimensional space as \(\vec{x}_1, \vec{x}_2, \dots, \vec{x}_N \in \mathbb{R}^D\). For a point \(\vec{x}_i, \quad 1 \leq i \leq N\), its neighbor set is defined as \(N_i^k \subseteq \{1, 2, 3, \dots, i-1, i+1, \dots, N \}\), where \(N_i^k\) can also be called as the indices of \(k\) nearest neighbors of \(\vec{x}_i\). The tuning parameter \(k\) is chosen small enough so that the patch around \(\vec{x}_i\) is flat. However, \(k\) should also be strictly larger than \(d\) so as to let the algorithm work.

As we can tell from these, LLE works well only if data points are dense and hopefully evenly distributed, which will be explained in detail in later examples. The parameter tuning of the appropriate number of neighbors, \(k\), faces challenges of complexity, non-linearity, and diversity of high-dimensional input samples. A larger \(k\) value might cause the algorithm to overlook or even lose the local nonlinear features on the manifold. This issue is exacerbated as neighbor selection, typically based on Euclidean distance, can result in distant neighbors when considering the intrinsic geometry of the data, akin to a short circuit. Conversely, an overly small \(k\) value may lead the LLE algorithm to fragment the continuous manifold into isolated local pieces, losing global characteristics.

5.5.2.2 Reconsruct with Linear Weights

As put before, we try to reconstruct each \(\vec{x}_i\) using an almost convex weighted combination of its neighbors. The respective weights of all its neighbors \(\vec{x}_j, \; j \neq i\) for each \(\vec{x}_i\) is quite essential in the later reconstruction of the underlying intrinsic configuration, as we consider these weights to remain invariant before and after mapping.

To explain it in mathematical formulas, the approximate of \(\vec{x}_i\): \(\tilde{x}_i\) is defined as \(\tilde{x}_i = \sum_{j=1}^N w_{ij} \vec{x}_j\). There are two constraints for this formula: First, \(w_{ij} \equiv 0\), if \(j \notin N_i^k\) (consistent with the assumption of \(k\) nearest neighbors); Second, the sum of weights for each \(\vec{x}_i\) is always one, i.e., \(\sum_{j=1}^N w_{ij}=1\).

Then, the problem of finding the optimal \(w_{ij}, \; 1 \leq i,j \leq N\) is equivalent to solving the following constrained Least Squares problem for \(\forall 1 \leq i \leq N\):

\[\begin{equation} \begin{aligned} & \min \left\| \vec{x}_i-\sum_{j \in N_i^k} w_{i j} \vec{x}_j\right\|^2 \\ & \text { s.t. } \quad \sum_{j \in N_i^k} w_{i j}=1 . \end{aligned} (\#eq:opt_prob) \end{equation}\]

It is worth noting that the weights can be negative theoretically, though in practice, we don’t expect that to happen.

Invariance to Rotation, Rescaling and Translation

Define \(\epsilon(w) = \sum_{i=1}^N \left\| \vec{x}_i-\sum_{j \in N_i^k} w_{i j} \vec{x}_j\right\|^2\), which is the cost function.

1. \(\epsilon(w)\) is unchanged by rotation or rescaling by common factor

Actually \(\sum_{i=1}^N \left\| a \text{U} \vec{x}_i-\sum_{j \in N_i^k} w_{i j} a \text{U} \vec{x}_j\right\|^2 = a^2 \epsilon(w)\), where \(a\) is a non-zero scaler and \(\text{U}\) is an orthonormal matrix.

2. \(\epsilon(w)\) is unchanged by translations

Thanks to the constraint that \(\sum_{j=1}^N w_{ij}=1\), for any translation \(\vec{x}_i \rightarrow \vec{x}_i + \vec{y}\), the cost function does not change.

\[\sum_{i=1}^N \left\| (\vec{x}_i + \vec{y}) -\sum_{j \in N_i^k} w_{i j} (\vec{x}_j + \vec{y}) \right\|^2 = \sum_{i=1}^N \left\| \vec{x}_i-\sum_{j \in N_i^k} w_{i j} \vec{x}_j\right\|^2 = \epsilon(w)\]

From the expressions, we develop a strategy that optimizes one row of matrix \(w\) at a time. Now let’s try to rewrite \(\epsilon(\vec{w}_i)=\left\| \vec{x}_i-\sum_{j \in N_i^k} w_{i j} \vec{x}_j\right\|^2\).

\[\begin{align} \epsilon(\vec{w}_i) &= \left\| \vec{x}_i-\sum_{j \in N_i^k} w_{i j} \vec{x}_j\right\|^2 \\ & = \left[ \sum_{j=1}^N w_{ij} (\vec{x}_i - \vec{w}_j) \right]^T \left[ \sum_{l=1}^N w_{il} (\vec{x}_i - \vec{w}_l) \right]^T \\ & = \sum_{j=1}^N \sum_{l=1}^N w_{ij} w_{il} (\vec{x}_i -\vec{x}_j)^T (\vec{x}_i - \vec{x}_l) \\ & = \vec{w}_i^T G_i \vec{w}_i \end{align}\]

\(\vec{w}_i^T = (w_{i1}, w_{i2}, \dots w_{iN})\) is the \(i^{th}\) row of W. Here \(G_i \in \mathbb{R}^{N \times N}\), where entry \(G_{i}(j.l), \; 1 \leq j,l \leq N\) can be represented as:

\[ G_{i}(j,l) = \begin{cases} (\vec{x}_i - \vec{x}_j)^T (\vec{x}_i - \vec{x}_l) & j,l \in N_i^k \\ 0 & j \; or \; l \notin N_i^k \end{cases} \]

The \((j,l)\) entry of \(G_i\) is actually the inner product of \(\vec{x}_j\) and \(\vec{x}_l\) when centered around \(\vec{x}_i\). From this expression, we know that actually \(G_i\) is a sparse matrix and can be reduced to a compact matrix \(\tilde{G}_i \in \mathbb{R}^{k \times k}\) that eliminates those empty columns and rows.

\[\begin{align} \tilde{G}_i & = (\vec{x}_{i[1]} - \vec{x}_i, \dots, \vec{x}_{i[k]} - \vec{x}_i)^T (\vec{x}_{i[1]} - \vec{x}_i, \dots, \vec{x}_{i[k]} - \vec{x}_i) \\ & = Q_i^T Q_i \end{align}\]

where \([1]\) denotes the first entry in \(N_i^k\). So \(\tilde{G}_i\) is actually a real symmetric and positive semi-definite matrix.

Now let’s go back to deal with the optimization function — Equation @ref(eq:opt_prob) can be solved with Lagrange multiplier given that it has only equality constraints. (More details about the use of Lagrange multiplier can be found in Lagrange multiplier)

Optimizing Equation @ref(eq:opt_prob) is equivalent to minimizing (for \(\forall 1 \leq i \leq N\)) \[ f(\vec{w}_i, \lambda) = \vec{w}^T_i G_i \vec{w}_i - \lambda (\vec{w}^T_i \mathbf{1}_k -1) \] which has the result: \[ \vec{w}_i^{\star} = \frac{\tilde{G}^{-}_i \mathbf{1}_k}{\mathbf{1}^T_k \tilde{G}^{-}_i \mathbf{1}_k} \]

Complement:

As discussed before, we can only prove that \(\tilde{G}_i\) is positive semi-definite, however, we cannot ensure that it is positive definite, which means \(\tilde{G}_i\) is not necessarily invertible. That is why we use the generalized inverse sign here. In practice, it can be done through performing SVD on \(\tilde{G}_i\) and select the first few large singular values and eliminate the rest. Then when computing \(\tilde{G}^{-}_i\), just do the reciprocal of these retained singular values.

Actually \(\tilde{G}_i\) only has \(d\) (the intrinsic original dimension if you forget) relatively large eigenvalues. The rest are either very small or zeros. So it is very likely that \(\tilde{G}_i\) is singular, making the computation result highly unstable. In [28], the authors proposed to address this issue through regularizing \(\tilde{G}_i\).

\[ \tilde{G}_i \leftarrow \tilde{G}_i+ \left(\frac{\Delta^2}{k}\right) \operatorname{Tr}\left(\tilde{G}_i\right) \mathbf{I} \] Here \(Tr(\tilde{G}_i)\) denotes the trace of \(\tilde{G}_i\) and \(\Delta \ll 1\).

5.5.2.3 Embedding

In the previous step, we have recovered the optimal weight matrix

\[ \mathbf{W} = \left(\begin{array}{l} \vec{w}_1^T \\ \vdots \\ \vec{w}_N^T \end{array}\right) \]

The optimal weights \(\mathbb{W}\) reflects local, linear geometry around each \(\vec{x}_i\), thus if the configuration \(\vec{y}_1, \vec{y}_2, \dots, \vec{y}_N \in \mathbb{R}^d\) are the lower dimensional representation, they should also “match” the local geometry.

Following this idea, the objective function is:

\[\begin{equation} \begin{aligned} & \underset{\mathbf{Y}}{\text{argmin}} \sum_{i=1}^N \left\| \vec{y}_i - \sum_{j=1}^{N} w_{ij} \vec{y}_j^T \right\|^2 \\ = & \underset{\mathbf{Y}}{\text{argmin}} \sum_{i=1}^N \left\| \sum_{j=1}^N w_{ij} (\vec{y}_i - \vec{y}_j)^T \right\|^2 \\ = & \underset{\mathbf{Y}}{\text{argmin}} \left\| \mathbf{Y - WY} \right\|_F^2 \\ = & \underset{\mathbf{Y}}{\text{argmin}} \left\| (\mathbf{I}_N - \mathbf{W}) \mathbf{Y} \right\|_F^2 \\ = & \underset{\mathbf{Y}}{\text{argmin}} \; \text{Tr} \left[ \mathbf{Y}^T (\mathbf{I}_N - \mathbf{W})^T (\mathbf{I}_N - \mathbf{W}) \mathbf{Y} \right] \end{aligned} \tag{5.2} \end{equation}\] where \(\mathbf{Y}=(\vec{y}_1 | \vec{y}_2 | \dots | \vec{y}_N)^T\)

There are two constraints:

1. \(\mathbf{1}_N^T \mathbf{Y} = \vec{0}\). This forces \(\mathbf{Y}\) to be centered

2. \(\frac{1}{N} \mathbf{Y}^T \mathbf{Y} = \mathbf{I}_d\). This fixes rotation and scaling.

Key Observation

Considering the final expression of Equation (5.2), the optimization function is now equivalent to finding \(\vec{y}_i\)s that minimizes \(\mathbf{Y}^T (\mathbf{I}_N - \mathbf{W})^T (\mathbf{I}_N - \mathbf{W}) \mathbf{Y}\).

Here we introduce \(\mathbf{M} = (\mathbf{I}_N - \mathbf{W})^T (\mathbf{I}_N - \mathbf{W})\), which is a positive semi-definite matrix. Since \(\mathbf{M} \mathbf{1}_N = (\mathbf{I} - \mathbf{W})^T (\mathbf{1}_N - \mathbf{W} \mathbf{1}_N) = \vec{0}\), \(\mathbf{1}_N\) is an eigen-vector of \(\mathbf{M}\) with eigenvalue zero.

\[ \begin{aligned} \mathbf{Y}^T (\mathbf{I}_N - \mathbf{W})^T (\mathbf{I}_N - \mathbf{W}) \mathbf{Y} = \mathbf{Y}^T \mathbf{M} \mathbf{Y} \end{aligned} \] From constraint (b), we know that columns of \(\mathbf{Y}\) are orthogonal to each other. As a result, this whole problem can be simplified to finding the eigen-vectors of \(\mathbf{M}\) with the smallest eigenvalues.

Compute eigen-vectors with the smallest \(d+1\) eigenvalues \(0=\lambda_1 \leq \lambda_2 < \dots < \lambda_{d+1}\), eliminate \(\mathbf{1}_N\) (the first one). The remaining \(d\) vectors are respectively \(\vec{v}_2, \vec{v}_3, \dots \vec{v}_{d+1} \in \mathbb{R}^N\). So \(\mathbf{Y} = (\vec{v}_2 | \vec{v}_2 | \dots | \vec{v}_{d+1})\), we successfully recover the corresponding \(\vec{y}_1, \vec{y}_2, \dots, \vec{y}_N \in \mathbb{R}^d\).

An illustration of the algorithm

In the original paper [26], the authors provide a very intuitive plot that summarizes the above three steps.

LLE_illustration

Parameter Tuning

There are two parameters to tune in LLE, i.e. (the number of neighbors: \(k \,\); the dimension of the recovered configuration: \(d \,\)).

1. For selection of \(d\), we usually use a reverse scree plot and find the elbow point. It is worth noting that we are choosing the smallest \(d+1\) eigenvalues and compute their corresponding eigen-vectors here. Since the eigen-vectors and eigenvalues of a particular matrix is super sensitive to any sort of noises or perturbations, especially for those small eigenvalues, it is hard to accurately derive the corresponding eigen-vectors \(\vec{v}_2, \dots, \vec{v}_{d+1}\). This is called ill-conditioned eigen-problem.

2. Choose the optimal \(k\)

LLE seeks to preserve local structure through nearest neighbor connections. This is the key point to LLE. As a result, we may use the neighbor set of the original \(\vec{x}_1, \vec{x}_2, \dots, \vec{x}_N \in \mathbb{R}^D\) and \(\vec{y}_1, \vec{y}_2, \dots, \vec{y}_N \in \mathbb{R}^d\) as a criteria.

As explained before, we use \(N_i^k\) to denote the indices of k-nearest neighbors to \(\vec{x}_i\). Similarly, we can also use \(V_i^k\) to denote the indices of k-nearest neighbors to \(\vec{y}_i\). They should be as close as possible.

So our objective function here is: \[ Q(k)= \frac{\sum_{i=1}^N \left| N_i^k \cap V_i^k \right|}{Nk} \] Plot \(Q(k)\) against \(k\), select \(k^{\star}\) where the increase of \(Q(k)\) becomes negligible.

5.5.3 Strengths and Weaknesses of LLE

5.5.3.1 Strengths

1. High Computation Efficiency

The low computation cost of LLE algorithm may be its most shining advantage over other manifold learning methods, and it is actually one of its biggest selling points when it was first introduced. The LLE algorithm involves solving a sparse eigen problem, with computational complexity of roughly \(O(N^2 d^2 + N d^3)\) where \(N\) is the number of data points and \(d\) is the dimension of the recovered configuration.

In comparison, ISOMAP requires computing shortest paths between all pairs of points, which is typically done using Dijkstra’s or Floyd-Warshall algorithm, leading to a complexity of \(O(N^2 log N)\) or \(O(N^3)\) respectively. Then, it involves eigen decomposition similar to classical MDS which is \(O(N^3)\).

In practice, \(d \ll N\), hence the computation cost of LLE is lower than that of ISOMAP in most cases.

2. Few parameters to tune

There are only two parameters to tune, respectively the number of neighbors included in the map: \(k\), and the dimension of the recovered configuration: \(d\). In addition, there exist clear methods to find the optimal \(k\) and \(d\), as stated in the previous part. This makes LLE algorithm easy to find the optimal parameters.

5.5.3.2 Weaknesses

1. Sensitivity to tuning parameters

The result of LLE is quite sensitive to its two control parameters: the number of neighbors \(k \,\) and the dimension of the recovered configuration: \(d\).

Here we use the Swiss Roll example to illustrate this. LLE is optimal at \(k=45\). However, when \(k=40\), the recovered lower-dimensional configuration is wrong (Green points and yellow points overlap, which is not the case in the Swiss Roll); and when we slightly increase \(k\) to 50, the recovered two-dimensional expression is not necessarily a rectangle.

2. Vulnerable to sparse or unevenly-distributed samples

The vulnerability towards sparsity and uneven distribution exists in almost all manifold learning methods, including ISOMAP, as we have illustrated in the previous section. LLE is not immune to this either. When a data set is unevenly distributed, since LLE relies on the original Euclidean distance metric, it tends to select neighbors from a singular direction where these neighbors are densely clustered. Clearly, using these selected neighbors to reconstruct the reference point results in significant redundancy in that specific direction. Concurrently, essential information from other directions or regions is not retained for the reconstruction of the reference point. As a result, these selected neighbors are inadequate for accurately representing and reconstructing the reference point. Consequently, much of the intrinsic structure and internal features will be lost after dimension reduction using LLE.

3. Sensitivity to noise

LLE is extremely sensitive to noise. Even a small noise would cause failure in deriving low dimensional configurations. Justin Wang, et al. utilize various visualization examples to illustrate this drawback in their paper [29], you may take a look if you are interested. Various algorithms have been developed to address this issue, i.e., Robustly Locally Linear Embedding (RLLE) [30], and Locally Linear Embedding with Additive Noise (LLEAN) [29]. The former works well when outliers exist, while the latter has a satisfactory performance when the original points are distorted with noises.

5.6 Laplacian Eigenmap (LEs)

The Laplacian eigenmap [31] is method of manifold learning with algorithmic and geometric similarities to LLEs. Like LLEs and ISOMAP, Laplacian eigenmaps make use of \(k\)-nearest neighbor relationships and the solution of an eigenvalue problem to reconstruct the low-dimensional manifold. As suggested by the name, we will be using the graph Laplacian matrix and emphasize the preservation of nearby points on the manifold making Laplacian Eigenmaps a local method with a different emphasis than LLEs.

The graph Laplacian is an important matrix representation of our data which we will revisit later when discussing spectral clustering. In practice, Laplacian eigenmaps use sparse version of the graph Laplacian. For now, we will briefly introduce this matrix without sparsity and an important identity relating the graph Laplacian and the loss function we will minimize when constructing our low-dimensional representation.

Given data \(\vec{x}_1,\dots,\vec{x}_N\), we are going to build a weighted graph \(\mathcal{G}=(V,\mathcal{E}, {\bf W})\), with one node per sample and weighted, undirected edges connecting the nodes. The weights, \({\bf W}_{ij}\ge 0, \, 1\le i,j\le N\), correspond to a notion of affinity, which is typically a decreasing function of the distance between our data. Two common choices are (i) binary weights: \[{\bf W}_{ij} = \begin{cases} 1 & \|\vec{x}_i-\vec{x}_j\| \le \epsilon \\ 0 & \text{ else} \end{cases}\] and (ii) weights based on the radial basis function: \[{\bf W}_{ij} = \exp\left(-\frac{\|\vec{x}_i-\vec{x}_j\|^2}{2\sigma^2}\right). \] The symmetric matrix \({\bf W}\) is called the (weighted) adjacency matrix of the graph, \(\mathcal{G}\), encodes of all the pairwise relationships. We always set the diagonal entries of \({\bf W}\) to zero to preclude self-connections in the graph.

From the adjacency matrix, we can compute the total affinity of each node (sample) to all other nodes (samples) by summing along the rows. Specifically, the \(i\)th entry of the vector \({\bf W}\vec{1}_N\), has the total affinity of the \(i\)th node (\(\vec{x}_i\)) to all other nodes (data). Using these summed affinities, we then create the graph Laplacian \[\begin{equation} {\bf L} = {\bf D} - {\bf W} \end{equation}\] where \({\bf D}\) is a diagonal matrix with entries \({\bf W}\vec{1}_N\) along its diagonal. The graph Laplacian is a symmetric matrix, and while not obvious at first glance, it is also positive semidefinite thanks to the following important identity. Given any vector \(\vec{y}\in\mathbb{R}^N\), \[\begin{equation} \vec{y}^T{\bf L} \vec{y} = \frac{1}{2} \sum_{i=1}^N\sum_{j=1}^N {\bf W}_{ij} (y_i-y_j)^2 \tag{5.3} \end{equation}\] However, it is not full rank. The previous equation shows that \(\vec{1}\) is an eigenvector with eigenvalue 0.

In the previous statement we can view the entries of vector \(\vec{y}\) as \(N\) separate scalars. However, we can also extend the preceding expression to include Euclidean distances between vectors \(\vec{y}_1,\dots,\vec{y}_N\in\mathbb{R}^t\) to give the following important expression \[\begin{equation} \sum_{i=1}^N\sum_{j=1}^N {\bf W}_{ij} \|\vec{y}_i-\vec{y}_j\|^2 = \frac{1}{2}tr\left({\bf Y}^T{\bf LY}\right) \end{equation}\] where \({\bf Y}\in\mathbb{R}^{N\times t}\) has rows \(\vec{y}_1^T,\dots,\vec{y}_N^T.\) We will return to this identity and its implications for Laplacian Eigenmaps after discussing the algorithmic details of the method

5.6.1 Algorithm

5.6.1.1 Compute neighbor relationships

Fix either a number of nearest neighbors \(k > 0\) or maximum distance \(\epsilon > 0\). If using \(\epsilon\), then \(\vec{x}_i\) and \(\vec{x}_j\) (correspondingly nodes \(i\) and \(j\) in the graph) are neighbors if \(\|\vec{x}_i-\vec{x}_j\| \le \epsilon\). Alternatively, if using the nearest neighbor parameter \(k\), then we consider \(\vec{x}_i,\vec{x}_j\) to be neighbors if \(\vec{x}_i\) is one the \(k\) closest points to \(\vec{x}_j\) and \(\vec{x}_j\) is one of the \(k\) closest points to \(\vec{x}_i.\) The construction of neighbors, like the use of pairwise distance alone, results in symmetric neighbor relationship. We then connect nodes \(i\) and \(j\) with an edge if \(\vec{x}_i\) and \(\vec{x}_j\) are neighbors.

5.6.1.2 Compute weights and build graph Laplacian

Nodes which are not connected immediately receive an edge weight equal to zero. For all connected nodes, we compute the edge weight \[{\bf W}_{ij} = \exp\left(-\frac{\|\vec{x}_i-\vec{x}_j\|^2}{2\sigma^2}\right).\] Here we have shown weights based on the radial basis function, which is motivated by theoretical connections to the heat kernel and an approximation of the Laplacian on the manifold \(\mathcal{X}\) [31]. As such, this method is the default in most implementation of Laplacian eigenmaps. The parameter \(\sigma^2\) does require tuning which can have a large impact on the performance of the algorithm.

When \(k\) or \(\epsilon\) are small, which is typically the case in practice, the (weighted) adjacency matrix \({\bf W}\) will be sparse (most entries equal to 0). From this adjacency matrix, we then construct the graph Laplacian as above. The graph Laplacian built from these weights will also be sparse. By preserving only those connections between nearest points, we have only maintained the pairwise, local relationships on the manifold. We will now use \({\bf L}\) to construct a lower-dimensional representation of the data.

5.6.1.3 Solve generalized eigenvalue problem

Consider the loss function \[\begin{equation} \mathcal{L}(\vec{y}_1,\dots,\vec{y}_N) = \sum_{i=1}^N\sum_{j=1}^N {\bf W}_{ij} \|\vec{y}_i-\vec{y}_j\|^2. \end{equation}\] This loss function is most sensitive to large pairwise distance \(\|\vec{y}_i-\vec{y}_j\|\) when \({\bf W}_{ij}\) is also large (our original data were close). Thus, minimizing the preceding penalty prioritizes keeping \(\vec{y}_i\) and \(\vec{y}_j\) close when \(\vec{x}_i\) and \(\vec{x}_j\) have a high affinity (weight). As a result, Laplacian eigenmaps emphasize local geometry.

Vectors \(\vec{y}_1,\dots,\vec{y}_N\) which minimizes this loss function are not unique. First, there is an issue of translation. To address this issue, we will add a constraint that \({\bf D Y}\) is a centered data matrix, i.e. \({\bf Y}^T{\bf D}\vec{1} = \vec{0}\). We can view the matrix \({\bf DY}\) as a reweighting of the data matrix \({\bf Y}\) with higher weights \({\bf D}_{ii}\) for data \(\vec{x}_i\) which are closer to more points. Here \({\bf D}\) is the diagonal matrix used in the definition of the graph Laplacian. Note that \[{\bf DY} = \begin{bmatrix} {\bf D}_{11} \vec{y}_1^T \\ \vdots \\ {\bf D}_{NN}\vec{y}_N\end{bmatrix}.\] Requiring \({\bf DY}\) to be centered results in configurations where those points with highest affinity a constrained close to the origin in our lower dimensional representation.

However, solving the optimization problem with this modified centering constraint is still ill-posed, namely we could take \(\vec{y}_1=\dots=\vec{y}_N = \vec{0}\) giving a configuration which is collapsed onto the origin. In fact, given any configuration \(\vec{y}_1,\dots,\vec{y}_N\) we could decrease the loss by rescaling all of our data by some constant scalar scalar \(0 < c < 1\) since \({\bf D}(c{\bf Y}) = c {\bf DY}\) will still be centered. To address this scaling issue and give a meaningful \(t\)-dimensional configuration, we also add the constraint \({\bf Y}^T {\bf D Y} = {\bf I}_t\). This constraint eliminates the collapse of the \(t\)-dimensional configuration onto a \(t-1\) dimensional hyperplane, and in particular eliminates the cases where the \(1-\)dimensional configuration collapses onto a point.

Thus, we seek a data matrix \({\bf Y}\in\mathbb{R}^{N\times t}\) solving the following constrained optimization problem \[\begin{equation} \text{argmin}_{{\bf Y}^T {\bf D Y} = {\bf I}_t, {\bf Y}^T{\bf D}\vec{1}^T = \vec{0} } = tr({\bf Y}^T {\bf L Y}). \end{equation}\] To solve this problem, we first introduce the change of variable \(\tilde{\bf Y} = {\bf D}^{1/2}{\bf Y}\) so that the constraints become \[{\bf Y}^T{\bf D}\vec{1}^T = ({\bf D}^{1/2} {\bf Y})^T{\bf D}^{1/2}\vec{1}^T = \tilde{\bf Y}^T {\bf D}^{1/2} \vec{1}^T = \vec{0}\] and
\[{\bf Y}^T{\bf D Y}= {\bf Y}^T{\bf D}^{1/2} {\bf D}^{1/2}{\bf Y} = ({\bf D}^{1/2}{\bf Y})^T ({\bf D}^{1/2}{\bf Y}) = \tilde{\bf Y}^T\tilde{\bf Y} = {\bf I}\] implying that the columns of \(\tilde{\bf Y}\) are orthonormal After the change of variable, our optimization problem becomes \[\begin{equation} \text{argmin}_{\tilde{\bf Y}^T\tilde{\bf Y} = {\bf I}_t, \tilde{\bf Y}^{1/2} {\bf D}^{1/2}\vec{1} = \vec{0} } tr\left(\tilde{\bf Y}^T {\bf D}^{-1/2}{\bf L}{\bf D}^{-1/2}\tilde{\bf Y} \right). \end{equation}\] We can minimize this equation by making use of the eigenvalues and eigenvectors of the (symmetric) normalized graph Laplacian \[{\bf L}_{sym} = {\bf I} - {\bf D}^{-1/2}{\bf W}{\bf D}^{-1/2}.\]

Importantly, note that \({\bf L}_{sym}\) is symmetric. Furthermore, it is positive semidefinite since it can be viewed as the graph Laplacian of a graph with weights \({\bf W}_{ij}/\sqrt{{\bf D}_{ii}{\bf D}_{jj}}\) thus subject to the identity (5.3). Thus, it is diagonalizable with orthonormal eigenvectors \(\vec{v}_1,\dots,\vec{v}_N\in\mathbb{R}^N\) and associated nonnegative eigenvalues \(\lambda_1 \le \dots \le \lambda_N\) which we list in increasing order in this case. This leads to the first important observation

1. We could use any \(t\) of these vectors as the columns of \(\tilde{Y}\) and immediately satisfy the constraint \(\tilde{\bf Y}^T\tilde{\bf Y} = {\bf I}.\)

However, we have an additional constraint and the minimization to consider. Note that \(\vec{1}_N\) is an eigenvector of original graph Laplacian with eigenvalue \(0\). Thus, \({\bf L}_{sym}\) also has eigenvalue 0 with associated eigenvector \(\vec{v}_1 = {\bf D}^{1/2}\vec{1}\) since \[{\bf L}_{sym} ({\bf D}^{1/2}\vec{1}) ={\bf D}^{-1/2}({\bf D}^{1/2} - {\bf W D}^{-1/2}){\bf D}^{1/2}\vec{1} = {\bf D}^{-1/2}({\bf D}-{\bf W})\vec{1} = {\bf D}^{-1/2} {\bf L}\vec{1} = \vec{0}.\] As a result, all other eigenvectors of \({\bf L}_{sym}\) must be orthogonal to \(\vec{v}_1 = {\bf D}^{1/2}\vec{1}.\)

2. This suggests that if we drop the first eigenvector associated with eigenvalue \(\lambda_1 = 0\) and use \(t\) the remaining eigenvectors of \({\bf L}_{sym}\) as the columns of \(\tilde{\bf Y}\) we will satisfy the constraint \(\tilde{\bf Y}^T {\bf D}^{1/2}\vec{1} = \vec{0}.\)

The final observation is that we should choose the eigenvectors to minimize the objective.

3. We make use of the eigenvalues themselves and take \[\tilde{Y} = \begin{bmatrix} \vec{v}_2 & \dots & \vec{v}_{t+1} \end{bmatrix}\] so that \[tr\left( \tilde{\bf Y}^T{\bf L}_{sym}\tilde{\bf Y}\right) = \lambda_2+\dots+\lambda_{t+1}\] is minimized.

After undoing the change of variables, we use the rows of \[{\bf Y} = {\bf D}^{-1/2}\tilde{\bf Y}\in \mathbb{R}^{N\times t}\] as our \(t\)-dimensional configuration.

5.7 Hessian Eigenmaps (HLLEs)

5.7.1 Introduction

This method derives its name from the Hessian of a scalar valued function, properties of which we’ll briefly revisit in the context of two examples \[\begin{align*} g_1(x,y) & = x^2+y^2 \\ g_2(x,y) &=\sin(2\pi x)\sin(2\pi y). \end{align*}\] Let’s visualize these functions using surface plots.

Figure 5.2: Surface plots of \(g_1\) (left) and \(g_2\)(right) over the rectangle \([-2,2]^2\)

Which of these two functions is more curved? Informally, which function is wigglier? The answer is clearly \(g_2\), but let’s revisit ideas from calculus to quantify the difference.

Suppose \(g:\mathbb{R}^t\to\mathbb{R}\). The second degree Taylor expansion of a function \(g\) around \(\vec{x}\) is \[g(\vec{z}) \approx g(\vec{x}) + \nabla g(\vec{x})^T(\vec{z}-\vec{x}) + \frac{1}{2}(\vec{z}-\vec{x})^T(\mathcal{H}g(\vec{x}))(\vec{z}-\vec{x}).\] The gradient provides information for orienting the approximating tangent plane of the surface whereas the Hessian provides information on additional variations from the tangent plane which captures curvature of the surface near \(\vec{x}\). For \(g_1\) and \(g_2\) above we have \[\nabla g_1 = \begin{bmatrix}2x \\ 2y \end{bmatrix} \text{ and } \nabla g_2 = \begin{bmatrix} 2\pi\cos(2\pi x)\sin(2\pi y) \\ 2\pi\sin(2\pi x)\cos(2\pi y) \end{bmatrix}\] and \[\mathcal{H} g_1 = \begin{bmatrix} 2 & 0 \\ 0 & 2\end{bmatrix} \text{ and } \mathcal{H} g_2 = \begin{bmatrix} -4\pi^2\sin(2\pi x)\sin(2\pi y) & 4\pi^2 \cos(2\pi x)\cos(2\pi y) \\4\pi^2 \cos(2\pi x)\cos(2\pi y) & -4\pi^2 \sin(2\pi x)\sin(2\pi y)\end{bmatrix}.\] The Hessian of \(g_1\) is constant whereas the Hessian of \(g_2\) is both varying but generally comprised of larger values.

Subsequently, the squared Frobenius norms of the Hessians are \[\|\mathcal H g_1\|^2_F = 8 \text{ and } \mathcal{H} g_2 = 16\pi^4\left[1 + \cos(4\pi x)\cos(4\pi y)\right]\] where the later expression follows after some trig identities. Importantly, the squared Frobenius norm of the Hessian for \(g_2\) ranges from \(0\) to \(32\pi^4\). On average (over the rectangle \([-2,2]^2\)) \(\|\mathcal{H}g_1\|_F^2\) is much larger than \(\|\mathcal{H} g_1\|_F^2\)! Mathematically, we have \[\int_{-2}^2 \|\mathcal{H} g_2 \|_F^2 dx dy = 16\pi^2 > 8 = \int_{-2}^2\|\mathcal{H} g_1\|_F^2 dxdy.\]

To adapt the notion of gradients and Hessians to manifolds, let us revisit the idea of derivatives of real valued functions on a manifold. Suppose we have a manifold \(\mathcal{M} \subset \mathbb{R}^d\) and function \(f:\mathcal{M}\to\mathbb{R}\). Assuming that \(\mathcal{M}\) has intrinsic dimension \(t < d\), to each point \(\vec{x}\in\mathcal{M}\), we can associate a \(t\) dimensional tangent space \(T_{\vec{x}}(\mathcal{M})\) which we equip with a choice of orthonormal basis vectors and a corresponding chart \((U,\phi\) which maps a neighborhood, \(U\) containing \(\vec{x}\) to \(\mathbb{R}^t\). For now, it is sufficient to think of the \(\phi\) as the projection of neighboring points of \(\vec{x}\) onto the approximating tangent plane (if we then subtract \(\vec{x}\) from these projections, we get projections onto the tangent space). For \(\vec{z}\in U\), we referred to \(\phi(\vec{z})\) as the local coordinates of \(\vec{z}\).

From the perspective of calculus, the function \(f \circ \phi^{-1}\) is a function mapping an open subset of \(\mathbb{R}^t\) to \(\mathbb{R}\). We can then compute gradients and Hessians of \(f \circ \phi^{-1}\) using standard multivariable calculus and \(\nabla^{tan}f:=\nabla (f\circ \phi^{-1})\) and \(\mathcal{H}^{tan}(f) := \mathcal{H}(f\circ \phi^{-1})\). We’ll refer to \(\nabla^{tan}f\) and \(\mathcal{H}^{tan}f\) as the gradient and Hessian of \(f\) in local coordinates. For a choice of \(f\), the specific values in gradient or Hessian depends on our choice of \(\phi^{-1}\) (or equivalently the basis we choose for the tangent space). However, the (squared) Euclidean length of the gradient or the (squared) Frobenius norm of the Hessian do not depend on this choice since these quantities are unchanged by orthonormal transformations. We’ll being using this fact for the remainder of this section.

In principal, we could take every point \(\vec{x}\in\mathcal{M}\) and, using this approach, compute the local Hessian of \(f\) at \(\vec{x}\). If we take an average of the local Hessian over the manifold, denoted \[\mathbb{H}(f) = \int_{\mathcal{M}}\|\mathcal{H}^{tan} f \|_F^2\, d\mu(\vec{x})\] we have a measure of the average curvature of the function \(f\). Here, \(\mu(\vec{x})\) is the probability distribution of data on the manifold fold (a measure for those with the requisite background). The utility of this quantity of \(f\) is given in the following theorem.

Theorem 5.2 (Null Space of Hessian Operator on Manifolds) Suppose \(\mathcal{M} = \Psi(A)\) where \(A\) is an open, connected subset of \(\mathbb{R}^t\), and \(\Psi\) is a locally isometric embedding of \(A\) into \(\mathbb{R}^d\). Then \(\mathbb{H}(f)\) has a \(t+1\) dimensional null space spanned by the constant function and a \(t\)-dimensional space of functions spanned by the original coordinates in \(A\).

In short, if \(\mathbb{H}(f) = 0\) then there is some \(\alpha\in\mathbb{R}\) and \(\vec{\beta}\in\mathbb{R}^t\) such that \(f(\vec{x}) = \alpha + \beta^T\vec{z}\) where \(\vec{z}\in A\) is the preimage of \(\vec{x}\), e.g. \(\Psi(\vec{z}) = \vec{x}.\) By looking for functions in the null space of \(\mathbb{H}(f)\), we can then find functions which help us recover the original low-dimensional coordinates in \(A\)!

Unfortunately, we cannot calculate \(\mathbb{H}f\) for a general function \(f\). However, we can estimate \(\mathbb{H}(f)\) using observed data \(\vec{x}_1,\dots,\vec{x}_N\in\mathcal{M}\), which is the central idea behind HLLEs. Specifically, we will construct a discretized estimate of \(\mathbb{H}(f)\) which we can write in the quadratic form \[\vec{f}^T{\bf H}\vec{f} \text{ where } \vec{f} = \left(f(\vec{x}_1),\dots,f(\vec{x}_N)\right)^T\in\mathbb{R}^N.\] The matrix \({\bf H}\in\mathbb{R}^{N\times N}\) is symmetric and positive semidefinite and depends only on the observed data \(\vec{x}_1,\dots,\vec{x}_N\). Its construction will be discussed in greater detail in the algorithms section below.

For now, assume that we have access to \({\bf H}\). Suppose that \(f\) is in the null space of \(\mathcal{H}(f)\) then i) \(\vec{f}^T{\bf H}\vec{f} \approx 0\) and ii) from Theorem 5.2 it must be the case that \(f(\vec{x}_i) = \alpha + \beta^T\vec{z}_i\) where \(\vec{z}_i\) is a point in the lower dimensional configuration \(A\). Importantly, if we can identify the null space of the matrix \({\bf H}\) then we can use this to recover the original low-dimensional coordinates (up to rigid motion and coordinate rescaling).

Similar to LLEs and ISOMAP, this algorithm begins with a k-nearest neighbor search and finishes with a eigenvalue decomposition. As inputs we provide the data, the intrinsic dimension, \(t\), of the manifold, and a specified number of nearest neighbors \(k\). For reasons we’ll discuss later, we must select \(k > t(t+3)/2\). This constraint can be problematic in practice. For example, if we believe the intrinsic dimension of \(\mathcal{M}\) is 100, then we need to pick \(k > 5150\) meaning we need at least 5151 samples in our dataset! Such a constraint may be difficult to satisfy in practice, but for dimension reduction focused on visualization with the ambition assumption that \(t=1,2,\) or \(3\), \(k \ge 10\) is sufficient.

5.7.1.1 Compute nearest neighbors and local coordinates

For each \(\vec{x}_i\) compute its \(k\)-nearest neighbors which we’ll denote \(\mathcal{N}_i\). We’ll use each of these neighbors to estimate the Hessian of \(f\) at \(\vec{x}_i\). Next we’ll estimate local coordinates of the neighbors in \(\mathcal{N}_i\) by using SVD to approximate the tangent plane \(T_{\vec{x}_i}(\mathcal{M})\). One important observation to recall. We expect the points in \(\mathcal{N}_i\) to all be close together and centered around \(\vec{x}_i\). In this neighborhood around \(\vec{x}_i\) we expect the manifold \(\mathcal{M}\) to look like a small patch of \(\mathbb{R}^d\). As such, they should be (nearly) contained in a \(t\)-dimensional affine subspace of \(\mathbb{R}^d.\)

Let \({\bf M}_i \in \mathbb{R}^{k\times d}\) be the matrix of centered nearest neighbors of \(\vec{x}_i\). Specifically, the \(j\)th row of \({\bf M}_i\) is \((\vec{x}_{i_j} - \vec{x}_i)^T\) where \(\vec{x}_{i_j}\) is the \(j\)th nearest neighbor of \(\vec{x}_i\). We apply a singular value decomposition to \({\bf M}_i\) giving factorization \[{\bf M}_i = {\bf U}_i {\bf S}_i {\bf V}_i^T\] where \({\bf U}_i\in\mathbb{R}^{k\times k}\) and \({\bf V}_i\in\mathbb{R}^{d\times d}\) have orthonormal columns and \({\bf S}_i\in\mathbb{R}^{k\times d}\) is diagonal. If the neighbors in \(\mathcal{N}_i\) were fully contained in a t-dimensional affine subspace we expect that \({\bf S}_i\) has \(t\) large singular values with the remainder equal to zero. Additionally in this case, the tangent space \(T_{\vec{x}_i}\) is \(\vec{x}_i + \text{span}\{\vec{v}_1,\dots,\vec{v}_t\}\) where \(\vec{v}_1,\dots,\vec{v}_t\) are the first \(t\) columns of \({\bf V}\).

In practice, we should only expect the first \(t\) singular values to be large with the remainder smaller but non-zero. However, we’ll still use the first \(t\) columns of \({\bf V}\) as an (approximate) orthonormal basis for the tangent space. The first \(t\) columns of \({\bf U}_i{\bf S}_i\) give approximations of the local coordinates of the neighbors \(\mathcal{N}_i\) (the \(j\)th row gives the local coordinates for the \(j\)th nearest neighbor). Hereafter, we’ll let \(\vec{u}^i_j = (u^i_{j1},\dots,u^i_{jt})^T\) denote the local coordinates of the \(j\)th nearest neighbor.

5.7.1.2 Estimate the Hessian

Recall from section 5.3, that we can define derivatives for a function \(f\) at \(\vec{x}_i\) using the local coordinates.

Now, we can use the local gradients and Hessians to construct a Taylor approximation to \(f(\vec{x}_j)\) for each point in the neighbor \(\mathcal{N}_j\). In the tangent space \(T_{\vec{x}_i}(\mathcal{M})\), we associate the origin with \(\vec{x}_i\). Taylor expanding around the origin to second order we have the following approximation \[\begin{align*} f(\vec{x}_{i_j}) &\approx f(\vec{x}_i) + \left[\nabla^{tan} f(\vec{x}_i)\right]^T \cdot \vec{u}^i_j + \frac{1}{2}(\vec{u}^i_j)^T \left(H^{tan}f (\vec{x}_i)\right)\vec{u}^i_j\\ &= f(\vec{x}_i) + \sum_{\ell=1}^t \left[\nabla^{tan} f(\vec{x}_i)\right]_\ell \vec{u}^i_{j\ell} + \frac{1}{2}\sum_{\ell=1}^t \left(H^{tan}f (\vec{x}_i)\right)_{\ell\ell} (u^i_{j\ell})^2 + \sum_{\ell < s} \left(H^{tan}f (\vec{x}_i)\right)_{\ell,s} u^i_{j\ell}u^i_{js} \end{align*}\] To estimate the entries of the Hessian, we use quadratic regression. Build design matrix \({\bf X}_i \in \mathbb{R}^{k\times (1+d+d(d+1)/2)}\) including all terms up to second order of the local coordinates such that \[{\bf X}_i = \begin{bmatrix} 1 & u^i_{11} & \dots u^i_{1t} & \frac{1}{2}(u^i_{11})^2 & \dots& \frac{1}{2}(u^i_{1t})^2 & \frac{\sqrt{2}}{2}u^i_{11}u^i_{12}& \frac{\sqrt{2}}{2} u^i_{11}u^i_{13} & \dots &\frac{\sqrt{2}}{2}u^i_{1,t-1}u^i_{1t} \\ \vdots & \vdots & \vdots & \vdots & \dots & \vdots & \vdots & \vdots & \dots & \vdots\\ 1 & u^i_{k1} & \dots u^i_{kt} & \frac{1}{2}(u^i_{k1})^2 & \dots &\frac{1}{2}(u^i_{kt})^2 & \frac{\sqrt{2}}{2}u^i_{k1}u^i_{k2}& \frac{\sqrt{2}}{2}u^i_{k1}u^i_{k3} & \dots & \frac{\sqrt{2}}{2}u^i_{k,t-1}u^i_{kt} \end{bmatrix} \] If we let \(\vec{f}_i = \left(f(\vec{x}_{i_1}),\dots,f(\vec{x}_{i_k})\right)^T\) be a vector containing the function values at the nearest neighbors of \(\vec{x}_i\), then we can estimate the coefficients of the Taylor expansion using the following regression formula \[\vec{f}_i = {\bf X}_i \vec{\beta}\] where the last \(t + t(t+1)/2\) terms of \(\vec{\beta}_i\) correspond to the entries of the \(H^{tan}(f)\) at \(\vec{x}_i\). We then estimate \(\vec{\beta}_i\) using the Moore-Penrose pseudoinverse \[\vec{\beta}_i \approx ({\bf X}_i^T{\bf X}_i)^{-1}{\bf X}_i^T \vec{f}_i.\] For the pseudoinverse to be invertible, \({\bf X}_i\) must have at least as many rows as columns. Thus, we need \(k \ge 1+t + t(t+1)/2 = 1 + t(t+3)/2 > t(t+3)/2\).

We can drop the first \(t+1\) entries which contain the intercept and the first order regression terms corresponding to the gradient. Let \(\vec{\beta}_{i,drop}\) denote the final \(t(t+1)/2\) entries of \(\vec{\beta}\). If we square then sum the entries of \(\vec{\beta}_{i,drop}\), we now have an estimate of \(\| H^{tan}(f)\|_F^2\) at \(\vec{x}_i\).

One final note on this issue. The \(1/2\) and \({\sqrt{2}}/2\) terms in \({\bf X}_i\) are introduced to ensure we are correctly estimating (and single counting) the diagonal entries of the Hessian while double counting its upper triangular elements. This detail is missing from the original paper [32]. As such, many implementations may also replicate the mistake. Be cautious when choosing a package for implementing HLLEs!

5.7.1.3 The Eigenvalue Problem

Define \({\bf S}_i \in \mathbb{R}^{k \times N}\) as follows \[({\bf S}_i)_{j\ell} = \begin{cases} 1 & \vec{x}_\ell \text{ is the } j \text{th nearest neighbor of } \vec{x}_i \\ 0 & \text{else.} \end{cases}\] Using this notation, we can express \(\vec{f}_i\) (the function values from the k nearest neighbors of \(\vec{x}_i\)) as the product of \({\bf S}_i\) and \(\vec{f}\) (the vector containing the function value at all samples \(\vec{x}_1,\dots,\vec{x}_N\). Specifically, \(\vec{f}_i= {\bf S}_i \vec{f}\).

Now, we will use an empirical average of our Hessian estimates at each data point to approximate the operator \(\mathcal{H}(f)\) as follows \[\begin{align*} \mathcal{H}(f) &= \int_{\mathcal{M}} \|H^{tan}(f)\|_F^2 dm \\ &\approx \frac{1}{N}\sum_{i=1}^N \|\|H^{tan}\big(f\big)(\vec{x}_i)\|_F^2 \\ &\approx \frac{1}{N}\sum_{i=1}^N \|{\bf H}_i \vec{f}_i\|^2 \end{align*}\]

Now, we can make use of the identity \(\vec{f}_i={\bf S}_i \vec{f}\) to write this approximation as a quadratic function of \(\vec{f}.\) \[\begin{align*} \approx \frac{1}{N}\sum_{i=1}^N \|{\bf H}_i \vec{f}_i\|^2 &= \frac{1}{N}\sum_{i=1}^N \vec{f_i}^T{\bf H}_i^T{\bf H}_i \vec{f}_i\\ &= \frac{1}{N} \sum_{i=1}^N \vec{f}^T{\bf S}_i^T{\bf H}_i^T{\bf H}_i {\bf S}_i \vec{f} \\ &= \vec{f}^T\left(\frac{1}{N} \sum_{i=1}^N{\bf S}_i^T {\bf H}_i^T{\bf H}_i{\bf S}_i\right)\vec{f} \end{align*}\]

Letting \[{\bf H} = \frac{1}{N} \sum_{i=1}^N{\bf S}_i^T {\bf H}_i^T{\bf H}_i{\bf S}_i\] our approximation becomes \[\vec{f}^T{\bf H}\vec{f}.\] Here is the final critical observation. If \(f\) is a vector in the null space of \(\mathcal{H}\) then it must be in the span of the coordinate functions, and we would expect the approximation to be small (close to zero except for sampling variability and estimation error). As such, we can look eigenvector(s) associated with the smallest eigenvalue(s) of \({\bf H}\) to give approximations to the coordinate functions!

Like Laplacian Eigenmaps and LLEs, \({\bf H}\) will have one eigenvalue which is exactly zero corresponding to the constant function. We’ll take the next \(t\) eigenvectors associated with the next smallest \(t\) eigenvalues and use them to build a data matrix. The rows of this data matrix give the corresponding low dimensional coordinates of our data.

Since \({\bf H}\) is symmetric and positive semidefinite (it is the sum of symmetric, psd matrices), its eigenvectors will be orthogonal. Selecting these eigenvectors as approximates for functions \(f_1,\dots,f_t\) is equivalent to imposing an orthogonality constraint on these vectors. We’ll also take the vectors to be unit length. As a result, we can expect HLLE to recover original coordinates up to rigid motion (scaling and rotation/reflection) and coordinate rescaling.

5.7.2 Assumptions, Strength, and Weaknesses

HLLEs rely on two assumptions: the manifold map is locally isometric and the lower-dimensional set \(A\) is connected. These are far more flexible than ISOMAP, which requires an isometric manifold map and convex \(A\). Additionally, HLLEs have fewer tuning parameters than Laplacian eigenmaps, which also require one to specify the bandwidth \(\sigma^2\) is the rbf kernel.

However, there are also drawbacks. The theoretical justification for HLLEs relies on a quality approximation to the local Hessian, which is difficult without a large number of samples. Furthermore, for a manifold with intrinsically high dimension \(t\), a prohibitively large number of nearest neighbors is required.

5.8 Uniform Manifold Approximation and Projection (UMAP)

5.8.1 Introduction

UMAP was introduced by McInnes, Healy, and Melville in 2018 [33] and has quickly become one of the most widely used tools for visualizing high-dimensional data. Like Laplacian eigenmaps, UMAP is a local method built from a \(k\)-nearest neighbor graph rather than the pairwise (geodesic) distances used by ISOMAP. Unlike every method we have discussed so far, however, UMAP does not produce its embedding by solving a single eigenvalue problem. Instead, it builds a weighted graph representing the neighborhood structure of the data in \(\mathbb{R}^d\), builds a comparable weighted graph on a set of candidate points in \(\mathbb{R}^t\), and then directly optimizes the positions of those points so that the two graphs are as similar as possible. This shift from an eigenvalue problem to a direct optimization problem is the central idea that distinguishes UMAP from ISOMAP, LLE, Laplacian eigenmaps, and HLLEs.

5.8.2 Algorithm

5.8.2.1 Constructing a fuzzy neighborhood graph

As with Laplacian eigenmaps, we begin by fixing a number of nearest neighbors \(k>0\) and, for each \(\vec{x}_i\), finding its \(k\)-nearest neighbors \(\mathcal{N}_i\). Rather than a single bandwidth \(\sigma^2\) shared across the entire dataset, UMAP calibrates a local bandwidth \(\sigma_i\) for each point \(\vec{x}_i\), allowing the method to adapt to regions of the manifold with different sampling density. Let \(\rho_i = \min_{j\in\mathcal{N}_i} \|\vec{x}_i-\vec{x}_j\|\) be the distance from \(\vec{x}_i\) to its nearest neighbor. We then choose \(\sigma_i\) so that \[\sum_{j\in\mathcal{N}_i} \exp\left(-\frac{\max(0,\|\vec{x}_i-\vec{x}_j\|-\rho_i)}{\sigma_i}\right) = \log_2(k),\] which can be solved numerically for each \(i\). Subtracting \(\rho_i\) guarantees that \(\vec{x}_i\) is connected to its closest neighbor with full weight, while the choice of \(\sigma_i\) ensures that every point contributes the same total weight to the graph regardless of how densely sampled its neighborhood happens to be. These choices give a directed, weighted graph with entries \[\tilde{\bf W}_{ij} = \begin{cases} \exp\left(-\dfrac{\max(0,\|\vec{x}_i-\vec{x}_j\|-\rho_i)}{\sigma_i}\right) & \vec{x}_j \in \mathcal{N}_i \\ 0 & \text{else.}\end{cases}\]

5.8.2.2 Symmetrization

The matrix \(\tilde{\bf W}\) need not be symmetric, since \(\vec{x}_j\) being a neighbor of \(\vec{x}_i\) does not guarantee \(\vec{x}_i\) is a neighbor of \(\vec{x}_j\). UMAP symmetrizes \(\tilde{\bf W}\) using the fuzzy set union \[{\bf W}_{ij} = \tilde{\bf W}_{ij} + \tilde{\bf W}_{ji} - \tilde{\bf W}_{ij}\tilde{\bf W}_{ji},\] so that \(\vec{x}_i\) and \(\vec{x}_j\) receive a high weight \({\bf W}_{ij}\) if either point considers the other a close neighbor. The resulting \({\bf W}\) plays the same role as the adjacency matrix \({\bf W}\) from Laplacian eigenmaps, just calibrated with this locally adaptive procedure.

5.8.2.3 Optimizing the low-dimensional embedding

Given the high-dimensional affinities \({\bf W}\), UMAP next defines an analogous affinity between candidate low-dimensional points \(\vec{y}_1,\dots,\vec{y}_N\in\mathbb{R}^t\), \[v_{ij} = \left(1+a\|\vec{y}_i-\vec{y}_j\|^{2b}\right)^{-1},\] where \(a,b>0\) are constants chosen so that \(v_{ij}\) closely matches a target curve controlled by the tuning parameter min_dist, which specifies how tightly points are allowed to pack together in the embedding. Initializing \(\vec{y}_1,\dots,\vec{y}_N\) using the spectral embedding from Laplacian eigenmaps (a convenient reuse of a method we have already discussed), UMAP then adjusts the points to minimize the cross-entropy between \({\bf W}\) and \({\bf V}\), \[C(\vec{y}_1,\dots,\vec{y}_N) = \sum_{i \ne j} \left[{\bf W}_{ij} \log\left(\frac{{\bf W}_{ij}}{v_{ij}}\right) + (1-{\bf W}_{ij})\log\left(\frac{1-{\bf W}_{ij}}{1-v_{ij}}\right)\right].\] Minimizing this loss via stochastic gradient descent produces an attractive force which pulls \(\vec{y}_i\) and \(\vec{y}_j\) together whenever \({\bf W}_{ij}\) is large, and a repulsive force which pushes \(\vec{y}_i\) and \(\vec{y}_j\) apart otherwise. Since most entries of \({\bf W}\) are (close to) zero, evaluating the repulsive force for every pair at every iteration would be wasteful; in practice, UMAP approximates it using a small number of randomly sampled non-neighbor pairs at each step, a trick known as negative sampling. The final coordinates \(\vec{y}_1,\dots,\vec{y}_N\) after optimization are the UMAP embedding.

Example 5.8 (UMAP applied to the Swiss Roll) Let’s revisit the swiss roll from earlier in this chapter and see how UMAP recovers its two-dimensional structure.

Figure 5.3: The Swiss Roll (left) and its UMAP embedding (right)

UMAP successfully unrolls the swiss roll, and the coloring shows that nearby points on the manifold remain nearby in the embedding. Notice, however, that the embedding is not a faithful rectangle in the way that ISOMAP’s embedding is – distances within the embedding, and the overall aspect ratio, need not match the original data closely. UMAP prioritizes preserving neighborhoods, not distances.

5.8.3 Strengths and Weaknesses

5.8.3.1 Strengths

1. Speed and scalability

Because UMAP relies on approximate nearest neighbor search and stochastic gradient descent with negative sampling rather than the eigendecomposition of a dense \(N\times N\) matrix, it scales far more gracefully to large datasets than ISOMAP, Laplacian eigenmaps, or HLLEs.

2. Balances local and global structure

The spectral initialization (itself a Laplacian eigenmap of \({\bf W}\)) gives the optimization a sensible global starting point before the attractive and repulsive forces refine local neighborhoods, so UMAP tends to preserve more of the manifold’s coarse, global organization than LLE or Laplacian eigenmaps alone.

3. Few tuning parameters

The two main tuning parameters, the number of neighbors \(k\) and min_dist, play a similar role to the tuning parameters of LLE and Laplacian eigenmaps and have reasonably interpretable effects on the resulting embedding.

5.8.3.2 Weaknesses

1. No uniqueness guarantee

Because the embedding is the result of a nonconvex optimization rather than an eigenvalue problem, different random initializations or optimization runs can produce noticeably different embeddings, unlike the (essentially) unique solutions returned by ISOMAP, LLE, Laplacian eigenmaps, and HLLEs.

2. Distances between (and sizes of) clusters are not meaningful

UMAP is explicitly designed to preserve neighborhoods, not distances. As a result, the relative sizes of clusters and the gaps between them in a UMAP embedding often do not reflect the true geometry of the data and should not be over-interpreted.

3. No natural inverse mapping

Like ISOMAP, LLE, and Laplacian eigenmaps, UMAP does not directly provide a mapping from the embedding back to the original space. We turn to a class of methods, autoencoders, which are built explicitly to provide such an inverse.

5.9 Autoencoders (AEs)

5.9.1 Introduction

In each of the previously discussed manifold learning techniques, we have never made an attempt to explicity model the manifold mapping (or its inverse). Rather, we have focused on methods which recover lower-dimensional representations of the data assuming certain properties of the manifold map (e.g. isometry) and its preimage (convexity). As a result, the journey from high-dimensional data to lower-dimensional representations has been a one-way trip.

Suppose now, that we take a more ambitious approach and try to model the manifold mapping and its inverse. In practice, we know nothing about these functions – their mere existence is assumed via the manifold hypothesis – but at minimum we can expect them to be complicated nonlinear functions. Enter autoencoders (AEs) which are specially designed neural networks that simultaneously model the manifold map and its inverse.

Generally, an autoencoder is comprised of three main components.

1. Encoder: A neural network \(\mathcal{E}(\vec{x})\) that compresses the input \(\vec{x}\) into a latent representation. The encoder is acting like the inverse manifold map.

2. Latent Space: The reduced dimensionality representation, often denoted as \(\vec{z}\), where \(\vec{z} = \mathcal{E}(\vec{x})\).

3. Decoder: A neural network \(\mathcal{D}(\vec{z})\) that aims to reconstruct the original input from the latent representation. This function acts like the manifold map.

Figure 5.4: Informal schematic of an autoencoder. The encoder and decoder are neural networks which map the original data into a lower dimensional latent space.

5.9.2 Neural networks: design and training

In their simplest form, neural networks iterate between affine transformations and elementwise nonlinear operations. If \(\vec{x}\in\mathbb{R}^d\) is the input to a neural network, the affine transformation is determined by the equation \({\bf W}\vec{x}+\vec{b}\) where \({\bf W}\in\mathbb{R}^{k_1\times d}\), called weights, and \({\bf b}\in\mathbb{R}^{k_1}\), called the bias, are tuneable parameters. We then apply a nonlinear function, \(\sigma_1:\mathbb{R}\to \mathbb{R}\), to \({\bf W}\in\mathbb{R}^{k_1\times d}\) resulting in a \(k_1\) dimensional vector, \(\sigma_1\left({\bf W}\in\mathbb{R}^{k_1\times d}\right)\) where the application of \(\sigma_1\) is interpreted elementwise. The nonlinear function is also called an activation. Together, the affine transformation and nonlinear transformation are called a layer of the network with width \(k_1\) reflecting the dimension of the output of the layer.

We can then iterate this idea, using the output of the previous layer as the input to the subsequent layers resulting in the sequence of operations, \[\begin{align*} \vec{\ell}_0 &= \vec{x} \\ \vec{\ell}_1 &=\sigma_1({\bf W}_1\vec{\ell}_0+\vec{b}_1) \\ \vec{\ell}_1 &=\sigma_2({\bf W}_2\vec{\ell}_1+\vec{b}_2) = \sigma_2\left({\bf W}_2\sigma_1({\bf W}_1\vec{\ell}_0+\vec{b}_1)+\vec{b}_2\right) \\ \vdots \\ \vec{\ell}_t &=\sigma_t({\bf W}_t\vec{\ell}_{t-1}+\vec{b}_t). \end{align*}\] In practice, we are free to choose the number of layers (depth of the network) and the width and activation within each layer providing innumerable configurations. Choosing the structure of the network – its architecture – has an enormous impact on how well this model can perform as a method of dimension reduction. Given the importance of these choices, building a neural network can appear overwhelming.

And while there is still some art to the practice there are several conventions specific to autoencoders which offer guidance. First, we want to compress our high dimensional data into a lower dimensional space. To this end, we’ll have an internal layer which maps to a lower dimensional space, referred to as the bottleneck. We refer to all layers between the input and the bottleneck as the encoder. The layers following the bottleneck then map back to the original high-dimensional space. We call this portion of the network the decoder. By convention, the network is symmetric around the bottleneck so the encoder and decoder have the same number of layers with widths in the opposite order. See the MNIST example below for a specific case.

There are also other data specific insights from geometric deep learning that motivate more complicated neural network architecture beyond the traditional, dense, feed-forward approach we have outlined here. For now, we’ll focus common activations and methods for training the model (tuning the weights and biases in layers) to accomplish dimension reduction.

5.9.2.1 Activation Functions

Activation functions play a vital role in neural networks, determining the output of a neuron based on its input. They introduce non-linear properties into the network, enabling it to learn from the error and make adjustments, which is essential for learning complex, nonlinear patterns.

Two common activation functions we can choose from.

1. Sigmoid (Logistic) Function:

Equation: \(f(x) = \frac{1}{1 + e^{-x}}\)

Range: Between 0 and 1

Pros: Smooth gradient, preventing “jumps” in output values

Cons: Can cause vanishing gradient problems in deep networks

2. Rectified Linear Unit (ReLU) Function:

Equation: \(f(x) = max(0, x)\)

Range: From 0 to infinity

Pros: Helps mitigate the vanishing gradient problem, leading to faster convergence

Cons: Neurons can sometimes “die”, especially with a large learning rate

There are many other choices, but the ReLu activation is a fantastic first choice since it results in a network which can be trained more easily.

5.9.2.2 Training

At its core, an autoencoder seeks to learn a compressed, efficient representation of input data by nonlinearly encoding the data into a lower dimension space and subsequently decoding it to reconstruct the original data [34]. This translates the training of a neural network into a self-supervised problem, wherein we want to approximate \(\vec{x}\) by its reconstruction \(\mathcal{D}\left(\mathcal{E}(\vec{x})\right).\) The primary objective during the training phase of an autoencoder is to minimize the reconstruction error, often quantified using metrics such as Mean Squared Error (MSE) between the input data and its reconstructed counterpart. Mathematically, this gives rise to the loss function \[L(\theta_e,\theta_d) = \frac{1}{N}\sum_{i=1}^N\|\vec{x}_i - \mathcal{D}_{\theta_d}\left( \mathcal{E}_{\theta_e}(\vec{x_i})\right)\|^2.\] Here \(\theta_e\) and \(\theta_d\) are the weights and biases of the encoder and decoder respectively.

Importantly, since the network is comprised of iterated affine transformations and nonlinear operations, it is a differentiable function of its parameters. Computing these derivatives is called backpropagation in the ML literature and relies on repeated applications of the chain rule. With gradients in hand, gradient based optimization is a natural choice for tuning the parameters of the model to minimize the MSE. At present, the most common choice of gradient based optimization is an algorithm called Adam [35] which uses an accelerated version of stochastic gradient descent to fit the autoencoder. The details of Adam and stochastic gradient descent are left to the curious reader.

Once the autoencoder has been trained and we have identified optimal parameters, we can use the encoder for dimension reduction. As a motivating example, let’s revisit MNIST.

Example 5.9 (AEs and MNIST) Compared to other manifold learning methods we have introduced, the implementation and training of an Autoencoder is much more complicated as we’ve outlined above. Here, we focus again on MNIST. Following the setup of other methods, we continue to flatten each \(28\times 28\) image into a \(784\)-dimensional vector, but we take one more step. Pixel intensities take values in the interval [0,255], which we divide by 255 so that every entry is between zero and one. We then apply the following dense neural network, which we train on the full MNIST data set. This network has 10 layers (we don’t count the input layer). All layers use a ReLu activation except for the last layer of the encoder (so the bottleneck is in \([0,1]^2\)) and the output (so that pixel intensities in the reconstruction are bounded between [0,1] matching the division by 255 done in preprocessing)

Figure 5.5: A deep, dense neural network for modeling MNIST data

Between the encoder and decoder, this model has 1,892,954 parameters, which we train via the Adam optimizer [35] using 200 epochs, batch size = 256, and a 80%/20% test/train split. After training, we can run each image through the encoder portion of the network to view where each image is mapped in the latent space. As we can see from the image below, the neural network does a far better job of separating the images of different digits into separate groups.

Figure 5.6: Latent space representation of MNIST

For comparison, here is the two-dimensional configuration of points obtained by PCA.

Figure 5.7: 2D Representation of MNIST via PCA

5.9.3 Strengths and weaknesses

Neural networks are incredibly powerful tools for approximating nonlinear functions, so it is natural to hope that we might be able to approximate the manifold map using the decoder. Like previous methods, however, the lower dimensional space is not unique so we can hope to learn the manifold map up to some invertible transformation without a lot more information. Furthermore, AEs are focused on minimizing the reconstruction error rather than directly using properties of the manifold containing our data. The greedy gradient based approach means we relinquish control on which relationships within the dataset we would like to preserve and what other relationships that might exist between the bottleneck and manifold (isometry for example). In particular, there are no constraints that an AE prioritize the neighborhood structure between points. Samples \(\vec{x}_i\) and \(\vec{x}_j\) may be close originally, but the nonlinear nature of a neural network could result in \(\mathcal{E}(\vec{x}_i)\) and \(\mathcal{E}(\vec{x}_j)\) being far apart.

Nonetheless, the utility of Autoencoders has been demonstrated in a wide array of applications beyond dimensionality reduction. Anomaly detection, denoising [36], and generative modeling [37] are a few exciting extensions and applications. The capacity of an AE to generalize and generate new samples is by far its greatest strength, which delineates it from the preceding methods. LLMs and many other generative models are based on the self-supervised approach central to AEs. The specifics that go into training state of the art language models are beyond the scope of this text. However, a basic mechanism is built into the decoder portion of the network. In principle, we could choose a point in the bottleneck (perhaps one between the lower-dimensional representation of two sample arising from the encoder), then pass it through the decoder to generate a new data point. For an interactive demonstration of this idea applied to MNIST, check out this java applet hosted by n8programs. The autoencoder architecture in that demo is different than the method we used in the example above hence the different organization of points in the bottleneck.

5.10 Comparison of Manifold Learning Methods

As an initial comparison, we’ll look at several 3-dimensional examples with known lower dimensional latent spaces where we can visualize the data. Importantly, these examples fail to satisfy some of the conditions of the preceding methods so we can assess how susceptible the methods we have discussed to violations in their assumptions. Table 5.1 summarizes the results.

Table 5.1: Manifold learning examples and their properties

Manifold	Intrinsic Dimension	Preimage	Preimage connected?	Preimage convex?	Isometric Manifold Map?	Locally Isometric Manifold Map?
Helix	One	Line segment	Yes	Yes	Yes	Yes
Swiss Roll	Two	Rectangle	Yes	Yes	No	No
Folded Washer	Two	Annulus	Yes	No	No	No
Rolled Washer	Two	Annulus	Yes	No	No	No

For comparison, we provide scatterplots of 2000 samples from each manifold which we will use for analysis. The points have been color coded to improve visualization and to aid comparison with recovered lower dimensional coordinates.

5.10.1 Visual comparison of results

5.10.1.1 Helix

The helix is the simplest example and satisfies all assumptions we have discussed in the preceding sections. As such, it is not surprising that each of the methods discussed does a reasonable job compressing the data to a line segment. For the helix, we have used the isometric manifold mapping \[\Psi(t) = \left(\frac{\sqrt{2}}{2}\cos (t), \frac{\sqrt{2}}{2}\sin t, \frac{\sqrt{2}}{2} t\right)^T\] where \(t \in [0,15]\) is the original 1-dimensional coordinate for the data. Below we plot the recovered 1-d representation of the data against the corresponding value of \(t\). As we have discussed, a perfect recovery of the original \(t\) is impossible. The best we can hope for is an affine relationship reflecting a one-to-one correspondence (up to translation, rescaling, and rigid motion) between the original low dimension coordinates and those recovered by our manifold learning methods.

Figure 5.8: Recovered latent space vs original coordinates

5.10.1.2 Swiss Roll

For the Swiss roll, we use the manifold map \[\Psi: \left[\frac{3\pi}{2},\frac{9\pi}{2}\right] \times [0,15] \to \mathbb{R}^3\] given by the equation \[\Psi(s,t) = \left(s\cos s, t, s\sin s\right)^T.\] Below, we compare the two-dimensional configurations recovered by each of the methods discussed in this chapter.

Figure 5.9: Recovered latent space vs original coordinates

5.10.1.3 Folded Washer

Figure 5.10: Recovered latent space vs original coordinates

5.10.1.4 Rolled Washer

Figure 5.11: Recovered latent space vs original coordinates

5.11 A Unifying View: Manifold Learning as Kernel PCA

Having covered ISOMAP, LLE, and Laplacian eigenmaps, it is worth stepping back to notice a common thread. Each method feels different – one preserves geodesic distances, one preserves local linear reconstructions, one preserves a weighted graph structure – yet each method ultimately produces its embedding the same way: by eigendecomposing a particular symmetric matrix built from the data and taking the leading eigenvectors (suitably scaled) as coordinates. This is exactly the recipe used by kernel PCA. It turns out that this similarity is not a coincidence. ISOMAP, LLE, and Laplacian eigenmaps can each be written as kernel PCA applied to a specific, data-dependent “kernel” matrix in place of the doubly centered kernel matrix \({\bf HKH}\) from earlier in this chapter [38].

ISOMAP. Recall that ISOMAP applies classical scaling directly to the matrix of squared geodesic distances, forming \(B = -\frac{1}{2}{\bf H}\Delta^{(2)}{\bf H}\) where \(\Delta^{(2)}_{ij} = \left(d^{\mathcal{G}}_{ij}\right)^2\). By the duality of PCA and classical scaling (4.4), this is precisely kernel PCA using the kernel matrix \({\bf K}_{\text{ISOMAP}} = B\). Since geodesic distances are only estimated from a finite neighborhood graph, \({\bf K}_{\text{ISOMAP}}\) is not guaranteed to be positive semidefinite in practice (unlike a kernel matrix built from a valid Mercer kernel); small negative eigenvalues are typically discarded.

LLE. LLE embeds the data using the smallest nontrivial eigenvectors of \({\bf M} = ({\bf I}-{\bf W})^T({\bf I}-{\bf W})\), whereas kernel PCA always uses the largest eigenvectors of a kernel matrix. The two views are reconciled by a simple eigenvalue flip. With \(\lambda_{\max}\) the largest eigenvalue of \({\bf M}\), the top nontrivial eigenvectors of \({\bf K}_{\text{LLE}} = \lambda_{\max}{\bf I} - {\bf M}\) coincide with the bottom nontrivial eigenvectors of \({\bf M}\), so the LLE embedding is exactly the kernel PCA embedding using \({\bf K}_{\text{LLE}}\).

Laplacian eigenmaps. By the same flip argument, the Laplacian eigenmap embedding built from the smallest nontrivial eigenvectors of \({\bf L}_{sym}\) is the kernel PCA embedding using the (Moore-Penrose) pseudoinverse \({\bf K}_{\text{LE}} = {\bf L}_{sym}^{\dagger}\) as the kernel matrix.

Two consequences follow from this unifying view. First, it explains why these methods behave so differently from kernel PCA with a fixed kernel like the rbf kernel. The matrices \({\bf K}_{\text{ISOMAP}}\), \({\bf K}_{\text{LLE}}\), and \({\bf K}_{\text{LE}}\) are all built from the entire training sample (through a \(k\)-nearest-neighbor graph), rather than from a function \(k(\vec{x},\vec{y})\) that can be evaluated on two arbitrary points. This is precisely why ordinary kernel PCA readily extends to new, out-of-sample points while ISOMAP, LLE, and Laplacian eigenmaps do not – doing so requires a Nyström-type approximation of the relevant kernel [38]. Second, it foreshadows a connection we revisit in the next chapter. The graph Laplacian \({\bf L}_{sym}\) used above to build \({\bf K}_{\text{LE}}\) is the same graph Laplacian at the heart of spectral clustering (6.5), and kernel \(k\)-means with a graph-based kernel recovers the same spectral clustering objective. Dimension reduction and clustering, at least for these methods, are two sides of the same eigenvalue problem.

5.12 Exercises

1. Let \(k_1\) and \(k_2\) be valid kernels on \(\mathbb{R}^d\times\mathbb{R}^d\) (i.e. each satisfies Mercer’s Condition 5.1).

   1. Show that \(k(\vec{x},\vec{y}) = k_1(\vec{x},\vec{y}) + k_2(\vec{x},\vec{y})\) is also a valid kernel. Hint: use the fact that the sum of two positive semidefinite matrices is positive semidefinite.

   2. For three points \(\vec{x}_1,\vec{x}_2,\vec{x}_3\in\mathbb{R}^d\), suppose a candidate kernel gives the matrix \[{\bf K} = \begin{bmatrix} 1 & 0.9 & 0.1 \\ 0.9 & 1 & 0.9 \\ 0.1 & 0.9 & 1\end{bmatrix}.\] Is \({\bf K}\) a valid kernel matrix? Justify your answer using Mercer’s Condition.

2. Consider the polynomial kernel \(k(\vec{x},\vec{y}) = (1+\vec{x}^T\vec{y})^2\) for \(\vec{x},\vec{y}\in\mathbb{R}^2\). Find an explicit feature map \(\varphi:\mathbb{R}^2\to\mathbb{R}^m\) (for some finite \(m\)) such that \(k(\vec{x},\vec{y}) = \varphi(\vec{x})^T\varphi(\vec{y})\).

3. Recall that kernel PCA computes scores from the eigendecomposition of the doubly centered kernel matrix \({\bf HKH}\) where \({\bf K}_{ij} = k(\vec{x}_i,\vec{x}_j)\). Suppose we use the linear kernel \(k(\vec{x},\vec{y}) = \vec{x}^T\vec{y}\).

   1. Show that \({\bf K} = {\bf XX}^T\) where \({\bf X}\) is the (uncentered) data matrix with rows \(\vec{x}_1^T,\dots,\vec{x}_N^T\).

   2. Using the duality of PCA and classical scaling, explain why kernel PCA with the linear kernel recovers exactly the same scores as ordinary PCA on \(\vec{x}_1,\dots,\vec{x}_N\).

4. Let \(\vec{x}_1 = -1,\, \vec{x}_2 = 0,\, \vec{x}_3 = 1 \in \mathbb{R}\), and consider the rbf kernel \(k(x,y) = \exp\left(-(x-y)^2\right)\) (i.e. \(\sigma = 1\)).

   1. Compute the \(3\times 3\) kernel matrix \({\bf K}\) by hand.

   2. Using R, compute the doubly centered matrix \({\bf HKH}\) and its eigendecomposition. Report the first kernel PC scores for \(\vec{x}_1,\vec{x}_2,\vec{x}_3\).

5. For each of the following subsets of \(\mathbb{R}^2\), state whether it is a manifold. If it is, give its intrinsic dimension; if it is not, explain which point(s) fail Definition 5.1 and why.

   1. Two lines through the origin with different slopes, i.e. \(\{(t,t): t\in\mathbb{R}\} \cup \{(t,-t):t\in\mathbb{R}\}\).

   2. The open unit disk \(\{\vec{x}\in\mathbb{R}^2: \|\vec{x}\| < 1\}\).

   3. A “T”-shaped union of two segments, \(\{(0,y): -1\le y \le 1\} \cup \{(x,0): -1 \le x \le 1\}\).

   4. The unit circle \(\{\vec{x}\in\mathbb{R}^2: \|\vec{x}\|=1\}\).

6. Let \(A = \mathbb{R}\) and define the manifold map \(\Psi:A\to\mathbb{R}^3\) by \(\Psi(z) = (z,z^2,z^3)^T\).

   1. Compute the Jacobian \({\bf J}_\Psi\) as a function of \(z\).

   2. Find the tangent space \(T_{\vec{p}}(\mathcal{M})\) at the point \(\vec{p} = \Psi(1) = (1,1,1)^T\).

   3. Find the approximating tangent plane to \(\mathcal{M}\) at \(\vec{p}\).

7. Consider the map \(\Psi(t) = (\cos t,\sin t)^T\) for \(t\in(0,2\pi)\), which traces out (most of) the unit circle.

   1. Show that \(\|d\Psi/dt\| = 1\) for all \(t\), and conclude that \(\Psi\) is a local isometry, i.e. for \(t_1,t_2\) close together, the manifold distance between \(\Psi(t_1)\) and \(\Psi(t_2)\) equals \(|t_1-t_2|\).

   2. Show that \(\Psi\) is not a global isometry by comparing the manifold distance between \(\Psi(0.1)\) and \(\Psi(2\pi-0.1)\) (measured by traveling around the circle the short way) to \(|0.1 - (2\pi - 0.1)|\).

8. Consider four points connected in a path, with edge weights \(d_{12}=1,\, d_{23}=2,\, d_{34}=1\) (there is no direct edge between any other pair of points).

   1. Write down the initial \(4\times 4\) distance matrix \({\bf \Delta}^{(0)}\) used to initialize the Floyd-Warshall algorithm (using \(\infty\) for pairs with no direct edge).

   2. Run the Floyd-Warshall algorithm by hand to find the geodesic distance matrix \(\Delta\) with entries \(\Delta_{ij} = d^{\mathcal{G}}_{ij}\).

9. Recall from Table 5.1 that the preimage of the folded washer is an annulus, which is not convex. Explain, in terms of the neighborhood graph construction in ISOMAP’s first step, why a non-convex preimage can cause the graph distance \(d^{\mathcal{G}}_{ij}\) to badly overestimate the true manifold distance \(d^{\mathcal{M}}_{ij}\) for some pairs of points.

10. Suppose \(\vec{x}_i\in\mathbb{R}\) has exactly two neighbors, \(\vec{x}_j\) and \(\vec{x}_l\), so that \(N_i^k = \{j,l\}\).

    1. Using the constraint \(w_{ij}+w_{il}=1\), show that the LLE reconstruction problem \(\min (\vec{x}_i - w_{ij}\vec{x}_j - w_{il}\vec{x}_l)^2\) has the closed-form solution \[w_{ij} = \frac{\vec{x}_i - \vec{x}_l}{\vec{x}_j-\vec{x}_l}, \qquad w_{il} = \frac{\vec{x}_j - \vec{x}_i}{\vec{x}_j-\vec{x}_l}.\]

    2. For \(\vec{x}_i = 2,\, \vec{x}_j = 0,\,\vec{x}_l = 5\), compute \(w_{ij}\) and \(w_{il}\) numerically and verify that \(\vec{x}_i = w_{ij}\vec{x}_j + w_{il}\vec{x}_l\).

11. Let \(\vec{x}_1,\dots,\vec{x}_N\in\mathbb{R}^D\) have optimal LLE reconstruction weight matrix \({\bf W}\). Show directly (i.e. without invoking the general argument for \(\epsilon(w)\) in the text) that if every point is reflected through the origin, i.e. \(\vec{x}_i \mapsto -\vec{x}_i\) for all \(i\), the optimal reconstruction weights are unchanged.

12. Let \(\vec{x}_1 = 0,\,\vec{x}_2=1,\,\vec{x}_3=5,\,\vec{x}_4=6 \in\mathbb{R}\). Use \(\epsilon = 2\) and binary weights to build the neighborhood graph for Laplacian eigenmaps.

    1. Write down the adjacency matrix \({\bf W}\), the degree matrix \({\bf D}\), and the graph Laplacian \({\bf L} = {\bf D}-{\bf W}\).

    2. How many eigenvalues of \({\bf L}\) are exactly zero? Relate your answer to the number of connected components of the graph.

    3. For \(\vec{y} = (1,2,3,4)^T\), verify the identity \(\vec{y}^T{\bf L}\vec{y} = \frac{1}{2}\sum_{i=1}^4\sum_{j=1}^4 {\bf W}_{ij}(y_i-y_j)^2\) ((5.3)) by computing both sides directly.

13. In HLLEs, the local Hessian is estimated using quadratic regression with design matrix \({\bf X}_i \in \mathbb{R}^{k\times(1+t+t(t+1)/2)}\), and we require \(k\) to be at least the number of columns of \({\bf X}_i\) for the pseudoinverse to be well defined.

    1. Explain why \({\bf X}_i\) has \(1+t+t(t+1)/2\) columns by counting the intercept, linear, and quadratic/cross terms in the local quadratic Taylor expansion.

    2. Show that the requirement \(k \ge 1+t+t(t+1)/2\) is equivalent to \(k > t(t+3)/2\).

    3. For an intrinsic dimension of \(t=5\), what is the minimum number of neighbors \(k\) required?

14. Let \(g(x,y) = xy\).

    1. Compute \(\nabla g\) and \(\mathcal{H}g\).

    2. Compute \(\|\mathcal{H}g\|_F^2\). Is it constant over \(\mathbb{R}^2\)?

    3. Compare your answer in (b) to \(\|\mathcal{H}g_1\|_F^2 = 8\) from the text’s example with \(g_1(x,y)=x^2+y^2\). Which function is “more curved” on average over \([-2,2]^2\)?

15. Laplacian eigenmaps use a single bandwidth \(\sigma^2\) shared by every point, while UMAP calibrates a local bandwidth \(\sigma_i\) for each point \(\vec{x}_i\) individually. Suppose a dataset has one very densely sampled region and one very sparsely sampled region.

    1. Explain what is likely to go wrong if a single, shared \(\sigma^2\) is tuned to work well in the dense region and then applied to the sparse region (or vice versa).

    2. Explain how UMAP’s local bandwidth calibration, via the condition \(\sum_{j\in\mathcal{N}_i} \exp\left(-\max(0,\|\vec{x}_i-\vec{x}_j\|-\rho_i)/\sigma_i\right) = \log_2(k)\), addresses this issue.

16. Consider an autoencoder for data \(\vec{x}\in\mathbb{R}^{500}\) with encoder layer widths \(500 \to 100 \to 20 \to 3\) and a symmetric decoder \(3\to 20 \to 100 \to 500\) (each arrow indicates one layer with weights and a bias). Using the formula for the number of parameters in an affine layer, compute the total number of trainable parameters in the encoder, the decoder, and the full autoencoder.

17. Suppose the encoder and decoder of an autoencoder are both single linear layers with no bias and no activation function, i.e. \(\mathcal{E}(\vec{x}) = {\bf W}\vec{x}\) for \({\bf W}\in\mathbb{R}^{t\times d}\) and \(\mathcal{D}(\vec{z}) = {\bf V}\vec{z}\) for \({\bf V}\in\mathbb{R}^{d\times t}\), so that the reconstruction of \(\vec{x}_i\) is \({\bf VW}\vec{x}_i\).

    1. Explain why the reconstruction loss \[L({\bf W},{\bf V}) = \frac{1}{N}\sum_{i=1}^N \|\vec{x}_i - {\bf VW}\vec{x}_i\|^2\] can be written as \(\frac{1}{N}\|{\bf X} - {\bf X}({\bf VW})^T\|_F^2\) for the data matrix \({\bf X}\) with rows \(\vec{x}_1^T,\dots,\vec{x}_N^T\).

    2. The matrix \({\bf VW}\in\mathbb{R}^{d\times d}\) has rank at most \(t\). Using the Eckart-Young-Mirsky Theorem 4.3, explain why the optimal choice of \({\bf VW}\) is the best rank-\(t\) approximation to \({\bf X}\) (up to centering), and conclude that a linear autoencoder’s bottleneck spans the same subspace as the leading \(t\) PCA loadings.

18. Suppose \({\bf M} = ({\bf I}-{\bf W})^T({\bf I}-{\bf W})\) has eigenvalues \(0 = \mu_1 < \mu_2 = 1 < \mu_3 = 3 < \mu_4 = 4\) with corresponding eigenvectors \(\vec{v}_1,\vec{v}_2,\vec{v}_3,\vec{v}_4\).

    1. Using \(\lambda_{\max}=4\), compute the eigenvalues of \({\bf K}_{\text{LLE}} = \lambda_{\max}{\bf I}-{\bf M}\).

    2. Confirm that the ordering of the eigenvalues flips, i.e. \(\vec{v}_1\) (the trivial eigenvector of \({\bf M}\), associated with \(\mu_1=0\)) becomes the eigenvector associated with the largest eigenvalue of \({\bf K}_{\text{LLE}}\).

    3. Since \(\vec{v}_1\) corresponds to the trivial, constant-function direction, it must still be discarded even when working with \({\bf K}_{\text{LLE}}\). Which eigenvector should actually be used as the top coordinate of the LLE embedding, and to which eigenvalue of \({\bf K}_{\text{LLE}}\) does it correspond?

19. The S-curve is another standard manifold learning benchmark, defined by the map \(\Psi(s,t) = \left(\sin(s),\, t,\, \text{sign}(s)(\cos(s)-1)\right)^T\) for \(s\in(-3\pi/2,3\pi/2)\) and \(t\in(0,10)\).

    N <- 2000
    s <- runif(N, min = -3*pi/2, max = 3*pi/2)
    t <- runif(N, min = 0, max = 10)
    s_curve <- cbind(sin(s), t, sign(s)*(cos(s)-1))

    1. Generate \(N=2000\) samples from the S-curve using the code above, and produce a 3-dimensional scatterplot color-coded by \(s\).

    2. Apply ISOMAP, LLE, Laplacian eigenmaps, HLLEs, and UMAP to recover a two-dimensional embedding of these data.

    3. Which method(s) best recover the rectangular structure of the original \((s,t)\) coordinates? Which assumptions of each method are (or are not) satisfied by the S-curve, and how does this explain the differences you observe?

References

[22]

Mika, S., Schölkopf, B., Smola, A., Müller, K.-R., Scholz, M. and Rätsch, G. (1998). Kernel PCA and de-noising in feature spaces. In Advances in neural information processing systems vol 11, (M. Kearns, S. Solla and D. Cohn, ed). MIT Press.

[23]

Little, A., Lee, J., Jung, Y.-M. and Maggioni, M. (2009). Estimation of intrinsic dimensionality of samples from noisy low-dimensional manifolds in high dimensions with multiscale SVD. In IEEE Workshop on Statistical Signal Processing Proceedings pp 85–8.

[24]

Tenenbaum, J. B., Silva, V. de and Langford, J. C. (2000). A global geometric framework for nonlinear dimensionality reduction. Science 290 2319–23.

[25]

Bernstein, M., Silva, V., Langford, J. and Tenenbaum, J. (2001). Graph approximations to geodesics on embedded manifolds.

[26]

Roweis, S. T. and Saul, L. K. (2000). Nonlinear dimensionality reduction by locally linear embedding. Science 290 2323–6.

[27]

Chen, J. and Liu, Y. (2011). Locally linear embedding: A survey. Artif. Intell. Rev. 36 29–48.

[28]

Saul, L. K. and Roweis, S. T. (2003). Think globally, fit locally: Unsupervised learning of low dimensional manifolds. 4 119–55.

[29]

Anon. (2019). Locally linear embedding with additive noise. Pattern Recognition Letters 123 47–52.

[30]

Chang, H. and Yeung, D.-Y. (2006). Robust locally linear embedding. Pattern Recognition 39 1053–65.

[31]

Belkin, M. and Niyogi, P. (2001). Laplacian eigenmaps and spectral techniques for embedding and clustering. In Advances in neural information processing systems vol 14, (T. Dietterich, S. Becker and Z. Ghahramani, ed). MIT Press.

[32]

Donoho, D. L. and Grimes, C. (2003). Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data. Proceedings of the National Academy of Sciences 100 5591–6.

[33]

McInnes, L., Healy, J. and Melville, J. (2018). UMAP: Uniform manifold approximation and projection for dimension reduction. arXiv preprint arXiv:1802.03426.

[34]

Hinton, G. E. and Salakhutdinov, R. R. (2006). Reducing the dimensionality of data with neural networks. Science 313 504–7.

[35]

Kingma, D. P. and Ba, J. (2015). Adam: A method for stochastic optimization. In International conference on learning representations (ICLR).

[36]

Vincent, P., Larochelle, H., Bengio, Y. and Manzagol, P.-A. (2008). Extracting and composing robust features with denoising autoencoders. In ICML ’08 pp 1096–103. Association for Computing Machinery.

[37]

Kingma, D. P. and Welling, M. (2013). Auto-encoding variational bayes. CoRR abs/1312.6114.

[38]

Ham, J., Lee, D. D., Mika, S. and Schölkopf, B. (2004). A kernel view of the dimensionality reduction of manifolds. In Proceedings of the twenty-first international conference on machine learning (ICML) p 47. ACM.