A Riemannian geometric framework for manifold learning of non-Euclidean data

Abstract

A growing number of problems in data analysis and classification involve data that are non-Euclidean. For such problems, a naive application of vector space analysis algorithms will produce results that depend on the choice of local coordinates used to parametrize the data. At the same time, many data analysis and classification problems eventually reduce to an optimization, in which the criteria being minimized can be interpreted as the distortion associated with a mapping between two curved spaces. Exploiting this distortion minimizing perspective, we first show that manifold learning problems involving non-Euclidean data can be naturally framed as seeking a mapping between two Riemannian manifolds that is closest to being an isometry. A family of coordinate-invariant first-order distortion measures is then proposed that measure the proximity of the mapping to an isometry, and applied to manifold learning for non-Euclidean data sets. Case studies ranging from synthetic data to human mass-shape data demonstrate the many performance advantages of our Riemannian distortion minimization framework.

Appendices

Appendix

A: Further mathematical details of manifold learning algorithms

1.1 A.1 Proof of Proposition 1

Proof

The inverse metric $JG^{-1}J^{\top }$ at $x = x_i$ is obtained in (14) as

$$\begin{aligned} J G^{-1} J^{\top }(x_i) = \frac{1}{2} Y (\text {diag}(L_i) - e_i e_i^{\top } L - L^{\top } e_i e_i^{\top }) Y^{\top }, \end{aligned}$$

where $Y = \begin{bmatrix} y_1, \ldots , y_N \end{bmatrix} \in {\mathbb {R}}^{n\times N}$ is the matrix representation of the embeddings, $L = \frac{1}{ch}( {\tilde{D}}^{-1}{\tilde{K}} - I) \in {\mathbb {R}}^{N\times N}$ is the normalized graph Laplacian (${\tilde{D}}, {\tilde{K}}\in {\mathbb {R}}^{N\times N}$ are obtained from Algorithm 1 and both $c,h>0$), $L_i \in {\mathbb {R}}^N$ is the i-th row of L, and $e_i = (0,\ldots ,1,\ldots ,0) \in {\mathbb {R}}^N$ is a standard basis vector whose i-th component is one.

To see if $J G^{-1} J^{\top }(x_i)$ is positive semi-definite, it suffices to see if

$$\begin{aligned} M_i \equiv c \ h \left( \text {diag}(L_i) - e_i e_i^{\top } L - L^{\top } e_i e_i^{\top } \right) \in {\mathbb {R}}^{N\times N} \end{aligned}$$

(20)

is positive semi-definite. For any $v = (v_1, \ldots , v_N) \in {\mathbb {R}}^N$,

$$\begin{aligned} v^\top M_i v&= c \ h \left( \sum _{j,k=1}^N \left( \text {diag}(L_i) - e_i e_i^{\top } L - L^{\top } e_i e_i^{\top } \right) _{jk} v_j v_k \right) \end{aligned}$$

(21)

$$\begin{aligned}&= c \ h \left( \sum _{j=1}^N L_{ij} v_j^2 - 2 L_{ij} v_i v_j \right) \end{aligned}$$

(22)

$$\begin{aligned}&= v_i^2 + \sum _{j=1}^N ({\tilde{D}}_{ii})^{-1} {\tilde{K}}_{ij} (v_j^2 - 2v_i v_j) \end{aligned}$$

(23)

$$\begin{aligned}&= \sum _{j=1}^N ({\tilde{D}}_{ii})^{-1} {\tilde{K}}_{ij} (v_i - v_j)^2 \ge 0, \end{aligned}$$

(24)

where $L_{ij}$ denotes the (i, j) entry of L in (22). In deriving (23)-(24), we use the equalities $L_{ij} = \frac{1}{ch}(({\tilde{D}}_{ii})^{-1} {\tilde{K}}_{ij} - \delta _{ij})$ ($\delta _{ij} = 1$ if $i=j$ and 0 otherwise) and $\sum _{j=1}^N ({\tilde{D}}_{ii})^{-1} {\tilde{K}}_{ij} = 1$, and also the inequality ${\tilde{K}}_{ij} \ge 0$ for $i,j = 1,\ldots ,N$. Since the inequality $v^\top M_i v \ge 0$ holds for any $v \in {\mathbb {R}}^N$, $M_i$ is positive semi-definite; then $JG^{-1} J^\top (x_i) = \frac{1}{2ch} Y M_i Y^\top$ also becomes positive semi-definite for all $i = 1, \ldots , N$. $\square$

1.2 A.2 Riemannian relaxation

In the Riemannian relaxation method of McQueen et al. (2016), ${{\mathcal {M}}}$ is chosen to be an m-dimensional submanifold of Euclidean ambient space ${\mathbb {R}}^D$, with Riemannian metric G corresponding to the Euclidean metric on ${\mathbb {R}}^D$ projected to ${{\mathcal {M}}}$. The target manifold ${{\mathcal {N}}}$ is set to be ${\mathbb {R}}^n$ for some a priori chosen dimension $n \ge \text{ dim }({{\mathcal {M}}})$; the Riemannian metric on ${{\mathcal {N}}}$ is set to $H=I$.

Given Euclidean data points $u_i \in {\mathbb {R}}^D$, $i = 1, \ldots , N$ ($x_i \in {\mathbb {R}}^m$ in local coordinates), denote their n-dimensional embeddings by $y_i \in {\mathbb {R}}^n$. The embedding is then obtained as the solution to the following optimization:

$$\begin{aligned} \min _{y_i} \sum _{i=1}^N \Vert JG^{-1}J^{\top }(u_i) - I\Vert ^2 \alpha _i, \end{aligned}$$

(25)

where $JG^{-1}J^{\top }(u_i)$ denotes the $JG^{-1}J^{\top }$ estimated on $u_i$ using the method presented in Perrault-Joncas and Meila (2013), $\Vert \cdot \Vert$ denotes the matrix spectral norm, and $\alpha _i$ are weights. If $n > \text {dim}(\mathcal {M})$, I in (25) is replaced by $R_m R_m^{\top }$, where $R_m = [r_1, \ldots , r_m] \in {\mathbb {R}}^{n\times m}$ with $r_i \in {\mathbb {R}}^n$ the i-th singular vector of $JG^{-1}J^{\top }$.

From the perspective of our Riemannian distortion framework, assuming the rank of $J G^{-1} J^{\top }$ is m and the weights $\alpha _i$ in (25) are set to ${\tilde{d}}_{i}$ (obtained from Algorithm 1), the objective function in (25) can be expressed as

$$\begin{aligned} \min _{f} \mathcal {D}(f) = \int _{{\mathcal {M}}} \max _{i} (\lambda _i-1)^2 \sqrt{\det G} \, dx^1 \cdots dx^m, \end{aligned}$$

(26)

where the $\lambda _i$ are the m nonzero eigenvalues of $J G^{-1} J^{\top }$, which are identical to those of $J^{\top } J G^{-1}$. Since in practice the numerical estimation of $J G^{-1} J^{\top } \in {\mathbb {R}}^{n\times n}$ may yield a rank higher than m when $n > \text{ dim }({{\mathcal {M}}})$, one solution is to impose a soft constraint on the rank of $J G^{-1} J^{\top }$, e.g., in McQueen et al. (2016) the optimization is formulated as

$$\begin{aligned} \min _{f} \mathcal {D}(f) = \int _{{\mathcal {M}}} \max \left( \max _{i \in \mathcal {I}_m} (\lambda _i-1)^2, \ \max _{i \not \in \mathcal {I}_m} \left( \frac{\lambda _i}{\epsilon }\right) ^2\right) \sqrt{\det G} \, dx^1 \cdots dx^m, \end{aligned}$$

(27)

where $\lambda _i$ are the eigenvalues of $J G^{-1} J^{\top }$, $\mathcal {I}_m$ denotes the set of indices of the m largest eigenvalues, and $\epsilon >0$ is a scalar parameter intended to suppress the smaller $(n-m)$ eigenvalues.

1.3 A.3 Proof of Proposition 2

Proof

For $H = I$, the discretized formulation of the harmonic mapping distortion in the form of (15) is obtained as follows:

$$\begin{aligned} {{\mathcal {D}}}(Y)= & {} \sum _{i=1}^N \text {Tr}(JG^{-1}J^{\top }(x_i)) \ {\tilde{d}}_i \end{aligned}$$

(28)

$$\begin{aligned}= & {} \frac{1}{2} \sum _{i=1}^N \text {Tr} \left( Y(\text {diag} (L_i) - e_i e_i^{\top } L - L^{\top } e_i e_i^{\top })Y^{\top } \right) {\tilde{d}}_i \end{aligned}$$

(29)

$$\begin{aligned}= & {} \frac{1}{2}\text {Tr} \left( Y(\text {diag} ( \mathbb {1}_N^{\top } {\tilde{D}} L) - {\tilde{D}} L - L^{\top } {\tilde{D}} ) Y^{\top } \right) \end{aligned}$$

(30)

$$\begin{aligned}= & {} \frac{1}{c \ h} \text {Tr}(Y({\tilde{D}}-{\tilde{K}})Y^{\top }), \end{aligned}$$

(31)

where ${\tilde{K}}, {\tilde{d}}_i, {\tilde{D}}$ are obtained from Algorithm 1, $L = \frac{1}{c \ h} ({\tilde{D}}^{-1}{\tilde{K}} - I)\in {\mathbb {R}}^{N \times N}$ is the graph Laplacian from Algorithm 1, $L_i$ is the i-th row of L, $e_i = (0,\ldots ,1,\ldots ,0) \in {\mathbb {R}}^N$ is a standard basis vector whose i-th component is one, and $\mathbb {1}_N \in {\mathbb {R}}^N$ denotes an N-dimensional vector whose components are all one. In deriving (28)–(31), the estimate of $JG^{-1}J^{\top }$ at $x_i$ in (14), and the equalities $\text {Tr}(J^{\top } HJG^{-1}) = \text {Tr}(JG^{-1}J^{\top }H)$ and $\mathbb {1}_N^{\top } {\tilde{D}} L = 0$ are used.

Given a constant matrix $Y_b$ specified by the boundary condition, minimizing (31) for $Y_r$ reduces to

$$\begin{aligned} \min _{Y_r} \text{ Tr }(Y({\tilde{D}}-{\tilde{K}})Y^{\top }) = \text{ Tr }(Y_{b}({\tilde{D}}_{bb}-{\tilde{K}}_{bb})Y_{b}^{\top } - 2Y_{b}{\tilde{K}}_{br}Y_r^{\top } + Y_r({\tilde{D}}_{rr}-{\tilde{K}}_{rr})Y_r^{\top }). \end{aligned}$$

(32)

A closed-form solution for $Y_r$ is obtained as

$$\begin{aligned} Y_r = Y_b {\tilde{K}}_{br} ({\tilde{D}}_{rr}-{\tilde{K}}_{rr})^{-1} = Y_b W, \end{aligned}$$

(33)

where $W = {\tilde{K}}_{br} ({\tilde{D}}_{rr}-{\tilde{K}}_{rr})^{-1} \in {\mathbb {R}}^{N_b \times N_r}$.

Assume that ${\tilde{K}}_{ij} = {\tilde{K}}_{ji} \ge 0$ for all $i, j = 1, \ldots , N$, a graph with ${\tilde{K}}_{rr}$ as its adjacency matrix is connected, and ${\tilde{K}}_{br}$ is not a zero matrix. Then the matrix $({\tilde{D}}_{rr}-{\tilde{K}}_{rr})$ becomes positive-definite, so that W always exists. The positive-definiteness of $({\tilde{D}}_{rr}-{\tilde{K}}_{rr})$ can be shown from the following inequality: for any $v = (v_1, \ldots , v_{N_r})\ne 0 \in {\mathbb {R}}^{N_r}$,

$$\begin{aligned} v^\top ({\tilde{D}}_{rr}-{\tilde{K}}_{rr}) v&= \sum _{i,j = 1}^{N_r} ({\tilde{D}}_{rr}-{\tilde{K}}_{rr})_{ij} v_i v_j \end{aligned}$$

(34)

$$\begin{aligned}&= \sum _i = 1=^{N_r} ({\tilde{D}}_{rr})_{ii} v_i^2 - \sum _{i,j = 1}^{N_r} ({\tilde{K}}_{rr})_{ij} v_i v_j \end{aligned}$$

(35)

$$\begin{aligned}&= \sum _{i = 1}^{N_r} \left( \sum _{k = 1}^{N_b} ({\tilde{K}}_{br})_{ki} \right) v_i^2 + \frac{1}{2}\sum _{i,j = 1}^{N_r} ({\tilde{K}}_{rr})_{ij} (v_i - v_j)^2 > 0, \end{aligned}$$

(36)

where we use the fact that $({\tilde{D}}_{rr})_{ii} = \sum _{k=1}^{N_b} ({\tilde{K}}_{br})_{ki} + \sum _{k = 1}^{N_r} ({\tilde{K}}_{rr})_{ki}$ in deriving (36). From the direct application of Cramer’s rule, it can be shown that each entry of $({\tilde{D}}_{rr} - {\tilde{K}}_{rr})^{-1}$ is non-negative. Since every entry of ${\tilde{K}}_{br}$ is non-negative, all the entries of W are also non-negative. Furthermore, W satisfies the equation $\mathbb {1}_{N_r}^{\top } = \mathbb {1}_{N_b}^{\top } W$ from the equality ${\tilde{D}}_{rr} \mathbb {1}_{N_r} = {\tilde{K}}_{rr} \mathbb {1}_{N_r} + {\tilde{K}}_{br}^{\top } \mathbb {1}_{N_b}$; hence the entries of each column of W sum to one. $\square$

B: Experimental details for Section 4

1.1 B.1 Swiss roll

Here we explain further experimental details for the case study performed in Sect. 4.1. The data points are non-uniformly sampled; referring to the unfolded manifold in Fig. 3a, the density is set to oscillate along the horizontal axis, while uniform along the vertical axis. When choosing the initial parameter value $\theta _0$ for Algorithm 2, locality-preserving embeddings are preferable. As a choice for such an initial parameter value, we use two-dimensional embedding obtained from the Isomap (Tenenbaum et al. 2000). Any other embeddings that preserve locality can also be used as an initial guess, e.g., those from locally linear embedding (LLE; Roweis and Saul 2000), Laplacian eigenmap (LE; Belkin and Niyogi 2003), diffusion map (DM; Coifman and Lafon 2006), Hessian eigenmap (HLLE; Donoho and Grimes 2003), or local tangent space alignment (LTSA; Zhang and Zha 2004).

For the embedding obtained from the Isomap method, we test its five different scalings as the initial parameter value $\theta _0$ for Algorithm 3; we then choose the output embeddings that show the best match to the pairwise distances between ten nearest neighbors in the ambient space. Also note that the kernel bandwidth parameter h for the approximation of the graph Laplacian in Algorithm 2 is chosen to have the same order to the averaged nearest neighbor distance from each of the data points according to Lafon (2004).

1.2 B.2 Synthetic P(2) data

1.2.1 B.2.1 Details for the submanifold considered in Section 4.2

The tangent space of $\text {P(n)}$ at any $P \in \text {P(n)}$ can be identified with $\text {S(n)}$, the space of $n \times n$ symmetric matrices. Given $X, Y \in \text {S(n)}$, the affine-invariant Riemannian metric at P is defined by the inner product

$$\begin{aligned} \langle X, Y \rangle _{P} = \text {Tr}(P^{-1} X P^{-1} Y). \end{aligned}$$

(37)

Consider the following orthogonal decomposition of $P \in \text {P(2)}$:

$$\begin{aligned} P = RSR^{\top }, \end{aligned}$$

(38)

where $R = \left[ \begin{array}{rr} \cos \theta &{} \quad -\sin \theta \\ \sin \theta &{} \quad \cos \theta \end{array} \right] \in \text {SO(2)}$ with $\theta \in [0, \frac{\pi }{2})$, and $S = \text {diag}(e^{p}, e^{q})$ for scalar p, q. A local coordinate chart can be defined in terms of $(p,q,\theta )$ on an open subset $U = \{P \in \text {P(2)} \, | \, P \ne c I \ \ \text {for} \ \ c > 0\}$. The affine-invariant Riemannian metric in (37) is then represented in $(p,q,\theta )$-coordinates (at $p \ne q$) as

$$\begin{aligned} ds^2 = \text {Tr}((P^{-1}dP)^2) = {dp}^2 + {dq}^2 + 2\left( e^{p-q} + e^{q-p} - 2\right) d\theta ^2. \end{aligned}$$

(39)

For the case study in Sect. 4.2, the data set shown in Fig. 4a is generated by joining two cylinders (with a hole) $\mathcal {C}_1$ and $\mathcal {C}_2$ in $(p, q, \theta )$-coordinates, where $\mathcal {C}_1 = \{(p, q, \theta ) \ | \ p = \sin \theta _S, \; q = -1 + \cos \theta _S, \; \theta _S \in \left[ 0, \frac{4}{3}\pi \right] , \; \theta \in \left[ 0, \frac{\pi }{4} \right] \}$ and $\mathcal {C}_2 = \{(p, q, \theta ) \ | \ p = \sin \theta _S, \; q = 1 - \cos \theta _S, \; \theta _S \in \left[ -\frac{4}{3}\pi , 0 \right] , \; \theta \in \left[ 0, \frac{\pi }{4} \right] \}$ (see Fig. 4a; the backbone curve in the figure corresponds to the direction along which $\theta _S$ varies). The affine-invariant Riemannian metric on this submanifold (at $\theta _S \ne 0$) is obtained in terms of coordinates $(\theta , \theta _S)$, $\theta _S \ne 0$, as

$$\begin{aligned} ds^2 = d\theta _S^2 + 2\left( e^{1 + \sqrt{2}\sin (|\theta _S| - \frac{\pi }{4})} + e^{-1 - \sqrt{2}\sin (|\theta _S| - \frac{\pi }{4})} - 2\right) d\theta ^2. \end{aligned}$$

(40)

Because of the nonzero Riemannian curvature of this submanifold, isometric embeddings in two-dimensional Euclidean space do not exist.

1.2.2 B.2.2 Evaluation of the pairwise distance and tangent vector angle errors

For data points $x_i \in \text {P(2)}$ and corresponding embeddings $y_i \in {\mathbb {R}}^2$, $i = 1, \ldots , N$, the pairwise distance error for k nearest points is defined as

$$\begin{aligned} \text {Pairwise distance error} = \frac{1}{k N}\sum _{i=1}^N\sum _{j\in \{\text {NN}_k(i)\}}\left( \Vert y_i-y_j\Vert - \text {dist}(x_i, x_j)\right) ^2, \end{aligned}$$

(41)

where $y_i$ denotes the optimal embedding of $x_i$, $\text {NN}_{k}(i)$ denotes the set of indices of k nearest neighbor points to $x_i$, and $\text {dist}(x_i,x_j)$ denotes the ground truth distance between $x_i$ and $x_j$, i.e., the geodesic distance measured on the submanifold. To measure angles, the tangent vectors are approximated by the difference between nearest neighbors. The tangent vector angle error is defined as

$$\begin{aligned}&\text {Tangent vector angle error} \nonumber \\&\quad = \frac{2}{k(k-1)N}\sum _{l=1}^N\sum _{i,j \in \{\text {NN}_{k}(l)\}} \left( \text {acos}(\langle v_i, v_j\rangle ) - \text {acos}(\langle V_i, V_j\rangle )\right) ^2, \end{aligned}$$

(42)

where $v_i, V_i$ respectively denote the tangent vector from the l-th data point to the i-th data point in the optimal embeddings and the original data points, and $\langle \cdot , \cdot \rangle$ denotes the inner product.

When reporting the final manifold learning results in Table 1, for the reference values to evaluate the pairwise distance error, we numerically obtain the minimal geodesic distances on the submanifolds. Also, the inner product in (37) is used to calculate the reference values for the angles between tangent vectors.

1.3 B.3 Human mass-inertia data

1.3.1 B.3.1 Synthesizing human mass-inertia data

Since mass-inertial parameter data for humans are not readily available, we use human shape data from Yang et al. (2014) to synthesize this data set; specifically, assuming uniform mass density, we integrate the volumes of the human body shapes to construct mass-inertial parameter data for the corresponding $N_l=10$ links.

Fig. 8

Principal components (PCs) of the human mass-inertia data obtained from PGA and vector space PCA. a–b respectively depict the third and fourth PCs obtained from PGA for the range of $\pm 4$ standard deviations from the mean; c–d respectively depict the third and fourth PCs obtained from PCA for different range of standard deviations from the mean

1.3.2 B.3.2 Further principal components of human mass-inertia data

As a supplement of Fig. 5 in Sect. 4.3, here we provide the third and fourth principal components of the human mass-inertia data obtained from both principal geodesic analysis (PGA) and vector space principal component analysis (PCA). The variations corresponding to the third and fourth principal components of PGA are shown in Fig. 8a–b. Principal component 3 captures variations in the height and torso thickness, and principal component 4 captures variations mainly in the height.

In the case of vector space PCA shown in Fig. 8c–d, the variations near the mean are qualitatively similar to those obtained for PGA. However, the positive-definiteness requirement is violated even for data points just 0.5 standard deviations away from the mean. The ellipsoids for those inertial parameters collapse and can be observed in the dashed red ellipses of Fig. 8c–d.