# LDLE: Low Distortion Local Eigenmaps

Dhruv Kohli
<EMAIL_ADDRESS>
Department of Mathematics
University of California San Diego
CA 92093, USA  Alexander Cloninger
<EMAIL_ADDRESS>
Department of Mathematics
University of California San Diego
CA 92093, USA  Gal Mishne
<EMAIL_ADDRESS>
Halicioğlu Data Science Institute
University of California San Diego
CA 92093, USA

###### Abstract

We present Low Distortion Local Eigenmaps (LDLE), a manifold learning
technique which constructs a set of low distortion local views of a dataset in
lower dimension and registers them to obtain a global embedding. The local
views are constructed using the global eigenvectors of the graph Laplacian and
are registered using Procrustes analysis. The choice of these eigenvectors may
vary across the regions. In contrast to existing techniques, LDLE can embed
closed and non-orientable manifolds into their intrinsic dimension by tearing
them apart. It also provides gluing instruction on the boundary of the torn
embedding to help identify the topology of the original manifold. Our
experimental results will show that LDLE largely preserved distances up to a
constant scale while other techniques produced higher distortion. We also
demonstrate that LDLE produces high quality embeddings even when the data is
noisy or sparse.

Keywords: Manifold learning, graph Laplacian, local parameterization,
Procrustes analysis, closed manifold, non-orientable manifold

##  1 Introduction

Manifold learning techniques such as Local Linear Embedding [37], Diffusion
maps [17], Laplacian eigenmaps [3], t-SNE [30] and UMAP [32], aim at
preserving local information as they map a manifold embedded in higher
dimension into lower (possibly intrinsic) dimension. In particular, UMAP and
t-SNE follow a top-down approach as they start with an initial low-dimensional
global embedding and then refine it by minimizing a local distortion measure
on it. In contrast, similar to LTSA [49] and [40], a bottom-up approach for
manifold learning can be imagined to consist of two steps, first obtaining low
distortion local views of the manifold in lower dimension and then registering
them to obtain a global embedding of the manifold. In this paper, we take this
bottom-up perspective to embed a manifold in low dimension, where the local
views are obtained by constructing coordinate charts for the manifold which
incur low distortion.

###  1.1 Local Distortion

Let $(\mathcal{M},g)$ be a $d$-dimensional Riemannian manifold with finite
volume. By definition, for every $x_{k}$ in $\mathcal{M}$, there exists a
coordinate chart $(\mathcal{U}_{k},\Phi_{k})$ such that
$x_{k}\in\mathcal{U}_{k}$, $\mathcal{U}_{k}\subset M$ and $\Phi_{k}$ maps
$\mathcal{U}_{k}$ into $\mathbb{R}^{d}$. One can imagine $\mathcal{U}_{k}$ to
be a local view of $\mathcal{M}$ in the ambient space. Using rigid
transformations, these local views can be registered to recover $\mathcal{M}$.
Similarly, $\Phi_{k}(\mathcal{U}_{k})$ can be imagined to be a local view of
$\mathcal{M}$ in the $d$-dimensional embedding space $\mathbb{R}^{d}$. Again,
using rigid transformations, these local views can be registered to obtain the
$d$-dimensional embedding of $\mathcal{M}$.

As there may exist multiple mappings which map $\mathcal{U}_{k}$ into
$\mathbb{R}^{d}$, a natural strategy would be to choose a mapping with low
distortion. Multiple measures of distortion exist in literature [14]. The
measure of distortion used in this work is as follows. Let $d_{g}(x,y)$ denote
the shortest geodesic distance between $x,y\in\mathcal{M}$. The distortion of
$\Phi_{k}$ on $\mathcal{U}_{k}$ as defined in [25] is given by

$\displaystyle\text{Distortion}(\Phi_{k},\mathcal{U}_{k})=\left\|\Phi_{k}\right\|_{\text{Lip}}\left\|\Phi_{k}^{-1}\right\|_{\text{Lip}}$
(1)

where $\left\|\Phi_{k}\right\|_{\text{Lip}}$ is the Lipschitz norm of
$\Phi_{k}$ given by

$\displaystyle\left\|\Phi_{k}\right\|_{\text{Lip}}$
$\displaystyle=\sup_{\begin{subarray}{c}x,y\in\mathcal{U}_{k}\\\ x\neq
y\end{subarray}}\frac{\left\|\Phi_{k}(x)-\Phi_{k}(y)\right\|_{2}}{d_{g}(x,y)},$
(2)

and similarly,

$\displaystyle\left\|\Phi^{-1}_{k}\right\|_{\text{Lip}}$
$\displaystyle=\sup_{\begin{subarray}{c}x,y\in\mathcal{U}_{k}\\\ x\neq
y\end{subarray}}\frac{d_{g}(x,y)}{\left\|\Phi_{k}(x)-\Phi_{k}(y)\right\|_{2}}.$
(3)

Note that $\text{Distortion}(\Phi_{k},\mathcal{U}_{k})$ is always greater than
or equal to $1$. If $\text{Distortion}(\Phi_{k},\mathcal{U}_{k})=1$, then
$\Phi_{k}$ is said to have no distortion on $\mathcal{U}_{k}$. This is
achieved when the mapping $\Phi_{k}$ preserves distances between points in
$\mathcal{U}_{k}$ up to a constant scale, that is, when $\Phi_{k}$ is a
similarity on $\mathcal{U}_{k}$. It is not always possible to obtain a mapping
with no distortion. For example, there does not exist a similarity which maps
a locally curved region on a surface into a Euclidean plane. This follows from
the fact that the sign of the Gaussian curvature is preserved under similarity
transformation which in turn follows from the Gauss’s Theorema Egregium.

###  1.2 Our Contributions

This paper takes motivation from the work in [25] where the authors provide
guarantees on the distortion of the coordinate charts of the manifold
constructed using carefully chosen eigenfunctions of the Laplacian. However,
this only applies to the charts for small neighborhoods on the manifold and
does not provide a global embedding. In this paper, we present an approach to
realize their work in the discrete setting and obtain low-dimensional low
distortion local views of the given dataset using the eigenvectors of the
graph Laplacian. Moreover, we piece together these local views to obtain a
global embedding of the manifold. The main contributions of our work are as
follows:

  1. 1.

We present an algorithmic realization of the construction procedure in [25]
that applies to the discrete setting and yields low-dimensional low distortion
views of small metric balls on the given discretized manifold (See Section 2
for a summary of their procedure).

  2. 2.

We present an algorithm to obtain a global embedding of the manifold by
registering its local views. The algorithm is designed so as to embed closed
as well as non-orientable manifolds into their intrinsic dimension by tearing
them apart. It also provides gluing instructions for the boundary of the
embedding by coloring it such that the points on the boundary which are
adjacent on the manifold have the same color (see Figure 2).

LDLE consists of three main steps. In the first step, we estimate the inner
product of the Laplacian eigenfunctions’ gradients using the local correlation
between them. These estimates are used to choose eigenfunctions which are in
turn used to construct low-dimensional low distortion parameterizations
$\Phi_{k}$ of the small balls $U_{k}$ on the manifold. The choice of the
eigenfunctions depend on the underlying ball. A natural next step is to align
these local views $\Phi_{k}(U_{k})$ in the embedding space, to obtain a global
embedding. One way to align them is to use Generalized Procrustes Analysis
(GPA) [18, 20, 43]. However, we empirically observed that GPA is less
efficient and prone to errors due to large number of local views with small
overlaps between them. Therefore, motivated from our experimental observations
and computational necessity, in the second step, we develop a clustering
algorithm to obtain a small number of intermediate views
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ with low distortion, from the large
number of smaller local views $\Phi_{k}(U_{k})$. This makes the subsequent GPA
based registration procedure faster and less prone to errors.

Finally, in the third step, we register intermediate views
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ using an adaptation of GPA which
enables tearing of closed and non-orientable manifolds so as to embed them
into their intrinsic dimension. The results on a 2D rectangular strip and a 3D
sphere are presented in Figures 1 and 2, to motivate our approach.

The paper organization is as follows. Section 2 provides relevant background
and motivation. In Section 3 we present the construction of low-dimensional
low distortion local parameterizations. Section 4 presents our clustering
algorithm to obtain intermediate views. Section 5 registers the intermediate
views to a global embedding. In Section 6 we compare the embeddings produced
by our algorithm with existing techniques on multiple datasets. Section 7
concludes our work and discusses future directions.

###  1.3 Related Work

Laplacian eigenfunctions are ubiquitous in manifold learning. A large
proportion of the existing manifold learning techniques rely on a fixed set of
Laplacian eigenfunctions, specifically, on the first few non-trivial low
frequency eigenfunctions, to construct a low-dimensional embedding of a
manifold in high dimensional ambient space. These low frequency eigenfunctions
not only carry information about the global structure of the manifold but they
also exhibit robustness to the noise in the data [17]. Laplacian eigenmaps
[3], Diffusion maps [17] and UMAP [32] are examples of such top-down manifold
learning techniques. While there are limited bottom-up manifold learning
techniques in the literature, to the best of our knowledge, none of them makes
use of Laplacian eigenfunctions to construct local views of the manifold in
lower dimension.

##### LTSA

is an example of a bottom-up approach for manifold learning whose local
mappings project local neighborhoods onto the respective tangential spaces. A
local mapping in LTSA is a linear transformation whose columns are the
principal directions obtained by applying PCA on the underlying neighborhood.
These directions form an estimate of the basis for the tangential space.
Having constructed low-dimensional local views for each neighborhood, LTSA
then aligns all the local views to obtain a global embedding. As discussed in
their work and as we will show in our experimental results, LTSA lacks
robustness to the noise in the data. This further motivates our approach of
using robust low-frequency Laplacian eigenfunctions for the construction of
local views. Moreover, due to the specific constraints used in their
alignment, LTSA embeddings fail to capture the aspect ratio of the underlying
manifold (see Appendix F for details).

##### Laplacian eigenmaps

uses the eigenvectors corresponding to the $d$ smallest eigenvalues (excluding
zero) of the normalized graph Laplacian to embed the manifold in
$\mathbb{R}^{d}$. It can also be perceived as a top-down approach which
directly obtains a global embedding that minimizes Dirichlet energy under some
constraints. For manifolds with high aspect ratio, in the context of Section
1.1, the distortion of the local parameterizations based on the restriction of
these eigenvectors on local neighborhoods, could become extremely high. For
example, as shown in Figure 1, the Laplacian eigenmaps embedding of a
rectangle with an aspect ratio of $16$ looks like a parabola. This issue is
explained in detail in [38, 8, 19, 6].

##### UMAP,

to a large extent, resolves this issue by first computing an embedding based
on the $d$ non-trivial low-frequency eigenvectors of a symmetric normalized
Laplacian and then “sprinkling” white noise in it. It then refines the noisy
embedding by minimizing a local distortion measure based on fuzzy set cross
entropy. Although UMAP embeddings seem to be topologically correct, they
occasionally tend to have twists and sharp turns which may be unwanted (see
Figure 1).

##### t-SNE

takes a different approach of randomly initializing the global embedding,
defining a local t-distribution in the embedding space and local Gaussian
distribution in the high dimensional ambient space, and finally refining the
embedding by minimizing the Kullback–Leibler divergence between the two sets
of distributions. As shown in Figure 1, t-SNE tends to output a dissected
embedding even when the manifold is connected. Note that the recent work by
[26] showed that t-SNE with spectral initialization results in a similar
embedding as that of UMAP. Therefore, in this work, we display the output of
the classic t-SNE construction, with random initialization only.

Input |  LDLE |  LDLE with $\partial\mathcal{M}$ known apriori |  LTSA |  UMAP |  t-SNE |  Laplacian Eigenmaps
---|---|---|---|---|---|---
|  |  |  |  |  |
Figure 1: Embeddings of a rectangle ($4\times 0.25$) with high aspect ratio in
$\mathbb{R}^{2}$ into $\mathbb{R}^{2}$.

A missing feature in existing manifold learning techniques is their ability to
embed closed manifolds into their intrinsic dimensions. For example, a sphere
in $\mathbb{R}^{3}$ is a $2$-dimensional manifold which can be represented by
a connected domain in $\mathbb{R}^{2}$ with boundary gluing instructions
provided in the form of colors. We solve this issue in this paper (see Figure
2).

Input |  LDLE |  LTSA |  UMAP |  t-SNE |  Laplacian Eigenmaps
---|---|---|---|---|---
|  |  |  |  |
|  |  |  |  |
Figure 2: Embeddings of a sphere in $\mathbb{R}^{3}$ into $\mathbb{R}^{2}$.
The top and bottom row contain the same plots colored by the height and the
azimuthal angle of the sphere ($0-2\pi$), respectively. LDLE automatically
colors the boundary so that the points on the boundary which are adjacent on
the sphere have the same color. The arrows are manually drawn to help the
reader identify the two pieces of the boundary which are to be stitched
together to recover the original sphere. LTSA, UMAP and Laplacian eigenmaps
squeezed the sphere into different viewpoints of $\mathbb{R}^{2}$ (side or top
view of the sphere). t-SNE also tore apart the sphere but the embedding lacks
interpretability as it is “unaware” of the boundary.

##  2 Background and Motivation

Due to their global nature and robustness to noise, in our bottom-up approach
for manifold learning, we propose to construct low distortion (see Eq. (1))
local mappings using low frequency Laplacian eigenfunctions. A natural way to
achieve this is to restrict the eigenfunctions on local neighborhoods.
Unfortunately, the common trend of using first $d$ non-trivial low frequency
eigenfunctions to construct these local mappings fails to produce low
distortion on all neighborhoods. This directly follows from the Laplacian
Eigenmaps embedding of a high aspect-ratio rectangle shown in Figure 1. The
following example explains that even in case of unit aspect-ratio, a local
mapping based on the same set of eigenfunctions would not incur low distortion
on each neighborhood, while mappings based on different sets of eigenfunctions
may achieve that.

|
---|---
|
Figure 3: (Left) Distortion of $\Phi_{1}^{*}$ (top) and $\Phi_{2}^{*}$
(bottom) on discs of radius $0.01$ centered at $(x,y)$ for all
$x,y\in[0,1]\times[0,1]$. $\Phi_{2}^{*}$ produces close to infinite distortion
on the discs located in the white region. (Right) Mapping of the discs at
various locations in the square using $\Phi_{1}^{*}$ (top) and $\Phi_{2}^{*}$
(bottom).

Consider a unit square $[0,1]\times[0,1]$ such that for every point $x_{k}$ in
the square, $\mathcal{U}_{k}$ is the disc of radius $0.01$ centered at
$x_{k}$. Consider a mapping $\Phi_{1}^{*}$ based on the first two non-trivial
eigenfunctions $\cos(\pi x)$ and $\cos(\pi y)$ of the Laplace-Beltrami
operator on the square with Neumann boundary conditions, that is,

$\displaystyle\Phi_{1}^{*}(x,y)=(\cos(\pi x),\cos(\pi y)).$ (4)

As shown in Figure 3, $\Phi_{1}^{*}$ maps the discs along the diagonals to
other discs. The discs along the horizontal and vertical lines through the
center are mapped to ellipses. The skewness of these ellipses increases as we
move closer to the middle of the edges of the unit square. Thus, the
distortion of $\Phi_{1}^{*}$ is low on the discs along the diagonals and high
on the discs close to the middle of the edges of the square.

Now, consider a different mapping based on another set of eigenfunctions,

$\displaystyle\Phi_{2}^{*}(x,y)$ $\displaystyle=(\cos(5\pi x),\cos(5\pi y)).$
(5)

Compared to $\Phi_{1}^{*}$, $\Phi_{2}^{*}$ produces almost no distortion on
the discs of radius $0.01$ centered at $(0.1,0.5)$ and $(0.9,0.5)$ (see Figure
3). Therefore, in order to achieve low distortion, it seem to make sense to
construct local mappings for different regions based on different sets of
eigenfunctions.

The following result from [25] manifests the above claim as it shows that, for
a given small neighborhood on a Riemannian manifold, there always exist a
subset of Laplacian eigenfunctions such that a local parameterization based on
this subset is bilipschitz and has bounded distortion. A more precise
statement follows.

######  Theorem 1 ([25], Theorem 2.2.1).

Let $(\mathcal{M},g)$ be a $d$-dimensional Riemannian manifold. Let
$\Delta_{g}$ be the Laplace-Beltrami operator on it with Dirichlet or Neumann
boundary conditions and let $\phi_{i}$ be an eigenfunction of $\Delta_{g}$
with eigenvalue $\lambda_{i}$. Assume that $|\mathcal{M}|=1$ where
$|\mathcal{M}|$ is the volume of $\mathcal{M}$ and the uniform ellipticity
conditions for $\Delta_{g}$ are satisfied. Let $x_{k}\in\mathcal{M}$ and
$r_{k}$ be less than the injectivity radius at $x_{k}$ (the maximum radius
where the the exponential map is a diffeomorphism). Then, there exists a
constant $\kappa>1$ which depends on $d$ and the metric tensor $g$ such that
the following hold. Let $\rho\leq r_{k}$ and $B_{k}\equiv
B_{\kappa^{-1}\rho}(x_{k})$ where

$\displaystyle B_{\epsilon}(x)$ $\displaystyle=\\{y\in\mathcal{M}\ |\
d_{g}(x,y)<\epsilon\\}.$ (6)

Then there exist $i_{1},i_{2},\ldots,i_{d}$ such that, if we let

$\displaystyle\gamma_{ki}=\left(\frac{\int_{B_{k}}\phi_{i}^{2}(y)dy}{|B_{k}|}\right)^{-1/2}$
(7)

then the map

$\displaystyle\Phi_{k}:B_{k}$ $\displaystyle\rightarrow\mathbb{R}^{d}$
$\displaystyle x$
$\displaystyle\rightarrow(\gamma_{ki_{1}}\phi_{i_{1}}(x),\ldots,\gamma_{ki_{d}}\phi_{i_{d}}(x))$
(8)

is bilipschitz such that for any $y_{1},y_{2}\in B_{k}$ it satisfies

$\displaystyle\frac{\kappa^{-1}}{\rho}d_{g}(y_{1},y_{2})\leq\left\|\Phi_{k}(y_{1})-\Phi_{k}(y_{2})\right\|\leq\frac{\kappa}{\rho}d_{g}(y_{1},y_{2}),$
(9)

where the associated eigenvalues satisfy

$\displaystyle\kappa^{-1}\rho^{-2}\leq\lambda_{i_{1}},\ldots,\lambda_{i_{d}}\leq\kappa\rho^{-2},$
(10)

and the distortion is bounded from above by $\kappa^{2}$ i.e.

$\displaystyle\sup_{\begin{subarray}{c}y_{1},y_{2}\in B_{k}\\\ y_{1}\neq
y_{2}\end{subarray}}\frac{\left\|\Phi_{k}(y_{1})-\Phi_{k}(y_{2})\right\|}{d_{g}(y_{1},y_{2})}\sup_{\begin{subarray}{c}y_{1},y_{2}\in
B_{k}\\\ y_{1}\neq
y_{2}\end{subarray}}\frac{d_{g}(y_{1},y_{2})}{\left\|\Phi_{k}(y_{1})-\Phi_{k}(y_{2})\right\|}\leq\frac{\kappa}{\rho}\frac{\rho}{\kappa^{-1}}=\kappa^{2}.$
(11)

Motivated by the above result, we adopt the form of local paramterizations
$\Phi_{k}$ in Eq. (8) as local mappings in our work. The main challenge then
is to identify the set of eigenfunctions for a given neighborhood such that
the resulting parameterization produces low distortion on it. The existence
proof of the above theorem by the authors of [25] suggests a procedure to
identify this set in the continuous setting. Below, we provide a sketch of
their procedure and in Section 3 we describe our discrete realization of it.

###  2.1 Eigenfunction Selection in the Continuous Setting

Before describing the procedure used in [25] to choose the eigenfunctions, we
first provide some intuition about the desired properties for the chosen
eigenfunctions $\phi_{i_{1}},\ldots,\phi_{i_{d}}$ so that the resulting
parameterization $\Phi_{k}$ has low distortion on $B_{k}$.

Consider the simple case of $B_{k}$ representing a small open ball of radius
$\kappa^{-1}\rho$ around $x_{k}$ in $\mathbb{R}^{d}$ equipped with the
standard Euclidean metric. Then the first-order Taylor approximation of
$\Phi_{k}(x)$, $x\in B_{k}$, about $x_{k}$ is given by

$\displaystyle\Phi_{k}(x)$
$\displaystyle\approx\Phi_{k}(x_{k})+J(x-x_{k})\text{ where
}J=[\gamma_{ki_{1}}\nabla\phi_{i_{1}}(x_{k})\ldots\gamma_{ki_{d}}\nabla\phi_{i_{d}}(x_{k})]^{T}.$
(12)

Note that $\gamma_{ki_{s}}$ are positive scalars constant with respect to $x$.
Now, $\text{Distortion}(\Phi_{k},B_{k})=1$ if and only if $\Phi_{k}$ preserves
distances between points in $B_{k}$ up to a constant scale (see Eq. (1)). That
is,

$\displaystyle\left\|\Phi_{k}(x)-\Phi_{k}(y)\right\|_{2}=c\left\|x-y\right\|_{2}\
\forall x,y\in B_{k}\text{ and for some constant }c>0.$ (13)

Using the first-order approximation of $\Phi_{k}$ we get,

$\displaystyle\left\|J(x-y)\right\|_{2}\approx c\left\|x-y\right\|_{2}\
\forall x,y\in B_{k}\text{ and for some constant }c>0.$ (14)

Therefore, for low distortion $\Phi_{k}$, $J$ must approximately behave like a
similarity transformation and therefore, $J$ needs to be approximately
orthogonal up to a constant scale. In other words, the chosen eigenfunctions
should be such that $\gamma_{ki_{1}}\nabla\phi_{i_{1}}(x_{k}),$ $\ldots,$
$\gamma_{ki_{d}}\nabla\phi_{i_{d}}(x_{k})$ are close to being orthogonal and
have similar lengths.

The same intuition holds in the manifold setting too. The construction
procedure described in [25] aims to choose eigenfunctions such that

  1. (a)

they are close to being locally orthogonal, that is,
$\nabla\phi_{i_{1}}(x_{k}),\ldots,\nabla\phi_{i_{d}}(x_{k})$ are approximately
orthogonal, and

  2. (b)

that their local scaling factors
$\gamma_{ki_{s}}\left\|\nabla\phi_{i_{s}}(x_{k})\right\|_{2}$ are close to
each other.

Note. Throughout this paper, we use the convention
$\nabla\phi_{i}(x_{k})=\nabla(\phi_{i}\ \circ\ \exp_{x_{k}})(0)$ where
$\exp_{x_{k}}$ is the exponential map at $x_{k}$. Therefore,
$\nabla\phi_{i}(x_{k})$ can be represented by a $d$-dimensional vector in a
given $d$-dimensional orthonormal basis of $T_{x_{k}}\mathcal{M}$. Even though
the representation of these vectors depend on the choice of the orthonormal
basis, the value of the canonical inner product between these vectors, and
therefore the $2$-norm of the vectors, are the same across different basis.
This follows from the fact that an orthogonal transformation preserves the
inner product.

######  Remark 1.

Based on the above first order approximation, one may take our local mappings
$\Phi_{k}$ to also be projections onto the tangential spaces. However, unlike
LTSA [49] where the basis of the tangential space is estimated by the local
principal directions, in our case it is estimated by the locally orthogonal
gradients of the global eigenfunctions of the Laplacian. Therefore, LTSA
relies only on the local structure to estimate the tangential space while, in
a sense, our method makes use of both local and global structure of the
manifold.

A high level overview of the procedure presented in [25] to choose
eigenfunctions which satisfy the properties in (a) and (b) follows.

  1. 1.

A set $S_{k}$ of the indices of candidate eigenfunctions is chosen such that
$i\in S_{k}$ if the length of $\gamma_{ki}\nabla\phi_{i}(x_{k})$ is bounded
from above by a constant, say $C$.

  2. 2.

A direction $p_{1}\in T_{x_{k}}\mathcal{M}$ is selected at random.

  3. 3.

Subsequently $i_{1}\in S_{k}$ is selected so that
$\gamma_{ki_{1}}|\nabla\phi_{i_{1}}(x_{k})^{T}p_{1}|$ is sufficiently large.
This motivates $\gamma_{ki_{1}}\nabla\phi_{i_{1}}(x_{k})$ to be approximately
in the same direction as $p_{1}$ and the length of it to be close to the upper
bound $C$.

  4. 4.

Then, a recursive strategy follows. To find the $s$-th eigenfunction for
$s\in\\{2,\ldots,d\\}$, a direction $p_{s}\in T_{x_{k}}\mathcal{M}$ is chosen
such that it is orthogonal to
$\nabla\phi_{i_{1}}(x_{k}),\ldots,\nabla\phi_{i_{s-1}}(x_{k})$.

  5. 5.

Subsequently, $i_{s}\in S_{k}$ is chosen so that
$\gamma_{ki_{s}}|\nabla\phi_{i_{s}}(x_{k})^{T}p_{s}|$ is sufficiently large.
Again, this motivates $\gamma_{ki_{s}}\nabla\phi_{i_{s}}(x_{k})$ to be
approximately in the same direction as $p_{s}$ and the length of it to be
close to the upper bound $C$.

Since $p_{s}$ is orthogonal to
$\nabla\phi_{i_{1}}(x_{k}),\ldots,\nabla\phi_{i_{s-1}}(x_{k})$ and the
direction of $\gamma_{ki_{s}}\nabla\phi_{i_{s}}$ is approximately the same as
$p_{s}$, therefore $(a)$ is satisfied. Since for all $s\in\\{1,\ldots,d\\}$,
$\gamma_{ki_{s}}\nabla\phi_{i_{s}}(x_{k})$ has a length close to the upper
bound $C$, therefore $(b)$ is also satisfied. The core of their work lies in
proving that these $\phi_{i_{1}},\ldots,\phi_{i_{d}}$ always exist under the
assumptions of the theorem such that the resulting parameterization $\Phi_{k}$
has bounded distortion (see Eq. (11)). This bound depends on the intrinsic
dimension $d$ and the natural geometric properties of the manifold. The main
challenge in practically realizing the above procedure lies in the estimation
of $\nabla\phi_{i_{s}}(x_{k})^{T}p_{s}$. In Section 3, we overcome this
challenge.

##  3 Low-dimensional Low Distortion Local Parameterization

In the procedure to choose $\phi_{i_{1}},\ldots,\phi_{i_{d}}$ to construct
$\Phi_{k}$ as described above, the selection of the first eigenfunction
$\phi_{i_{1}}$ relies on the derivative of the eigenfunctions at $x_{k}$ along
an arbitrary direction $p_{1}\in T_{x_{k}}\mathcal{M}$, that is, on
$\nabla\phi_{i}(x_{k})^{T}p_{1}$. In our algorithmic realization of the
construction procedure, we take $p_{1}$ to be the gradient of an eigenfunction
at $x_{k}$ itself (say $\nabla\phi_{j}(x_{k})$). We relax the unit norm
constraint on $p_{1}$; note that this will neither affect the math nor the
output of our algorithm. Then the selection of $\phi_{i_{1}}$ would depend on
the inner products $\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k})$. The value
of this inner product does not depend on the choice of the orthonormal basis
for $T_{x_{k}}\mathcal{M}$. We discuss several ways to obtain a numerical
estimate of this inner product by making use of the local correlation between
the eigenfunctions [42, 16]. These estimates are used to select the subsequent
eigenfunctions too.

In Section 3.1, we first review the local correlation between the
eigenfunctions of the Laplacian. In Theorem 2 we show that the limiting value
of the scaled local correlation between two eigenfunctions equals the inner
product of their gradients. We provide two proofs of the theorem where each
proof leads to a numerical procedure described in Section 3.2, followed by
examples to empirically compare the estimates. Finally, in Section 3.3, we use
these estimates to obtain low distortion local parameterizations of the
underlying manifold.

###  3.1 Inner Product of Eigenfunction Gradients using Local Correlation

Let $(\mathcal{M},g)$ be a $d$-dimensional Riemannian manifold with or without
boundary, rescaled so that $|\mathcal{M}|\leq 1$. Denote the volume element at
$y$ by $\omega_{g}(y)$. Let $\phi_{i}$ and $\phi_{j}$ be the eigenfunctions of
the Laplacian operator $\Delta_{g}$ (see statement of Theorem 1) with
eigenvalues $\lambda_{i}$ and $\lambda_{j}$. Let $x_{k}\in\mathcal{M}$ and
define

$\displaystyle\Psi_{kij}(y)=(\phi_{i}(y)-\phi_{i}(x_{k}))(\phi_{j}(y)-\phi_{j}(x_{k})).$
(15)

Then the local correlation between the two eigenfunctions $\phi_{i}$ and
$\phi_{j}$ at the point $x_{k}$ at scale $t_{k}^{-1/2}$ as defined in [42, 16]
is given by

$\displaystyle
A_{kij}=\int_{\mathcal{M}}p(t_{k},x_{k},y)\Psi_{kij}(y)\omega_{g}(y),$ (16)

where $p(t,x,y)$ is the fundamental solution of the heat equation on
$(\mathcal{M},g)$. As noted in [42], for $(t_{k},x_{k})\in\mathbb{R}_{\geq
0}\times\mathcal{M}$ fixed, we have

$\displaystyle
p(t_{k},x_{k},y)\sim\left\\{\begin{matrix}[l]t_{k}^{-d/2}&d_{g}(x_{k},y)\leq
t_{k}^{-1/2}\\\
0&\text{otherwise}\end{matrix}\right.\qquad\text{and}\qquad\int_{M}p(t_{k},x_{k},y)\omega_{g}(y)=1.$
(17)

Therefore, $p(t_{k},x_{k},\cdot)$ acts as a local probability measure centered
at $x_{k}$ with scale $t_{k}^{-1/2}$ (see Eq. (67) in Appendix A for a precise
form of $p$). We define the scaled local correlation to be the ratio of the
local correlation $A_{kij}$ and a factor of $2t_{k}$.

######  Theorem 2.

Denote the limiting value of the scaled local correlation by
$\widetilde{A}_{kij}$,

$\displaystyle\widetilde{A}_{kij}=\lim_{t_{k}\rightarrow
0}\frac{A_{kij}}{2t_{k}}$ (18)

Then $\widetilde{A}_{kij}$ equals the inner product of the gradients of the
eigenfunctions $\phi_{i}$ and $\phi_{j}$ at $x_{k}$, that is,

$\displaystyle\widetilde{A}_{kij}=\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k}).$
(19)

Two proofs are provided in Appendix A and B. A brief summary is provided
below.

Proof 1. In the first proof we choose a sufficiently small $\epsilon_{k}$ and
show that

$\displaystyle\lim_{t_{k}\rightarrow 0}A_{kij}$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\int_{B_{\epsilon_{k}}(x_{k})}G(t_{k},x_{k},y)\Psi_{kij}(y)\omega_{g}(y)$
(20)

where $B_{\epsilon}(x)$ is defined in Eq. (6) and

$\displaystyle G(t,x,y)$ $\displaystyle=\frac{e^{-d_{g}(x,y)^{2}/4t}}{(4\pi
t)^{d/2}}.$ (21)

Then, by using the properties of the exponential map at $x_{k}$ and applying
basic techniques in calculus, we show that $\lim_{t_{k}\rightarrow
0}A_{kij}/2t_{k}$ evaluates to
$\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k})$.

Proof 2. In the second proof, as in [41, 42], we used the Feynman-Kac formula,

$\displaystyle A_{kij}$
$\displaystyle=[e^{-t_{k}\Delta_{g}}((\phi_{i}-\phi_{i}(x_{k}))(\phi_{j}-\phi_{j}(x_{k}))](x_{k})$
(22)

and note that

$\displaystyle\lim_{t_{k}\rightarrow
0}\frac{A_{kij}}{2t_{k}}=\left.\frac{1}{2}\frac{\partial A_{kij}}{\partial
t_{k}}\right|_{t_{k}=0}=\frac{-1}{2}\left\\{\Delta_{g}[(\phi_{i}-\phi_{i}(x_{k}))(\phi_{j}-\phi_{j}(x_{k}))](x_{k})\right\\}.$
(23)

Then, by applying the formula of the Laplacian of the product of two
functions, we show that the above equation equals
$\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k})$.

###  3.2 Estimate of $\widetilde{A}_{kij}$ in the Discrete Setting

To apply Theorems 1 and 2 in practice on data, we need an estimate of
$\widetilde{A}_{kij}$ in the discrete setting. There are several ways to
obtain this estimate. A generic way is by using the algorithms [9, 1] based on
Local Linear Regression (LLR) to estimate the gradient vector
$\nabla\phi_{i}(x_{k})$ itself from the values of $\phi_{i}$ in a neighbor of
$x_{k}$. An alternative approach is to use a finite sum approximation of Eq.
(20) combined with Eq. (18). A third approach is based on the Feynman-Kac
formula where we make use of Eq. (23) in the discrete setting. In the
following we explain the latter two approaches.

####  3.2.1 Finite sum approximation

Let $(x_{k})_{k=1}^{n}$ be uniformly distributed points on $(\mathcal{M},g)$.
Let $d_{e}(x_{k},x_{k}^{\prime})$ be the distance between $x_{k}$ and
$x_{k^{\prime}}$. The accuracy with which $\widetilde{A}_{kij}$ can be
estimated mainly depends on the accuracy of $d_{e}(\cdot\ ,\ \cdot)$ to the
local geodesic distances. For simplicity, we use $d_{e}(x_{k},x_{k}^{\prime})$
to be the Euclidean distance $\left\|x_{k}-x_{k^{\prime}}\right\|_{2}$. A more
accurate estimate of the local geodesic distances can be computed using the
method described in [29].

We construct a sparse unnormalized graph Laplacian $L$ using Algo. 1, where
the weight matrix $K$ of the graph edges is defined using the Gaussian kernel.
The bandwidth of the Gaussian kernel is set using the local scale of the
neighborhoods around each point as in self-tuning spectral clustering [47].
Let $\bm{\phi}_{i}$ be the $i$th non-trivial eigenvector of $L$ and denote
$\phi_{i}(x_{j})$ by $\bm{\phi}_{ij}$.

Input:
$d_{e}(x_{k},x_{k^{\prime}})_{k,k^{\prime}=1}^{n},k_{\textrm{nn}},k_{\textrm{tune}}\text{
where }k_{\textrm{tune}}\leq k_{\textrm{nn}}$

Output: $L$

1 $\mathcal{N}_{k}\leftarrow$ set of indices of $k_{\textrm{nn}}$ nearest
neighbours of $x_{k}$ based on $d_{e}(x_{k},\cdot)$;

2 $\sigma_{k}\leftarrow d_{e}(x_{k},x_{k^{*}})$ where $x_{k^{*}}$ is the
$k_{\textrm{tune}}$th nearest neighbor of $x_{k}$;

3 $K_{kk}\leftarrow 0,K_{kk^{\prime}}\leftarrow
e^{-d_{e}(x_{k},x_{k^{\prime}})^{2}/\sigma_{k}\sigma_{k^{\prime}}},k^{\prime}\in\mathcal{N}_{k}$;

4 $D_{kk}\leftarrow\sum_{k^{\prime}}K_{kk^{\prime}},\
D_{kk^{\prime}}\leftarrow 0,k\neq k^{\prime}$;

5 $L\leftarrow D-K$;

Algorithm 1 Sparse Unnormalized Graph Laplacian based on [47]

We estimate $\widetilde{A}_{kij}$ by evaluating the scaled local correlation
$A_{kij}/2t_{k}$ at a small value of $t_{k}$. The limiting value of $A_{kij}$
is estimated by substituting a small $t_{k}$ in the finite sum approximation
of the integral in Eq. (20). The sum is taken on a discrete ball of a small
radius $\epsilon_{k}$ around $x_{k}$ and is divided by $2t_{k}$ to obtain an
estimate of $\widetilde{A}_{kij}$.

We start by choosing $\epsilon_{k}$ to be the distance of $k_{\text{lv}}$th
nearest neighbor of $x_{k}$ where $k_{\text{lv}}$ is a hyperparameter with a
small integral value (subscript lv stands for local view). Thus,

$\displaystyle\epsilon_{k}=\text{distance to the }k_{\text{lv}}\text{th
nearest neighbor of }x_{k}.$ (24)

Then the limiting value of $t_{k}$ is given by

$\displaystyle\sqrt{\text{chi2inv}(p,d)}\sqrt{2t_{k}}=\epsilon_{k}\implies
t_{k}=\frac{1}{2}\frac{\epsilon_{k}^{2}}{\text{chi2inv}(p,d)},$ (25)

where chi2inv is the inverse cdf of the chi-squared distribution with $d$
degrees of freedom evaluated at $p$. We take $p$ to be $0.99$ in our
experiments. The rationale behind the above choice of $t_{k}$ is described in
Appendix C.

Now define the discrete ball around $x_{k}$ as

$\displaystyle U_{k}$ $\displaystyle=\\{x_{k^{\prime}}\ |\
d_{e}(x_{k},x_{k^{\prime}})\leq\epsilon_{k}\\}.$ (26)

Let $U_{k}$ denote the $k$th local view of the data in the high dimensional
ambient space. For convenience, denote the estimate of
$G(t_{k},x_{k},x_{k^{\prime}})$ by $G_{kk^{\prime}}$ where $G$ is as in Eq.
(21). Then

$\displaystyle G_{kk^{\prime}}$
$\displaystyle=\left\\{\begin{matrix}[l]\frac{\exp(-d_{e}(x_{k},x_{k^{\prime}})^{2}/4t_{k})}{\sum_{x\in
U_{k}}\exp(-d_{e}(x_{k},x)^{2}/4t_{k})}&,\ x_{k^{\prime}}\in
U_{k}-\\{x_{k}\\}\\\ 0&,\ \text{otherwise}.\end{matrix}\right.$ (27)

Finally, the estimate of $\widetilde{A}_{kij}$ is given by

$\displaystyle\widetilde{A}_{kij}$
$\displaystyle=\frac{1}{2t_{k}}G_{k}^{T}((\bm{\phi_{i}}-\bm{\phi_{ik}})\odot(\bm{\phi_{j}}-\bm{\phi_{jk}}))$
(28)

where $G_{k}$ is a column vector containing the $k$th row of the matrix $G$
and $\odot$ represents the Hadamard product.

####  3.2.2 Estimation based on Feynman-Kac formula

This approach to estimate $\widetilde{A}_{kij}$ is simply the discrete analog
of Eq. (23),

$\displaystyle\widetilde{A}_{kij}=\frac{-1}{2}L_{k}^{T}((\bm{\phi_{i}}-\bm{\phi_{ik}})\odot(\bm{\phi_{j}}-\bm{\phi_{jk}}))$
(29)

where $L_{k}$ is a column vector containing the $k$th row of $L$. A variant of
this approach which results in better estimates in the noisy case uses a low
rank approximation of $L$ using its first few eigenvectors (see Appendix H).

######  Remark 2.

It is not a coincidence that Eq. (28) and Eq. (29) look quite similar. In
fact, if we take $T$ to be a diagonal matrix with $(t_{k})_{k=1}^{n}$ as the
diagonal, then the matrix $T^{-1}(I-G)$ approximates $\Delta_{g}$ in the limit
of $(t_{k})_{k=1}^{n}$ tending to zero. Replacing $L$ with $T^{-1}(I-G)$ and
therefore $L_{k}$ with $(e_{k}-G_{k})/t_{k}$ reduces Eq. (29) to Eq. (28).
Here $e_{k}$ is a column vector with $k$th entry as $1$ and rest zeros.
Therefore the two approaches are the same in the limit.

######  Remark 3.

The above two approaches can also be generalized to compute the $\nabla
f_{i}(x_{k})^{T}\nabla f_{j}(x_{k})$ for arbitrary $\mathcal{C}^{2}$ mappings
$f_{i}$ and $f_{j}$ from $\mathcal{M}$ to $\mathbb{R}$ ( $\nabla
f_{i}(x_{k})=\nabla(f_{i}\ \circ\ \text{exp}_{x_{k}})(0)$ as per our
convention). To achieve this, simply replace $\phi_{i}$ and $\phi_{j}$ with
$f_{i}$ and $f_{j}$ in Eq. (28) and Eq. (29).

|
---
|
---
Figure 4: Comparison of different approaches to estimate $\widetilde{A}_{kij}$
in the discrete setting.

##### Example.

This example will follow us throughout the paper. Consider a square grid
$[0,1]\times[0,1]$ with a spacing of $0.01$ in both $x$ and $y$ direction.
With $k_{\textrm{nn}}=49$, $k_{\textrm{tune}}=7$ and
$d_{e}(x_{k},x_{k^{\prime}})=\left\|x_{k}-x_{k^{\prime}}\right\|_{2}$ as input
to the Algo. 1, we construct the graph Laplacian $L$. Using
$k_{\text{lv}}=25$, $d=2$ and $p=0.99$, we obtain the discrete balls $U_{k}$
and $t_{k}$. The $3$rd and $8$th eigenvectors of $L$ and the corresponding
analytical eigenfunctions are then obtained. The analytical value of
$\widetilde{A}_{k38}$ is displayed in Figure 4, followed by its estimate using
LLR [9], finite sum approximation and Feynman-Kac formula based approaches.
The analytical and the estimated values are normalized by
$\max_{k}\widetilde{A}_{kij}$ to bring them to the same scale. The absolute
error due to these approaches are shown below the estimates.

Even though, in this example, the Feynman-Kac formulation seem to have a
larger error, in our experiments, no single approach seem to be a clear winner
across all the examples. This becomes clear in Appendix H where we provided a
comparison of these approaches on a noiseless and a noisy Swiss Roll. The
results shown in this paper are based on finite sum approximation to estimate
$\widetilde{A}_{kij}$.

###  3.3 Low Distortion Local Parameterization from Laplacian Eigenvectors

We use $\nabla\phi_{i}\equiv\nabla\phi_{i}(x_{k})$ for brevity. Using the
estimates of $\widetilde{A}_{kij}$, we now present an algorithmic construction
of low distortion local parameterization $\Phi_{k}$ which maps $U_{k}$ into
$\mathbb{R}^{d}$. The pseudocode is provided below followed by a full
explanation of the steps and a note on the hyperparameters. Before moving
forward, it would be helpful for the reader to review the construction
procedure in the continuous setting in Section 2.1.

Input: $L,N,k_{\text{lv}},d,p,(\tau_{s},\delta_{s})_{s=1}^{d}$

Output: $(\Phi_{k},U_{k},\zeta_{kk})_{k=1}^{n}$

1 Compute $(\bm{\phi}_{i})_{i=1}^{N},\lambda_{1}\leq\ldots\leq\lambda_{N}$ by
eigendecomposition of $L$;

2 for $k\leftarrow 1$ to $n$ do

3 Compute $U_{k},(\widetilde{A}_{kij})_{i,j=1}^{N}$ (Eq. (26, 28));

4 Compute $(\gamma_{ki})_{i=1}^{N}$ (Eq. (30));

5 $\theta_{1}\leftarrow\tau_{1}$-percentile of
$(\widetilde{A}_{kii})_{i=1}^{N}$;

6 Compute $S_{k}$ (Eq. (31));

7 Compute $i_{1}$ (Eq. (35));

8 for $s\leftarrow 2$ to $d$ do

9 Compute $H^{s}_{kij}$ (Eq. (37));

10 $\theta_{s}\leftarrow\tau_{s}$-percentile of $(H^{s}_{kii})_{i\in S_{k}}$;

11 Compute $i_{s}$ (Eq. (42));

12

13 end for

14
$\Phi_{k}\leftarrow(\gamma_{ki_{1}}\bm{\phi}_{i_{1}},\ldots,\gamma_{ki_{d}}\bm{\phi}_{i_{d}})$
(Eq. (43));

15 Compute $\zeta_{kk}$ (Eq. (45));

16

17 end for

Algorithm 2 BiLipschitz-Local-Parameterization

An estimate of $\gamma_{ki}$ is obtained by the discrete analog of Eq. (7) and
is given by

$\displaystyle\gamma_{ki}=\text{Root-Mean-Square}(\\{\bm{\phi}_{ij}\ |\
x_{j}\in U_{k}\\})^{-1}.$ (30)

##### Step 1. Compute a set $S_{k}$ of candidate eigenvectors for $\Phi_{k}$.

Based on the construction procedure following Theorem 1, we start by computing
a set $S_{k}$ of candidate eigenvectors to construct $\Phi_{k}$ of $U_{k}$.
There is no easy way to retrieve the set $S_{k}$ in the discrete setting as in
the procedure. Therefore, we make the natural choice of using the first $N$
nontrivial eigenvectors $(\bm{\phi}_{i})_{i=1}^{N}$ of $L$ corresponding to
the $N$ smallest eigenvalues $(\lambda_{i})_{i=1}^{N}$, with sufficiently
large gradient at $x_{k}$, as the set $S_{k}$. The large gradient constraint
is required for the numerical stability of our algorithm. Therefore, we set
$S_{k}$ to be,

$\displaystyle S_{k}$ $\displaystyle=\\{i\in\\{1,\ldots,N\\}\ |\
\left\|\nabla\phi_{i}\right\|^{2}\geq\theta_{1}\\}=\\{i\in\\{1,\ldots,N\\}|\
\widetilde{A}_{kii}\geq\theta_{1}\\},$ (31)

where $\theta_{1}$ is $\tau_{1}$-percentile of the set
$(\widetilde{A}_{kii})_{i=1}^{N}$ and the second equality follows from Eq.
(19). Here $N$ and $\tau_{1}\in(0,100)$ are hyperparameters.

##### Step 2. Choose a direction $p_{1}\in T_{x_{k}}\mathcal{M}$.

The unit norm constraint on $p_{1}$ is relaxed. This will neither affect the
math nor the output of our algorithm. Since $p_{1}$ can be arbitrary we take
$p_{1}$ to be the gradient of an eigenvector $r_{1}$, that is
$\nabla\phi_{r_{1}}$. The choice of $r_{1}$ will determine
$\bm{\phi}_{i_{1}}$. To obtain a low frequency eigenvector, $r_{1}$ is chosen
so that the eigenvalue $\lambda_{r_{1}}$ is minimal, therefore

$\displaystyle r_{1}$ $\displaystyle=\mathop{\mathrm{argmin}}\limits_{j\in
S_{k}}\lambda_{j}.$ (32)

##### Step 3. Find $i_{1}\in S_{k}$ such that
$\gamma_{ki_{1}}|\nabla\phi_{i_{1}}^{T}p_{1}|$ is sufficiently large.

Since $p_{1}=\nabla\phi_{r_{1}}$, using Eq. (19), the formula for
$\nabla\phi_{i}^{T}p_{1}$ becomes

$\displaystyle\nabla\phi_{i}^{T}p_{1}$
$\displaystyle=\nabla\phi_{i}^{T}\nabla\phi_{r_{1}}=\widetilde{A}_{kir_{1}}.$
(33)

Then we obtain the eigenvector $\bm{\phi}_{i_{1}}$ so that
$\gamma_{ki_{1}}|\nabla\phi_{i_{1}}^{T}p_{1}|$ is larger than a certain
threshold. We do not know what the value of this threshold would be in the
discrete setting. Therefore, we first define the maximum possible value of
$\gamma_{ki_{1}}|\nabla\phi_{i}^{T}p_{1}|$ using Eq. (33) as

$\displaystyle\alpha_{1}=\underset{i\in S_{k}}{\max}\
\gamma_{ki}|\nabla\phi_{i}^{T}p_{1}|=\underset{i\in S_{k}}{\max}\
\gamma_{ki}|\widetilde{A}_{kir_{1}}|.$ (34)

Then we take the threshold to be $\delta_{1}\alpha_{1}$ where
$\delta_{1}\in(0,1]$ is a hyperparameter. Finally, to obtain a low frequency
eigenvector $\bm{\phi}_{i_{1}}$, we choose $i_{1}$ such that

$\displaystyle i_{1}$ $\displaystyle=\mathop{\mathrm{argmin}}\limits_{i\in
S_{k}}\\{\lambda_{i}:\gamma_{ki}|\nabla\phi_{i}^{T}p_{1}|\geq\delta_{1}\alpha_{1}\\}=\mathop{\mathrm{argmin}}\limits_{i\in
S_{k}}\\{\lambda_{i}:\gamma_{ki}|\widetilde{A}_{kir_{1}}|\geq\delta_{1}\alpha_{1}\\}.$
(35)

After obtaining $\bm{\phi}_{i_{1}}$, we use a recursive procedure to obtain
the $s$-th eigenvector $\bm{\phi}_{i_{s}}$ where $s\in\\{2,\ldots,d\\}$ in
order.

##### Step 4. Choose a direction $p_{s}\in T_{x_{k}}\mathcal{M}$ orthogonal to
$\nabla\phi_{i_{1}},\ldots,\nabla\phi_{i_{s}}$.

Again the unit norm constraint will be relaxed with no change in the output.
We are going to take $p_{s}$ to be the component of $\nabla\phi_{r_{s}}$
orthogonal to $\nabla\phi_{i_{1}},\ldots,\nabla\phi_{i_{s}}$ for a carefully
chosen $r_{s}$. For convenience, denote by $V_{s}$ the matrix with
$\nabla\phi_{i_{1}},\ldots,\nabla\phi_{i_{s-1}}$ as columns and let
$\mathcal{R}(V_{s})$ be the range of $V_{s}$. Let $\phi_{r_{s}}$ be an
eigenvector such that $\nabla\phi_{r_{s}}\not\in\mathcal{R}(V_{s})$. To find
such an $r_{s}$, we define

$\displaystyle H^{s}_{kij}$
$\displaystyle=\nabla\phi_{i}^{T}(I-V_{s}(V_{s}^{T}V_{s})^{-1}V_{s}^{T})\nabla\phi_{j}$
(36)
$\displaystyle=\widetilde{A}_{kij}-\begin{bmatrix}\widetilde{A}_{kii_{1}}\ldots\widetilde{A}_{kii_{s-1}}\end{bmatrix}\begin{bmatrix}\widetilde{A}_{ki_{1}i_{1}}&\widetilde{A}_{ki_{1}i_{2}}&\ldots&\widetilde{A}_{ki_{1}i_{s-1}}\\\
\widetilde{A}_{ki_{2}i_{1}}&\widetilde{A}_{ki_{2}i_{2}}&\ldots&\widetilde{A}_{ki_{2}i_{s-1}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\widetilde{A}_{ki_{s-1}i_{1}}&\widetilde{A}_{ki_{s-1}i_{2}}&\ldots&\widetilde{A}_{ki_{s-1}i_{s-1}}\end{bmatrix}^{-1}\begin{bmatrix}\widetilde{A}_{ki_{1}j}\\\
\widetilde{A}_{ki_{2}j}\\\ \vdots\\\ \widetilde{A}_{ki_{s-1}j}\end{bmatrix}$
(37)

Note that $H^{s}_{kii}$ is the squared norm of the projection of
$\nabla\phi_{i}$ onto the vector space orthogonal to $\mathcal{R}(V_{s})$.
Clearly $\nabla\phi_{i}\not\in\mathcal{R}(V_{s})$ if and only if
$H^{s}_{kii}>0$. To obtain a low frequency eigenvector $\bm{\phi}_{r_{s}}$
such that $H^{s}_{kr_{s}r_{s}}>0$ we choose

$\displaystyle r_{s}=\mathop{\mathrm{argmin}}\limits_{i\in
S_{k}}\\{\lambda_{i}:H^{s}_{kii}\geq\theta_{s}\\}$ (38)

where $\theta_{s}$ is the $\tau_{s}$-percentile of the set
$\\{H^{s}_{kii}:i\in S_{k}\\}$ and $\tau_{s}\in(0,100)$ is a hyperparameter.
Then we take $p_{s}$ to be the component of $\nabla\phi_{r_{s}}$ which is
orthogonal to $\mathcal{R}(V_{s})$,

$\displaystyle
p_{s}=(I-V_{s}(V_{s}^{T}V_{s})^{-1}V_{s}^{T})\nabla\phi_{r_{s}}.$ (39)

##### Step 5. Find $i_{s}\in S_{k}$ such that
$\gamma_{ki_{s}}|\nabla\phi_{i_{s}}^{T}p_{s}|$ is sufficiently large.

Using Eq. (36, 39), we note that

$\displaystyle\nabla\phi_{i}^{T}p_{s}=H^{s}_{kir_{s}}.$ (40)

To obtain $\bm{\phi}_{i_{s}}$ such that
$\gamma_{ki_{s}}|\nabla\phi_{i_{s}}^{T}p_{s}|$ is greater than a certain
threshold, as in step $3$, we first define the maximum possible value of
$\gamma_{ki_{s}}|\nabla\phi_{i}^{T}p_{s}|$ using Eq. (40) as,

$\displaystyle\alpha_{s}=\max_{i\in
S_{k}}\gamma_{ki}|\nabla\phi_{i}^{T}p_{s}|=\max_{i\in
S_{k}}\gamma_{ki}|H^{s}_{kir_{s}}|.$ (41)

Then we take the threshold to be $\delta_{s}\alpha_{s}$ where
$\delta_{s}\in[0,1]$ is a hyperparameter. Finally, to obtain a low frequency
eigenvector $\bm{\phi}_{i_{s}}$ we choose $i_{s}$ such that

$\displaystyle i_{s}$ $\displaystyle=\mathop{\mathrm{argmin}}\limits_{i\in
S_{k}}\\{\lambda_{i}:\gamma_{ki}|\nabla\phi_{i}^{T}p_{s}|\geq\delta_{s}\alpha_{s}\\}=\mathop{\mathrm{argmin}}\limits_{i\in
S_{k}}\\{\lambda_{i}:\gamma_{ki}|H^{s}_{kir_{s}}|\geq\delta_{s}\alpha_{s}\\}.$
(42)

In the end we obtain a $d$-dimensional parameterization $\Phi_{k}$ of $U_{k}$
given by

$\displaystyle\Phi_{k}$
$\displaystyle\equiv(\gamma_{ki_{1}}\bm{\phi}_{i_{1}},\ldots,\gamma_{ki_{d}}\bm{\phi}_{i_{d}})\
\text{where}$ $\displaystyle\Phi_{k}(x_{k^{\prime}})$
$\displaystyle=(\gamma_{k{i_{1}}}\bm{\phi}_{i_{1}k^{\prime}},\ldots,\gamma_{k{i_{d}}}\bm{\phi}_{i_{d}k^{\prime}})\
\text{and}$ (43) $\displaystyle\Phi_{k}(U_{k})$
$\displaystyle=(\Phi_{k}(x_{k^{\prime}}))_{x_{k^{\prime}}\in U_{k}}.$

We call $\Phi_{k}(U_{k})$ the $k$th local view of the data in the
$d$-dimensonal embedding space. It is a matrix with $|U_{k}|$ rows and $d$
columns. Denote the distortion of $\Phi_{k^{\prime}}$ on $U_{k}$ by
$\zeta_{kk^{\prime}}$. Using Eq. (1) we obtain

$\displaystyle\zeta_{kk^{\prime}}$
$\displaystyle=\text{Distortion}(\Phi_{k^{\prime}},U_{k})$ (44)
$\displaystyle=\sup_{\begin{subarray}{c}x_{l},x_{l^{\prime}}\in U_{k}\\\
x_{l}\neq
x_{l^{\prime}}\end{subarray}}\frac{\left\|\Phi_{k^{\prime}}(x_{l})-\Phi_{k^{\prime}}(x_{l^{\prime}})\right\|}{d_{e}(x_{l},x_{l^{\prime}})}\sup_{\begin{subarray}{c}x_{l},x_{l^{\prime}}\in
U_{k}\\\ x_{l}\neq
x_{l^{\prime}}\end{subarray}}\frac{d_{e}(x_{l},x_{l^{\prime}})}{\left\|\Phi_{k^{\prime}}(x_{l})-\Phi_{k^{\prime}}(x_{l^{\prime}})\right\|}.$
(45)

##### Postprocessing.

The obtained local parameterizations are post-processed so as to remove the
anomalous parameterizations having unusually high distortion. We replace the
local parameterization $\Phi_{k}$ of $U_{k}$ by that of a neighbor,
$\Phi_{k^{\prime}}$ where $x_{k^{\prime}}\in U_{k}$, if the distortion
$\zeta_{kk^{\prime}}$ produced by $\Phi_{k^{\prime}}$ on $U_{k}$ is smaller
than the distortion $\zeta_{kk}$ produced by $\Phi_{k}$ on $U_{k}$. If
$\zeta_{kk^{\prime}}<\zeta_{kk}$ for multiple $k^{\prime}$ then we choose the
parameterization which produces the least distortion on $U_{k}$. This
procedure is repeated until no replacement is possible. The pseudocode is
provided below.

Input:
$d_{e}(x_{k},x_{k^{\prime}})_{k,k^{\prime}=1}^{n},(I_{k},\Phi_{k},\zeta_{kk})_{k=1}^{n}$

Output: $(\Phi_{k},\zeta_{kk})_{k=1}^{n}$

1 $N_{\text{replaced}}\leftarrow 1$;

2 while $N_{\text{replaced}}>0$ do

3 $N_{\text{replaced}}\leftarrow 0$;

4 $\Phi^{\text{old}}_{k}\leftarrow\Phi_{k}$ for all $k\in\\{1,\ldots,n\\}$;

5 for $k\leftarrow 1$ to $n$ do

6 Compute $(\zeta_{kk^{\prime}})_{x_{k^{\prime}}\in U_{k}}$ (Eq. (45));

7 $k^{*}\leftarrow\mathop{\mathrm{argmin}}\limits_{x_{k^{\prime}}\in
U_{k}}\zeta_{kk^{\prime}}$;

8 if $k^{*}\neq k$ then

9 $\ \Phi_{k}\leftarrow\Phi_{k^{*}}^{\text{old}};\ \
\zeta_{kk}\leftarrow\zeta_{kk^{*}};\ \ N_{\text{replaced}}\leftarrow
N_{\text{replaced}}+1$;

10

11 end if

12

13 end for

14

15 end while

Algorithm 3 Postprocess-Local-Parameterization

A note on hyperparameters $N,(\tau_{s},\delta_{s})_{s=1}^{d}$. Generally, $N$
should be small so that the low frequency eigenvectors form the set of
candidate eigenvectors. In almost all of our experiments we take $N$ to be
$100$. The set of $(\tau_{s},\delta_{s})_{s=1}^{d}$ is reduced to two
hyperparameters, one for all $\tau_{s}$’s and one for all $\delta_{s}$’s. As
explained above, $\tau_{s}$ enforces certain vectors to be non-zero and
$\delta_{s}$ enforces certain directional derivatives to be large enough.
Therefore, a small value of $\tau_{s}$ in $(0,100)$ and a large value of
$\delta_{s}$ in $(0,1]$ is suitable. In most of our experiments, we used a
value of $50$ for all $\tau_{s}$ and a value of $0.9$ for all $\delta_{s}$.
Our algorithm is not too sensitive to the values of these hyperparameters.
Other values of $N$, $\tau_{s}$ and $\delta_{s}$ would also result in the
embeddings with high visual quality.

##### Example.

We now build upon the example of the square grid at the end of Section 3.2.
The values of the additional inputs are $N=100$, $\tau_{s}=50$ and
$\delta_{s}=0.9$ for all $s\in\\{1,\ldots,d\\}$. Using Algo. 2 and 3 we obtain
$10^{4}$ local views $U_{k}$ and $\Phi_{k}(U_{k})$ where $|U_{k}|=25$ for all
$k$. In the left image of Figure 5, we colored each point $x_{k}$ with the
distortion $\zeta_{kk}$ of the local parameterization $\Phi_{k}$ on $U_{k}$.
The mapped discrete balls $\Phi_{k}(U_{k})$ for some values of $k$ are also
shown in Figure 30 in the Appendix H.

|
---|---
Figure 5: Distortion of the obtained local parameterizations when the points
on the boundary are not known (left) versus when they are known apriori
(right). Each point $x_{k}$ is colored by $\zeta_{kk}$ (see Eq. (45)).

######  Remark 4.

Note that the parameterizations of the discrete balls close to the boundary
have higher distortion. This is because the injectivity radius at the points
close to the boundary is low and precisely zero at the points on the boundary.
As a result, the size of the balls around these points exceeds the limit
beyond which Theorem 1 is applicable.

At this point we note the following remark in [25].

######  Remark 5.

As was noted by L. Guibas, when M has a boundary, in the case of Neumann
boundary values, one may consider the “doubled” manifold, and may apply the
result in Theorem 1 for a possibly larger $r_{k}$.

Due to the above remark, assuming that the points on the boundary are known,
we computed the distance matrix for the doubled manifold using the method
described in [27]. Then we recomputed the local parameterizations $\Phi_{k}$
keeping all other hyperparameters the same as before. In the right image of
Figure 5, we colored each point $x_{k}$ with the distortion of the updated
parameterization $\Phi_{k}$ on $U_{k}$. Note the reduction in the distortion
of the paramaterizations for the neighborhoods close to the boundary. The
distortion is still high near the corners.

###  3.4 Time Complexity

The combined worst case time complexity of Algo. 1, 2 and 3 is
$O(n(N^{2}(k_{\text{lv}}+d)+k_{\text{lv}}^{3}N_{\text{post}}d))$ where
$N_{\text{post}}$ is the number of iterations it takes to converge in Algo. 3
which was observed to be less than $50$ for all the examples in this paper. It
took about a minute111Machine specification: MacOS version 11.4, Apple M1
Chip, $16$GB RAM. to construct the local views in the above example as well as
in all the examples in Section 6.

##  4 Clustering for Intermediate Views

Recall that the discrete balls $U_{k}$ are the local views of the data in the
high dimensional ambient space. In the previous section, we obtained the
mappings $\Phi_{k}$ to construct the local views $\Phi_{k}(U_{k})$ of the data
in the $d$-dimensional embedding space. As discussed in Section 1.2, one can
use the GPA [18, 20, 43] to register these local views to recover a global
embedding. In practice, too many small local views (high $n$ and small
$|U_{k}|$) result in extremely high computational complexity. Moreover, small
overlaps between the local views makes their registration susceptible to
errors. Therefore, we perform clustering to obtain $M\ll n$ intermediate
views, $\widetilde{U}_{m}$ and $\widetilde{\Phi}_{m}(\widetilde{U}_{m})$, of
the data in the ambient space and the embedding space, respectively. This
reduces the time complexity and increases the overlaps between the views,
leading to their quick and robust registration.

###  4.1 Notation

Our clustering algorithm is designed so as to ensure low distortion of the
parameterizations $\widetilde{\Phi}_{m}$ on $\widetilde{U}_{m}$. We first
describe the notation used and then present the pseudocode followed by a full
explanation of the steps. Let $c_{k}$ be the index of the cluster $x_{k}$
belongs to. Then the set of points which belong to cluster $m$ is given by

$\displaystyle\mathcal{C}_{m}=\\{x_{k}\ |\ c_{k}=m\\}.$ (46)

Denote by $c_{U_{k}}$ the set of indices of the neighboring clusters of
$x_{k}$. The neighboring points of $x_{k}$ lie in these clusters, that is,

$\displaystyle c_{U_{k}}=\\{c_{k^{\prime}}\ |\ x_{k^{\prime}}\in U_{k}\\}.$
(47)

We say that a point $x_{k}$ lies in the vicinity of a cluster $m$ if $m\in
c_{U_{k}}$. Let $\widetilde{U}_{m}$ denote the $m$th intermediate view of the
data in the ambient space. This constitutes the union of the local views
associated with all the points belonging to cluster $m$, that is,

$\displaystyle\widetilde{U}_{m}=\bigcup_{k:\ x_{k}\in\mathcal{C}_{m}}U_{k}.$
(48)

Clearly, a larger cluster means a larger intermediate view. In particular,
addition of $x_{k}$ to $\mathcal{C}_{m}$ grows the intermediate view
$\widetilde{U}_{m}$ to $\widetilde{U}_{m}\cup U_{k}$,

$\displaystyle\mathcal{C}_{m}\rightarrow\mathcal{C}_{m}\cup\\{x_{k}\\}\implies\widetilde{U}_{m}\rightarrow\widetilde{U}_{m}\cup
U_{k}$ (49)

Let $\widetilde{\Phi}_{m}$ be the $d$-dimensional parameterization associated
with the $m$th cluster. This parameterization maps $\widetilde{U}_{m}$ to
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$, the $m$th intermediate view of the
data in the embedding space. Note that a point $x_{k}$ generates the local
view $U_{k}$ (see Eq. (26)) which acts as the domain of the parameterization
$\Phi_{k}$. Similarly, a cluster $\mathcal{C}_{m}$ obtained through our
procedure, generates an intermediate view $\widetilde{U}_{m}$ (see Eq. (48))
which acts as the domain of the parameterization $\widetilde{\Phi}_{m}$.
Overall, our clustering procedure replaces the notion of a local view per an
individual point by an intermediate view per a cluster of points.

Input: $(U_{k},\Phi_{k})_{k=1}^{n},\eta_{\text{min}}$

Output:
$(\mathcal{C}_{m},\widetilde{U}_{m},\widetilde{\Phi}_{m})_{m=1}^{M},(c_{k})_{k=1}^{n}$

1 Initialize $c_{k}\leftarrow k$, $\mathcal{C}_{m}\leftarrow\\{x_{m}\\}$,
$\widetilde{\Phi}_{m}\leftarrow\Phi_{m}$ for all $k,m\in\\{1,\ldots,n\\}$;

2 for $\eta\leftarrow 2$ to $\eta_{\text{min}}$ do

3 Compute $b_{m\leftarrow x_{k}}$ for all $m,k\in\\{1,\ldots,n\\}$ (Eq. (47,
48, 50));

4
$m,k\leftarrow\mathop{\mathrm{argmax}}\limits_{m^{\prime},k^{\prime}}b_{m^{\prime}\leftarrow
x_{k^{\prime}}};\ \text{bid}^{*}\leftarrow b_{m\leftarrow x_{k}}$;

5 while $\text{bid}^{*}>0$ do

6 $s\leftarrow c_{k};\ \mathcal{C}_{s}\leftarrow\mathcal{C}_{s}-x_{k};\
c_{k}\leftarrow m;\ \mathcal{C}_{m}\leftarrow\mathcal{C}_{m}\cup x_{k}$;

7 Recompute $b_{m^{\prime}\leftarrow x_{k^{\prime}}}$ for all
$(m^{\prime},k^{\prime})\in\mathcal{S}$ (Eq. (51));

8
$m,k\leftarrow\mathop{\mathrm{argmax}}\limits_{m^{\prime},k^{\prime}}b_{m^{\prime}\leftarrow
x_{k^{\prime}}};\ \text{bid}^{*}\leftarrow b_{m\leftarrow x_{k}}$;

9

10 end while

11

12 end for

13$M\leftarrow$ the number of non-empty clusters;

14 Remove $\mathcal{C}_{m}$, $\widetilde{\Phi}_{m}$ when
$|\mathcal{C}_{m}|=0$, relabel clusters from $1$ to $M$ and update $c_{k}$
with new labels;

15 Compute $(\widetilde{U}_{m})_{m=1}^{M}$ (Eq. (48));

Algorithm 4 Clustering

###  4.2 Low Distortion Clustering

Initially, we start with $n$ singleton clusters where the point $x_{k}$
belongs to the $k$th cluster and the parameterization associated with the
$k$th cluster is $\Phi_{k}$. Thus, $c_{k}=k$, $\mathcal{C}_{m}=\\{x_{m}\\}$
and $\widetilde{\Phi}_{m}=\Phi_{m}$ for all $k,m\in\\{1,\ldots,n\\}$. This
automatically implies that initially $\widetilde{U}_{m}=U_{m}$. The
parameterizations associated with the clusters remain the same throughout the
procedure. During the procedure, each cluster $\mathcal{C}_{m}$ is perceived
as an entity which wants to grow the domain $\widetilde{U}_{m}$ of the
associated parameterization $\widetilde{\Phi}_{m}$ by growing itself (see Eq.
49), while simultaneously keeping the distortion of $\widetilde{\Phi}_{m}$ on
$\widetilde{U}_{m}$ low (see Eq. 45). To achieve that, each cluster
$\mathcal{C}_{m}$ places a careful bid $b_{m\leftarrow x_{k}}$ for each point
$x_{k}$. The global maximum bid is identified and the underlying point $x_{k}$
is relabelled to the bidding cluster, hence updating $c_{k}$. With this
relabelling, the bidding cluster grows and the source cluster shrinks. This
procedure of shrinking and growing clusters is repeated until all non-empty
clusters are large enough, i.e. have a size at least $\eta_{\text{min}}$, a
hyperparameter. In our experiments, we choose $\eta_{\text{min}}$ from
$\\{5,10,15,20,25\\}$. We iterate over $\eta$ which varies from $2$ to
$\eta_{\text{min}}$. In the $\eta$-th iteration, we say that the $m$th cluster
is small if it is non-empty and has a size less than $\eta$, that is, when
$|\mathcal{C}_{m}|\in(0,\eta)$. During the iteration, the clusters either
shrink or grow until no small clusters remain. Therefore, at the end of the
$\eta$-th iteration the non-empty clusters are of size at least $\eta$. After
the last ($\eta_{\text{min}}$th) iteration, each non-empty cluster will have
at least $\eta_{min}$ points and the empty clusters are pruned away.

##### Bid by cluster $m$ for $x_{k}$.

In the $\eta$-th iteration, we start by computing the bid $b_{m\leftarrow
x_{k}}$ by each cluster $m$ for each point $x_{k}$. The bid function is
designed so as to satisfy the following conditions. The first two conditions
are there to halt the procedure while the last two conditions follow
naturally. These conditions are also depicted in Figure 6.

  1. 1.

No cluster bids for the points in large clusters. Since $x_{k}$ belongs to
cluster $c_{k}$ therefore, if $|\mathcal{C}_{c_{k}}|>\eta$ then the
$b_{m\leftarrow x_{k}}$ is zero for all $m$.

  2. 2.

No cluster bids for a point in another cluster whose size is bigger than its
own size. Therefore, if $|\mathcal{C}_{m}|<|\mathcal{C}_{c_{k}}|$ then again
$b_{m\leftarrow x_{k}}$ is zero.

  3. 3.

A cluster bids for the points in its own vicinity. Therefore, if $m\not\in
c_{U_{k}}$ (see Eq. 47) then $b_{m\leftarrow x_{k}}$ is zero.

  4. 4.

Recall that a cluster $m$ aims to grow while keeping the distortion of
associated parameterization $\widetilde{\Phi}_{m}$ low on its domain
$\widetilde{U}_{m}$. If the $m$th cluster acquires the point $x_{k}$,
$\widetilde{U}_{m}$ grows due to the addition of $U_{k}$ to it (see Eq. (48)),
and so does the distortion of $\widetilde{\Phi}_{m}$ on it. Therefore, to
ensure low distortion, the natural bid by $\mathcal{C}_{m}$ for the point
$x_{k}$, $b_{m\leftarrow x_{k}}$, is
$\text{Distortion}(\widetilde{\Phi}_{m},U_{k}\cup\widetilde{U}_{m})^{-1}$ (see
Eq. 45).

Combining the above conditions, we can write the bid by cluster $m$ for the
point $x_{k}$ as,

$\displaystyle b_{m\leftarrow
x_{k}}=\left\\{\begin{matrix}\text{Distortion}(\widetilde{\Phi}_{m},U_{k}\cup\widetilde{U}_{m})^{-1}&\text{if
}|\mathcal{C}_{c_{k}}|\in(0,\eta)\land m\in
c_{U_{k}}\land|\mathcal{C}_{m}|\geq|\mathcal{C}_{c_{k}}|\\\
0&\text{otherwise}.\end{matrix}\right.$ (50)

In the practical implementation of above equation, $c_{U_{k}}$ and
$\widetilde{U}_{m}$ are computed on the fly using Eq. (47, 48).

Figure 6: Computation of the bid for a point in a small cluster by the
neighboring clusters in $\eta$-th iteration. (left) $x_{k}$ is a point
represented by a small red disc, in a small cluster $c_{k}$ enclosed by solid
red line. The dashed red line enclose $U_{k}$. Assume that the cluster $c_{k}$
is small so that $|\mathcal{C}_{c_{k}}|\in(0,\eta)$. Clusters $m_{1}$,
$m_{2}$, $m_{3}$ and $m_{4}$ are enclosed by solid colored lines too. Note
that $m_{1}$, $m_{2}$ and $m_{3}$ lie in $c_{U_{k}}$ (the nonempty overlap
between these clusters and $U_{k}$ indicate that), while $m_{4}\not\in
c_{U_{k}}$. Thus, the bid by $m_{4}$ for $x_{k}$ is zero. Since the size of
cluster $m_{3}$ is less than the size of cluster $c_{k}$ i.e.
$|\mathcal{C}_{m_{3}}|<|\mathcal{C}_{c_{k}}|$, the bid by $m_{3}$ for $x_{k}$
is also zero. Since clusters $m_{1}$ and $m_{2}$ satisfy all the conditions,
the bids by $m_{1}$ and $m_{2}$ for $x_{k}$ are to be computed. (right) The
bid $b_{m_{1}\leftarrow x_{k}}$, is given by the inverse of the distortion of
$\widetilde{\Phi}_{m_{1}}$ on $U_{k}\cup\widetilde{U}_{m_{1}}$, where the
dashed blue line enclose $\widetilde{U}_{m_{1}}$. If the bid
$b_{m_{1}\leftarrow x_{k}}$ is greater (less) than the bid $b_{m_{2}\leftarrow
x_{k}}$, then the clustering procedure would favor relabelling of $x_{k}$ to
$m_{1}$ ($m_{2}$).

##### Greedy procedure to grow and shrink clusters.

Given the bids by all the clusters for all the points, we grow and shrink the
clusters so that at the end of the current iteration $\eta$, each non-empty
cluster has a size at least $\eta$. We start by picking the global maximum
bid, say $b_{m\leftarrow x_{k}}$. Let $x_{k}$ be in the cluster $s$ (note that
$c_{k}$, the cluster of $x_{k}$, is $s$ before $x_{k}$ is relabelled). We
relabel $c_{k}$ to $m$, and update the set of points in clusters $s$ and $m$,
$\mathcal{C}_{s}$ and $\mathcal{C}_{m}$, using Eq. (46). This implicitly
shrinks $\widetilde{U}_{s}$ and grows $\widetilde{U}_{m}$ (see Eq. 48) and
affects the bids by clusters $m$ and $s$ or the bids for the points in these
clusters. Denote the set of pairs of the indices of all such clusters and the
points by

$\displaystyle\mathcal{S}=\\{(m^{\prime},k^{\prime})\in\\{1,\ldots,n\\}^{2}\
|\ m^{\prime}\in\\{m,s\\}\text{ or
}x_{k^{\prime}}\in\mathcal{C}_{s}\cup\mathcal{C}_{m}\\}.$ (51)

Then the bids $b_{m^{\prime}\leftarrow x_{k^{\prime}}}$ are recomputed for all
$(m^{\prime},k^{\prime})\in\mathcal{S}$. It is easy to verify that for all
other pairs, neither the conditions nor the distortion in Eq. (50) are
affected. After this computation, we again pick the global maximum bid and
repeat the procedure until the maximum bid becomes zero indicating that no
non-empty small cluster remains. This marks the end of the $\eta$-th
iteration.

##### Final intermediate views in the ambient and the embedding space.

At the end of the last iteration, all non-empty clusters have at least
$\eta_{\text{min}}$ points. Let $M$ be the number of non-empty clusters. Using
the pigeonhole principle one can show that $M$ would be less than or equal to
$n/\eta_{\text{min}}$. We prune away the empty clusters and relabel the non-
empty ones from $1$ to $M$ while updating $c_{k}$ accordingly. With this, we
obtain the clusters $(\mathcal{C}_{m})_{m=1}^{M}$ with associated
parameterizations $(\widetilde{\Phi}_{m})_{m=1}^{M}$. Finally, using Eq. (48),
we obtain the $M$ intermediate views $(\widetilde{U}_{m})_{m=1}^{M}$ of the
data in the ambient space. Then, the intermediate views of the data in the
embedding space are given by
$(\widetilde{\Phi}_{m}(\widetilde{U}_{m}))_{m=1}^{M}$. Note that
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ is a matrix with
$|\widetilde{U}_{m}|$ rows and $d$ columns (see Eq. (43)).

##### Example.

We continue with our example of the square grid which originally contained
about $10^{4}$ points. Therefore, before clustering we had about $10^{4}$
small local views $U_{k}$ and $\Phi_{k}(U_{k})$, each containing $25$ points.
After clustering with $\eta_{min}=10$, we obtained $635$ clusters and
therefore that many intermediate views $\widetilde{U}_{m}$ and
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ with an average size of $79$. When
the points on the boundary are known then we obtained $562$ intermediate views
with an average size of $90$. Note that there is a trade-off between the size
of the intermediate views and the distortion of the parameterizations used to
obtain them. For convenience, define $\tilde{\zeta}_{mm}$ to be the distortion
of $\widetilde{\Phi}_{m}$ on $\widetilde{U}_{m}$ using Eq. (45). Then, as the
size of the views are increased (by increasing $\eta_{min}$), the value of
$\tilde{\zeta}_{mm}$ would also increase. In Figure 7 we colored the points in
cluster $m$, $\mathcal{C}_{m}$, with $\tilde{\zeta}_{mm}$. In other words,
$x_{k}$ is colored by $\tilde{\zeta}_{c_{k}c_{k}}$. Note the increased
distortion in comparison to Figure 5.

|
---|---
Figure 7: Each point $x_{k}$ colored by $\tilde{\zeta}_{c_{k}c_{k}}$ when the
points on the boundary of the square grid are unknown (left) versus when they
are known apriori (right).

###  4.3 Time Complexity

Our practical implementation of Algo. 4 uses memoization for speed up. It took
about a minute to construct intermediate views using in the above example with
$n=10^{4}$, $k_{\text{lv}}=25$, $d=2$ and $\eta_{\text{min}}=10$, and it took
less than $2$ minutes for all the examples in Section 6. It was empirically
observed that the time for clustering is linear in $n$, $\eta_{\text{min}}$
and $d$ while it is cubic in $k_{\text{lv}}$.

##  5 Global Embedding using Procrustes Analysis

In this section, we present an algorithm based on Procrustes analysis to align
the intermediate views $\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ and obtain a
global embedding. The $M$ views $\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ are
transformed by an orthogonal matrix $T_{m}$ of size $d\times d$, a
$d$-dimensional translation vector $v_{m}$ and a positive scalar $b_{m}$ as a
scaling component. The transformed views are given by
$\widetilde{\Phi}^{g}_{m}(\widetilde{U}_{m})$ such that

$\displaystyle\widetilde{\Phi}^{g}_{m}(x_{k})=b_{m}\widetilde{\Phi}_{m}(x_{k})T_{m}+v_{m}\quad\textrm{for
all }x_{k}\in\widetilde{U}_{m}.$ (52)

First we state a general approach to estimate these parameters, and its
limitations in Section 5.1. Then we present an algorithm in Section 5.2 which
computes these parameters and a global embedding of the data while addressing
the limitations of the general procedure. In Section 5.3 we describe a simple
modification to our algorithm to tear apart closed manifolds. In Appendix F,
we contrast our global alignment procedure with that of LTSA.

---
Figure 8: (left) The intermediate views $\widetilde{U}_{m}$ and
$\widetilde{U}_{m^{\prime}}$ of a $2$d manifold in a possibly high dimensional
ambient space. These views trivially align with each other. The red star in
blue circles represent their overlap $\widetilde{U}_{mm^{\prime}}$. (middle)
The $m$th and $m^{\prime}$th intermediate views in the $2$d embedding space.
(right) Transformed views after aligning
$\widetilde{\Phi}_{m}(\widetilde{U}_{mm^{\prime}})$ with
$\widetilde{\Phi}_{m^{\prime}}(\widetilde{U}_{mm^{\prime}})$.

###  5.1 General Approach for Alignment

In general, the parameters $(T_{m},v_{m},b_{m})_{m=1}^{M}$ are estimated so
that for all $m$ and $m^{\prime}$, the two transformed views of the overlap
between $\widetilde{U}_{m}$ and $\widetilde{U}_{m^{\prime}}$, obtained using
the parameterizations $\widetilde{\Phi}^{g}_{m}$ and
$\widetilde{\Phi}^{g}_{m^{\prime}}$, align with each other. To be more
precise, define the overlap between the $m$th and the $m^{\prime}$th
intermediate views in the ambient space as the set of points which lie in both
the views,

$\displaystyle\widetilde{U}_{mm^{\prime}}=\widetilde{U}_{m}\cap\widetilde{U}_{m^{\prime}}.$
(53)

In the ambient space, the $m$th and the $m^{\prime}$th views are neighbors if
$\widetilde{U}_{mm^{\prime}}$ is non-empty. As shown in Figure 8 (left), these
neighboring views trivially align on the overlap between them. It is natural
to ask for a low distortion global embedding of the data. Therefore, we must
ensure that the embeddings of $\widetilde{U}_{mm^{\prime}}$ due to the $m$th
and the $m^{\prime}$th view in the embedding space, also align with each
other. Thus, the parameters $(T_{m},v_{m},b_{m})_{m=1}^{M}$ are estimated so
that $\widetilde{\Phi}^{g}_{m}(\widetilde{U}_{mm^{\prime}})$ aligns with
$\widetilde{\Phi}^{g}_{m^{\prime}}(\widetilde{U}_{mm^{\prime}})$ for all $m$
and $m^{\prime}$. However, due to the distortion of the parameterizations it
is usually not possible to perfectly align the two embeddings (see Figure 8).
We can represent both embeddings of the overlap as matrices with
$|\widetilde{U}_{mm^{\prime}}|$ rows and $d$ columns. Then we choose the
measure of the alignment error to be the squared Frobenius norm of the
difference of the two matrices. The error is trivially zero if
$\widetilde{U}_{mm^{\prime}}$ is empty. Overall, the parameters are estimated
so as to minimize the following alignment error

$\displaystyle\mathcal{L}((T_{m},v_{m},b_{m})_{m=1}^{M})=\frac{1}{2M}\sum_{\begin{subarray}{c}m=1\\\
m^{\prime}=1\end{subarray}}^{M}\left\|\widetilde{\Phi}^{g}_{m}(\widetilde{U}_{mm^{\prime}})-\widetilde{\Phi}^{g}_{m^{\prime}}(\widetilde{U}_{mm^{\prime}})\right\|^{2}_{F}.$
(54)

In theory, one can start with a trivial initialization of $T_{m}$, $v_{m}$ and
$b_{m}$ as $I_{d}$, $\mathbf{0}$ and $1$, and directly use GPA [18, 20, 43] to
obtain a local minimum of the above alignment error. This approach has two
issues.

  1. 1.

Like most optimization algorithms, the rate of convergence to a local minimum
and the quality of it depends on the initialization of the parameters. We
empirically observed that with a trivial initialization of the parameters, GPA
may take a great amount of time to converge and may also converge to an
inferior local minimum.

  2. 2.

Using GPA to align a view with all of its adjacent views would prevent us from
tearing apart closed manifolds; as an example see Figure 11.

These issues are addressed in subsequent Sections 5.2 and 5.3, respectively.

###  5.2 GPA Adaptation for Global Alignment

Input:
$(x_{k},c_{k},w)_{k=1}^{n},(\mathcal{C}_{m},\widetilde{\Phi}_{m},\widetilde{U}_{m})_{m=1}^{M}$,
to_tear, $\nu$, $N_{r}$

Output: $(T_{m},b_{m},v_{m})_{m=1}^{M}$

1 for $\text{Iter}\leftarrow 1$ to $N_{r}+1$ do

2 if $\text{Iter}=1$ then

3 Initialize $T_{m}\leftarrow I,v_{m}\leftarrow 0$;

4 Compute $b_{m}$ (Eq. (55));

5 Compute $(s_{m},p_{s_{m}})_{m=1}^{M}$ (Eq. (98, 100) in Appendix D);

6 $\mathcal{A}\leftarrow\\{s_{1}\\}$ %The set of already transformed views;

7

8 else

9 $(s_{m})_{m=2}^{M}\leftarrow$ random permutation of $(1,\ldots,M)$ excluding
$s_{1}$;

10

11 end if

12 for $m\leftarrow 2$ to $M$ do

13 $s\leftarrow s_{m}$, $p\leftarrow p_{s_{m}}$;

14 (Step R1) $T_{s},v_{s}\leftarrow$ Procrustes
($\widetilde{\Phi}^{g}_{p}(\widetilde{U}_{sp})$,$\widetilde{\Phi}^{g}_{s}(\widetilde{U}_{sp})$,
No scaling);

15 if $\text{to\\_tear}=$ False then

16 (Step R2) Compute $\mathcal{Z}_{s}$ (Eq. (56));

17

18 else

19 (Step R2) Compute $\mathcal{Z}_{s}$ (Eq. (58));

20

21 end if

22

23 (Step R3) $\mu_{s}\leftarrow$ Centroid of
$(\widetilde{\Phi}^{g}_{m^{\prime}}(\widetilde{U}_{sm^{\prime}}))_{m^{\prime}\in\mathcal{Z}_{s}}$;

24 (Step R4) $T_{s},v_{s}\leftarrow$ Procrustes
($\mu_{s},\widetilde{\Phi}^{g}_{s}(\cup_{m^{\prime}\in\mathcal{Z}}U_{sm^{\prime}})$,
No scaling);

25 (Step R5) $\mathcal{A}\leftarrow\mathcal{A}\cup\\{s\\}$;

26

27 end for

28

29 end for

Compute $(y_{k})_{k=1}^{n}$ (Eq. (57)).

Algorithm 5 Calculate-Global-Embedding

First we look for a better than trivial initialization of the parameters so
that the views are approximately aligned. The idea is to build a rooted tree
where nodes represent the intermediate views. This tree is then traversed in a
breadth first order starting from the root. As we traverse the tree, the
intermediate view associated with a node is aligned with the intermediate view
associated with its parent node (and with a few more views), thus giving a
better initialization of the parameters. Subsequently, we refine these
parameters using a similar procedure involving random order traversal over the
intermediate views.

|
---|---
(a.1). Nine intermediate views $(\widetilde{U}_{s_{m}})_{m=1}^{9}$ of a $2d$ manifold with boundary are shown. $\widetilde{U}_{s_{9}}$ has $\widetilde{U}_{s_{7}}$ and $\widetilde{U}_{s_{8}}$ as the neighboring views. |  (a.2). In combination with (a.1), nine intermediate views $(\widetilde{U}_{s_{m}})_{m=1}^{9}$ of a closed $2d$ manifold are shown. In addition to $\widetilde{U}_{s_{7}}$ and $\widetilde{U}_{s_{8}}$, $\widetilde{U}_{s_{9}}$ also has $\widetilde{U}_{s_{1}}$ as the neighboring view.
|
(b) The intermediate views $(\widetilde{\Phi}_{s_{m}}(\widetilde{U}_{s_{m}}))_{m=1}^{9}$ in the $2$d embedding space, as they were passed as input to Algo. 5. These views are scrambled in the embedding space and Algo. 5 will move them to the right location. |  (c) The transformed views after scaling them using $b_{m}$ as in Eq. (55).
Figure 9: An illustration of the intermediate views in the ambient and the
embedding space as they are passed as input to Algo. 5 and are scaled using
Eq. (55).

##### Initialization ($\text{Iter}=1$, $\text{to\\_tear}=$ False).

In the first outer loop of Algo. 5, we start with $T_{m}=I_{d}$, $v_{m}$ as
the zero vector and compute $b_{m}$ so as to bring the intermediate views
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ to the same scale as their
counterpart $\widetilde{U}_{m}$ in the ambient space. In turn this brings all
the views to similar scale (see Figure 9 (c)). We compute the scaling
component $b_{m}$ to be the ratio of the median distance between unique points
in $\widetilde{U}_{m}$ and in $\widetilde{\Phi}_{k}(\widetilde{U}_{m})$, that
is,

$\displaystyle b_{m}=\frac{\text{median}\left\\{d_{e}(x_{k},x_{k^{\prime}})\
|\ x_{k},x_{k^{\prime}}\in\widetilde{U}_{m},x_{k}\neq
x_{k^{\prime}}\right\\}}{\text{median}\left\\{\left\|\widetilde{\Phi}_{m}(x_{k})-\widetilde{\Phi}_{m}(x_{k^{\prime}})\right\|_{2}\
|\ x_{k},x_{k^{\prime}}\in\widetilde{U}_{m},x_{k}\neq
x_{k^{\prime}}\right\\}}.$ (55)

Then we transform the the views in a sequence $(s_{m})_{m=1}^{M}$. This
sequence corresponds to the breadth first ordering of a tree starting from its
root node (which represents $s_{1}$th view). Let the $p_{s_{m}}$th view be the
parent of the $s_{m}$th view. Here $p_{s_{m}}$ lies in
$\\{s_{1},\ldots,s_{m-1}\\}$ and it is a neighboring view of the $s_{m}$th
view in the ambient space, i.e. $\widetilde{U}_{s_{m}p_{s_{m}}}$ is non-empty.
Details about the computation of these sequences is provided in Appendix D.
Note that $p_{s_{1}}$ is not defined and consequently, the first view in the
sequence ($s_{1}$th view) is not transformed, therefore $T_{s_{1}}$ and
$v_{s_{1}}$ are not updated. We also define $\mathcal{A}$, initialized with
$s_{1}$, to keep track of visited nodes which also represent the already
transformed views. Then we iterate over $m$ which varies from $2$ to $M$. For
convenience, denote the current ($m$th) node $s_{m}$ by $s$ and its parent
$p_{s_{m}}$ by $p$. The following procedure updates $T_{s}$ and $v_{s}$ (refer
to Figure 9 and 10 for an illustration of this procedure).

Step R1 ($m=9$)
---
|
(d) The transformed intermediate views $(\widetilde{\Phi}^{g}_{s_{m}}(\widetilde{U}_{s_{m}}))_{m=1}^{9}$ before the start of the iteration $m=9$. The first eight views are approximately aligned and the ninth view is to be aligned. Inaccuracies occur due to distortion. |  (e) Assuming $p_{9}=s_{7}$, step R1 computed $T_{s_{9}}$ and $v_{s_{9}}$ so that $\widetilde{\Phi}^{g}_{s_{9}}(\widetilde{U}_{s_{9}s_{7}})$ aligns with $\widetilde{\Phi}^{g}_{s_{7}}(\widetilde{U}_{s_{9}s_{7}})$. The transformed view $\widetilde{\Phi}^{g}_{s_{9}}(\widetilde{U}_{s_{9}})$ is shown. Note that step R1 results in the same output for both cases in Fig. 9 (a)
Step R2 and R3 ($m=9$)
|
(f.1) For a manifold with boundary, $\widetilde{U}_{s_{9}}$ has non-empty overlaps with $\widetilde{U}_{s_{7}}$ and $\widetilde{U}_{s_{8}}$ only. Therefore, step R2 computed $\mathcal{Z}_{s_{9}}=\\{s_{7},s_{8}\\}$. The obtained $\mu_{s_{9}}$ in step R3 is also shown in black. |  (f.2) For a closed manifold, $\widetilde{U}_{s_{9}}$ has non-empty overlaps with $\widetilde{U}_{s_{1}}$, $\widetilde{U}_{s_{7}}$ and $\widetilde{U}_{s_{8}}$. Therefore, step R2 computed $\mathcal{Z}_{s_{9}}=\\{s_{1},s_{7},s_{8}\\}$. The obtained $\mu_{s_{9}}$ in step R3 is also shown in black.
Step R4 ($m=9$)
|
(g.1) For a manifold with boundary, step R4 updated $T_{s_{9}}$ and $v_{s_{9}}$ so that the view $\widetilde{\Phi}_{s_{9}}^{g}(\widetilde{U}_{s_{9}s_{7}}\ \cup\ \widetilde{U}_{s_{9}s_{8}})$ aligns with $\mu_{s_{9}}$ in (f.1). The resulting view $\widetilde{\Phi}_{s_{9}}^{g}(\widetilde{U}_{s_{9}})$ is shown. |  (g.2) For a closed manifold step R4 updated $T_{s_{9}}$ and $v_{s_{9}}$ so that view $\widetilde{\Phi}_{s_{9}}^{g}(\widetilde{U}_{s_{9}s_{1}}\cup\widetilde{U}_{s_{9}s_{7}}\cup\widetilde{U}_{s_{9}s_{8}})$ aligns with $\mu_{s_{9}}$ in (f.2). The resulting view $\widetilde{\Phi}_{s_{9}}^{g}(\widetilde{U}_{s_{9}})$ is shown. This is not a desired output as it distorts the global embedding. We resolve this issue in Section 5.3.
Figure 10: An illustration of steps R1 to R4 in Algo. 5, in continuation of
Figure 9.

##### Step R1.

We compute a temporary value of $T_{s}$ and $v_{s}$ by aligning the views
$\widetilde{\Phi}^{g}_{s}(\widetilde{U}_{sp})$ and
$\widetilde{\Phi}^{g}_{p}(\widetilde{U}_{sp})$ of the overlap
$\widetilde{U}_{sp}$, using Procrustes analysis [21] without modifying
$b_{s}$.

##### Step R2.

Then we identify more views to align the $s$th view with. We compute a subset
$\mathcal{Z}_{s}$ of the set of already visited nodes $\mathcal{A}$ such that
$m^{\prime}\in\mathcal{Z}_{s}$ if the $s$th view and the $m^{\prime}$th view
are neighbors in the ambient space. Note that, at this stage, $\mathcal{A}$ is
the same as the set $\\{s_{1},\ldots,s_{m-1}\\}$, the indices of the first
$m-1$ views. Therefore,

$\displaystyle\mathcal{Z}_{s}=\\{m^{\prime}|\
\widetilde{U}_{sm^{\prime}}\neq\emptyset\\}\cap\mathcal{A}.$ (56)

##### Step R3.

We then compute the centroid $\mu_{s}$ of the views
$(\widetilde{\Phi}^{g}_{m^{\prime}}(\widetilde{U}_{sm^{\prime}}))_{m^{\prime}\in\mathcal{Z}_{s}}$.
Here $\mu_{s}$ is a matrix with $d$ columns and the number of rows given by
the size of the set
$\cup_{m^{\prime}\in\mathcal{Z}_{s}}\widetilde{U}_{sm^{\prime}}$. A point in
this set can have multiple embeddings due to multiple parameterizations
$(\widetilde{\Phi}^{g}_{m^{\prime}})_{m^{\prime}\in\mathcal{Z}_{s}}$ depending
on the overlaps $(\widetilde{U}_{sm^{\prime}})_{m^{\prime}\in\mathcal{Z}_{s}}$
it lies in. The mean of these embeddings forms a row in $\mu_{s}$.

##### Step R4.

Finally, we update $T_{s}$ and $v_{s}$ by aligning the view
$\widetilde{\Phi}^{g}_{s}(\widetilde{U}_{sm^{\prime}})$ with
$\widetilde{\Phi}^{g}_{m^{\prime}}(\widetilde{U}_{sm^{\prime}})$ for all
$m^{\prime}\in\mathcal{Z}_{s}$. This alignment is based on the approach in
[18, 20] where, using the Procrustes analysis [21, 31], the view
$\widetilde{\Phi}^{g}_{s}(\cup_{m^{\prime}\in\mathcal{Z}_{s}}\widetilde{U}_{sm^{\prime}})$
is aligned with the centroid $\mu_{s}$, without modifying $b_{s}$.

##### Step R5.

After the $s$th view is transformed, we add it to the set of transformed views
$\mathcal{A}$.

##### Parameter Refinement ($\text{Iter}\geq 2$, $\text{to\\_tear}=$ False).

At the end of the first iteration of the outer loop in Algo. 5, we have an
initialization of $(T_{m},b_{m},v_{m})_{m=1}^{M}$ such that transformed
intermediate views are approximately aligned. To further refine these
parameters, we iterate over $(s_{m})_{m=2}^{M}$ in random order and perform
the same five step procedure as above, $N_{r}$ times. Besides the random-order
traversal, the other difference in a refinement iteration is that the set of
already visited nodes $\mathcal{A}$, contains all the nodes instead of just
the first $m-1$ nodes. This affects the computation of $\mathcal{Z}_{s}$ (see
Eq. (56)) in step R2 so that the $s$th intermediate view is now aligned with
all those views which are its neighbors in the ambient space. Note that the
step R5 is redundant during refinement.

In the end, we compute the global embedding $y_{k}$ of $x_{k}$ by mapping
$x_{k}$ using the transformed parameterization associated with the cluster
$c_{k}$ it belongs to,

$\displaystyle y_{k}=\widetilde{\Phi}^{g}_{c_{k}}(x_{k}).$ (57)

An illustration of the global embedding at various stages of Algo. 5 is
provided in Figure 11.

|  Input | First iteration of the outer loop and stages within inner loop |  End of outer loop
---|---|---|---
|  |  Before |  Half-way |  End |
Square |  |  |  |  |
Sphere |  |  |  |  |
Figure 11: $2$d embeddings of a square and a sphere at different stages of
Algo. 5. For illustration purpose, in the plots in the $2$nd and $3$rd columns
the translation parameter $v_{m}$ was manually set for those views which do
not lie in the set $\mathcal{A}$. Note that the embedding of the sphere is
fallacious. The reason and the resolution is provided in Section 5.3.

###  5.3 Tearing Closed Manifolds

When the manifold has no boundary, then the step R2 in above section may
result in a set $\mathcal{Z}_{s}$ containing the indices of the views which
are neighbors of the $s$th view in the ambient space but are far apart from
the transformed $s$th view in the embedding space, obtained right after step
R1. For example, as shown in Figure 10 (f.2), $s_{1}\in\mathcal{Z}_{s_{9}}$
because the $s_{9}$th view and the $s_{1}$th view are neighbors in the ambient
space (see Figure 9 (a.1, a.2)) but in the embedding space, they are far
apart. Due to such indices in $\mathcal{Z}_{s_{9}}$, the step R3 results in a
centroid, which when used in step R4, results in a fallacious estimation of
the parameters $T_{s}$ and $v_{s}$, giving rise to a high distortion
embedding. By trying to align with all its neighbors in the ambient space, the
$s_{9}$th view is misaligned with respect to all of them (see Figure 10
(g.2)).

##### Resolution ($\text{to\\_tear}=$ True).

We modify the step R2 so as to introduce a discontinuity by including the
indices of only those views in the set $\mathcal{Z}_{s}$ which are neighbors
of the $s$th view in both the ambient space as well as in the embedding space.
We denote the overlap between the $m$th and $m^{\prime}$th view in the
embedding space by $\widetilde{U}^{g}_{mm^{\prime}}$. There may be multiple
heuristics for computing $\widetilde{U}^{g}_{mm^{\prime}}$ which could work.
In the Appendix E, we describe a simple approach based on the already
developed machinery in this paper, which uses the hyperparameter $\nu$
provided as input to Algo. 5. Having obtained
$\widetilde{U}^{g}_{mm^{\prime}}$, we say that the $m$th and the
$m^{\prime}$th intermediate views in the embedding space are neighbors if
$\widetilde{U}^{g}_{mm^{\prime}}$ is non-empty.

##### Step R2.

Finally, we compute $\mathcal{Z}_{s}$ as,

$\displaystyle\mathcal{Z}_{s}=\\{m^{\prime}\ |\
\widetilde{U}_{sm^{\prime}}\neq\emptyset,\;\widetilde{U}^{g}_{sm^{\prime}}\neq\emptyset\\}\cap\mathcal{A}.$
(58)

Note that if it is known apriori that the manifold can be embedded in lower
dimension without tearing it apart then we do not require the above
modification. In all of our experiments except the one in Section 6.5, we do
not assume that this information is available.

With this modification, the set $\mathcal{Z}_{s_{9}}$ in Figure 10 (f.2) will
not include $s_{1}$ and therefore the resulting centroid in the step R3 would
be the same as the one in Figure 10 (f.1). Subsequently, the transformed
$s_{9}$th view would be the one in Figure 10 (g.1) rather than Figure 10
(g.2).

##### Gluing instruction for the boundary of the embedding.

Having knowingly torn the manifold apart, we provide at the output,
information on the points belonging to the tear and their neighboring points
in the ambient space. To encode the “gluing” instructions along the tear in
the form of colors at the output of our algorithm, we recompute
$\widetilde{U}^{g}_{mm^{\prime}}$. If $\widetilde{U}_{mm^{\prime}}$ is non-
empty but $\widetilde{U}^{g}_{mm^{\prime}}$ is empty, then this means that the
$m$th and $m^{\prime}$th views are neighbors in the ambient space but are torn
apart in the embedding space. Therefore, we color the global embedding of the
points on the overlap $\widetilde{U}_{mm^{\prime}}$ which belong to clusters
$\mathcal{C}_{m}$ and $\mathcal{C}_{m^{\prime}}$ with the same color to
indicate that although these points are separated in the embedding space, they
are adjacent in the ambient space (see Figures 19, 20 and 31).

An illustration of the global embedding at various stages of Algo. 5 with
modified step R2, is provided in Figure 12.

|  Input | First iteration of the outer loop and stages within inner loop |  End of outer loop
---|---|---|---
|  |  Before |  Half-way |  End |
Sphere |  |  |  |  |
Figure 12: $2$d embedding of a sphere at different stages of Algo. 5. For
illustration purpose, in the plots in the $2$nd and $3$rd columns the
translation parameter $v_{m}$ was manually set for those views which do not
lie in the set $\mathcal{A}$.

##### Example.

The obtained global embeddings of our square grid with
$\text{to\\_tear}=\text{True}$ and $\nu=3$, are shown in Figure 13. Note that
the boundary of the obtained embedding is more distorted when the points on
the boundary are unknown than when they are known apriori. This is because the
intermediate views near the boundary have higher distortion in the former case
than in the latter case (see Figure 7).

|
---|---
Figure 13: Global embedding of the square grid when the points on the boundary
are unknown (left) versus when they are known apriori (right).

###  5.4 Time Complexity

The worst case time complexity of Algo. 5 is
$O(N_{r}nk_{\text{lv}}^{2}d^{2}/\eta_{\text{min}})$ when to_tear is false. It
costs an additional time of $O(N_{r}n^{2}\text{max}(d,k_{\text{lv}}\log
n,n/\eta_{\text{min}}^{2})))$ when to_tear is true. In practice, one
refinement step took about $15$ seconds in the above example and between
$15$-$20$ seconds for all the examples in Section 6.

##  6 Experimental Results

We present experiments to compare LDLE222The python code is available at
https://github.com/chiggum/pyLDLE with LTSA [49], UMAP [32], t-SNE [30] and
Laplacian eigenmaps [3] on several datasets. First, we compare the embeddings
of discretized $2$d manifolds embedded in $\mathbb{R}^{2}$, $\mathbb{R}^{3}$
or $\mathbb{R}^{4}$, containing about $10^{4}$ points. These manifolds are
grouped based on the presence of the boundary and their orientability as in
Sections 6.2, 6.3 and 6.4. The inputs are shown in the figures themselves
except for the flat torus and the Klein bottle, as their $4$D
parameterizations cannot be plotted. Therefore, we describe their construction
below. A quantitative comparison of the algorithms is provided in Section
6.2.1. In Section 6.2.2 we assess the robustness of these algorithms to the
noise in the data. In Section 6.2.3 we assess the performance of these
algorithms on sparse data. Finally, in Section 6.5 we compare the embeddings
of some high dimensional datasets.

Flat Torus. A flat torus is a parallelogram whose opposite sides are
identified. In our case, we construct a discrete flat torus using a rectangle
with sides $2$ and $0.5$ and embed it in four dimensions as follows,

$\displaystyle X(\theta_{i},\phi_{j})$
$\displaystyle=\frac{1}{4\pi}(4cos(\theta_{i}),4\sin(\theta_{i}),\cos(\phi_{j}),\sin(\phi_{j}))$
(59)

where $\theta_{i}=0.01i\pi$, $\phi_{j}=0.04j\pi$, $i\in\\{0,\ldots,199\\}$ and
$j\in\\{0,\ldots,49\\}$.

Klein bottle. A Klein bottle is a non-orientable two dimensional manifold
without boundary. We construct a discrete Klein bottle using its $4$D Möbius
tube representation as follows,

$\displaystyle X(\theta_{i},\phi_{j})$
$\displaystyle=(R(\phi_{j})\cos\theta_{i},R(\phi_{j})\sin\theta_{i},r\sin\phi_{j}\cos\frac{\theta_{i}}{2},r\sin\phi_{j}\sin\frac{\theta_{i}}{2})$
(60) $\displaystyle R(\phi_{j})$ $\displaystyle=R+r\cos\phi_{j}$ (61)

where $\theta_{i}=i\pi/100$, $\phi_{j}=j\pi/25$, $i\in\\{0,\ldots,199\\}$ and
$j\in\\{0,\ldots,49\\}$.

###  6.1 Hyperparameters

To embed using LDLE, we use the Euclidean metric and the default values of the
hyperparameters and their description are provided in Table 1. Only the value
of $\eta_{\text{min}}$ is tuned across all the examples in Sections 6.2, 6.3
and 6.4 (except for Section 6.2.3), and is provided in Appendix G. For high
dimensional datasets in Section 6.5, values of the hyperaparameters which
differ from the default values are again provided in Appendix G.

Hyper- parameter |  Description |  Default value
---|---|---
$k_{\text{nn}}$ |  No. of nearest neighbors used to construct the graph Laplacian | $49$
$k_{\text{tune}}$ |  The nearest neighbor, distance to which is used as a local scaling factor in the construction of graph Laplacian | 7
$N$ |  No. of nontrivial low frequency Laplacian eigenvectors to consider for the construction of local views in the embedding space | 100
$d$ |  Intrinsic dimension of the underlying manifold | 2
$p$ |  Probability mass for computing the bandwidth $t_{k}$ of the heat kernel | 0.99
$k_{\text{lv}}$ |  The nearest neighbor, distance to which is used to construct local views in the ambient space | 25
$(\tau_{s})_{s=1}^{d}$ |  Percentiles used to restrict the choice of candidate eigenfunctions | 50
$(\delta_{s})_{s=1}^{d}$ |  Fractions used to restrict the choice of candidate eigenfunctions | 0.9
$\eta_{\text{min}}$ |  Desired minimum number of points in a cluster | 5
to_tear |  A boolean for whether to tear the manifold or not | True
$\nu$ |  A relaxation factor to compute the neighborhood graph of the intermediate views in the embedding space | 3
$N_{r}$ |  No. of iterations to refine the global embedding | 100
Table 1: Default values of LDLE hyperparameters.

For UMAP, LTSA, t-SNE and Laplacian eigenmaps, we use the Euclidean metric and
select the hyperparameters by grid search, choosing the values which result in
best visualization quality. For LTSA, we search for optimal n_neighbors in
$\\{5,10,25,50,75,100\\}$. For UMAP, we use $500$ epochs and search for
optimal n_neighbors in $\\{25,50,100,200\\}$ and min_dist in
$\\{0.01,0.1,0.25,0.5\\}$. For t-SNE, we use $1000$ iterations and search for
optimal perplexity in $\\{30,40,50,60\\}$ and early exaggeration in
$\\{2,4,6\\}$. For Laplacian eigenmaps, we search for $k_{\text{nn}}$ in
$\\{16,25,36,49\\}$ and $k_{\text{tune}}$ in $\\{3,7,11\\}$. The chosen values
of the hyperparameters are provided in Appendix G. We note that the Laplacian
eigenmaps fails to correctly embed most of the examples regardless of the
choice of the hyperparameters.

###  6.2 Manifolds with Boundary

In Figure 14, we show the $2$d embeddings of $2$d manifolds with boundary, in
$\mathbb{R}^{2}$ or $\mathbb{R}^{3}$, three of which have holes. To a large
extent, LDLE preserved the shape of the holes. LTSA perfectly preserved the
shape of the holes in the square but deforms it in the Swiss Roll. This is
because LTSA embedding does not capture the aspect ratio of the underlying
manifold as discussed in Section F. UMAP and Laplacian eigenmaps distorted the
shape of the holes and the region around them, while t-SNE produced dissected
embeddings. For the sphere with a hole which is a curved $2$d manifold with
boundary, LTSA, UMAP and Laplacian eigenmaps squeezed it into $\mathbb{R}^{2}$
while LDLE and t-SNE tore it apart. The correctness of the LDLE embedding is
proved in Figure 31. In the case of noisy swiss roll, LDLE and UMAP produced
visually better embeddings in comparison to the other methods.

We note that the boundaries of the LDLE embeddings in Figure 14 are usually
distorted. The cause of this is explained in Remark 4. When the points in the
input which lie on the boundary are known apriori then the distortion near the
boundary can be reduced using the double manifold as discussed in Remark 5 and
shown in Figure 4. The obtained LDLE embeddings when the points on the
boundary are known, are shown in Figure 15.

|  Barbell |  Square with two holes |  Sphere with a hole |  Swiss Roll with a hole |  Noisy Swiss Roll
---|---|---|---|---|---
Input |  |  |  |  |
LDLE |  |  |  |  |
LTSA |  |  |  |  |
UMAP |  |  |  |  |
t-SNE |  |  |  |  |
Laplacian Eigenmaps |  |  |  |  |
Figure 14: Embeddings of $2$d manifolds with boundary into $\mathbb{R}^{2}$. The noisy Swiss Roll is constructed by adding uniform noise in all three dimensions, with support on $[0,0.05]$. |  Barbell |  Square with two holes |  Swiss Roll with a hole
---|---|---|---
LDLE with $\partial\mathcal{M}$ known apriori |  |  |
Figure 15: LDLE embeddings when the points on the boundary are known apriori.

####  6.2.1 Quantitative comparison

To compare LDLE with other techniques in a quantitative manner, we compute the
distortion $\mathcal{D}_{k}$ of the embeddings of the geodesics originating
from $x_{k}$ and then plot the distribution of $\mathcal{D}_{k}$ (see Figure
16). The procedure to compute $\mathcal{D}_{k}$ follows. In the discrete
setting, we first define the geodesic between two given points as the shortest
path between them which in turn is computed by running Dijkstra algorithm on
the graph of $5$ nearest neighbors. Here, the distances are measured using the
Euclidean metric $d_{e}$. Denote the number of nodes on the geodesic between
$x_{k}$ and $x_{k^{\prime}}$ by $n_{kk^{\prime}}$ and the sequence of nodes by
$(x_{s})_{s=1}^{n_{kk^{\prime}}}$ where $x_{1}=x_{k}$ and
$x_{n_{kk^{\prime}}}=x_{k^{\prime}}$. Denote the embedding of $x_{k}$ by
$y_{k}$. Then the length of the geodesic in the latent space between $x_{k}$
and $x_{k^{\prime}}$, and the length of the embedding of the geodesic between
$y_{k}$ and $y_{k^{\prime}}$ are given by

$\displaystyle L_{kk^{\prime}}$
$\displaystyle=\sum_{s=2}^{n_{kk^{\prime}}}d_{e}(x_{s},x_{s-1}).$ (62)
$\displaystyle L^{g}_{kk^{\prime}}$
$\displaystyle=\sum_{s=2}^{n_{kk^{\prime}}}d_{e}(y_{s},y_{s-1}).$ (63)

Finally, the distortion $\mathcal{D}_{k}$ of the embeddings of the geodesics
originating from $x_{k}$ is given by the ratio of maximum expansion and
minimum contraction, that is,

$\displaystyle\mathcal{D}_{k}$
$\displaystyle=\sup_{k^{\prime}}\frac{L^{g}_{kk^{\prime}}}{L_{kk^{\prime}}}/\inf_{k^{\prime}}\frac{L^{g}_{kk^{\prime}}}{L_{kk^{\prime}}}=\sup_{k^{\prime}}\frac{L^{g}_{kk^{\prime}}}{L_{kk^{\prime}}}\sup_{k^{\prime}}\frac{L_{kk^{\prime}}}{L^{g}_{kk^{\prime}}}.$
(64)

A value of $1$ for $\mathcal{D}_{k}$ means the geodesics originating from
$x_{k}$ have the same length in the input and in the embedding space. If
$\mathcal{D}_{k}=1$ for all $k$ then the embedding is geometrically, and
therefore topologically as well, the same as the input up to scale. Figure 16
shows the distribution of $\mathcal{D}_{k}$ due to LDLE and other algorithms
for various examples. Except for the noisy Swiss Roll, LTSA produced the least
maximum distortion. Specifically, for the square with two holes, LTSA produced
a distortion of $1$ suggesting its strength on manifolds with unit aspect
ratio. In all other examples, LDLE produced the least distortion except for a
few outliers. When the boundary is unknown, the points which result in high
$\mathcal{D}_{k}$ are the ones which lie on and near the boundary. When the
boundary is known, these are the points which lie on or near the corners (see
Figures 4 and 5). We aim to fix this issue in future work.

####  6.2.2 Robustness to noise

To further analyze the robustness of LDLE under noise we compare the
embeddings of the Swiss Roll with Gaussian noise of increasing variance. The
resulting embeddings are shown in Figure 17. Note that certain points on LDLE
embeddings have a different colormap than the one used for the input. As
explained in Section 5.3, the points which have the same color under this
colormap are adjacent on the manifold but away in the embedding. To be
precise, these points lie close to the middle of the gap in the Swiss Roll,
creating a bridge between those points which would otherwise be far away on a
noiseless Swiss Roll. In a sense, these points cause maximum corruption to the
geometry of the underlying noiseless manifold. One can say that these points
are have adversarial noise, and LDLE embedding can automatically recognize
such points. We will further explore this in future work. LTSA, t-SNE and
Laplacian Eigenmaps fail to produce correct embeddings while UMAP embeddings
also exhibit high quality.

####  6.2.3 Sparsity

A comparison of the embeddings of the Swiss Roll with decreasing resolution
and increasing sparsity is provided in Figure 18. Unlike LTSA and Laplacian
Eigenmaps, the embeddings produced by LDLE, UMAP and t-SNE are of high
quality. Note that when the resolution is $10$, LDLE embedding of some points
have a different colormap. Due to sparsity, certain points on the opposite
sides of the gap between the Swiss Roll are neighbors in the ambient space as
shown in Figure 32 in Appendix I. LDLE automatically tore apart these
erroneous connections and marked them at the output using a different
colormap. A discussion on sample size requirement for LDLE follows.

The distortion of LDLE embeddings directly depend on the distortion of the
constructed local parameterizations, which in turn depends on reliable
estimates of the graph Laplacian and its eigenvectors. The work in [4, 22, 45,
13] provided conditions on the sample size and the hyperparameters such as the
kernel bandwidth, under which the graph Laplacian and its eigenvectors would
converge to their continuous counterparts. A similar analysis in the setting
of self-tuned kernels used in our approach (see Algo. 1) is also provided in
[12]. These imply that, for a faithful estimation of graph Laplacian and its
eigenvectors, the hyperparameter $k_{\text{tune}}$ (see Table 1) should be
small enough so that the local scaling factors $\sigma_{k}$ (see Algo. 1) are
also small, while the size of the data $n$ should be large enough so that
$n\sigma_{k}^{d+2}/\log(n)$ is sufficiently large for all
$k\in\\{1,\ldots,n\\}$. This suggests that $n$ needs to be exponential in $d$
and inversely related to $\sigma_{k}$. However, in practice, the data is
usually given and therefore $n$ is fixed. So the above mainly states that to
obtain accurate estimates, the hyperparameter $k_{\text{tune}}$ must be
decreased. This indeed holds as we had to decrease $k_{\text{tune}}$ from $7$
to $2$ (see Appendix G) to produce LDLE embeddings of high quality for
increasingly sparse Swiss Roll in Figure 18.

###  6.3 Closed Manifolds

In Figure 19, we show the $2$d embeddings of $2$d manifolds without a
boundary, a curved torus in $\mathbb{R}^{3}$ and a flat torus in
$\mathbb{R}^{4}$. LDLE produced similar representation for both the inputs.
None of the other methods do that. The main difference in the LDLE embedding
of the two inputs is based on the boundary of the embedding. It is composed of
many small line segments for the flat torus, and many small curved segments
for the curved torus. This is clearly because of the difference in the
curvature of the two inputs, zero everywhere for the flat torus and non-zero
almost everywhere on the curved torus. The mathematical correctness of the
LDLE embeddings using the cut and paste argument is shown in Figure 31. LTSA,
UMAP and Laplacian eignemaps squeezed both the manifolds into $\mathbb{R}^{2}$
while the t-SNE embedding is non-interpretable.

###  6.4 Non-Orientable Manifolds

In Figure 20, we show the $2$d embeddings of non-orientatble $2$d manifolds, a
Möbius strip in $\mathbb{R}^{3}$ and a Klein bottle in $\mathbb{R}^{4}$.
Laplacian eigenmaps produced incorrect embeddings, t-SNE produced dissected
and non-interpretable embeddings and LTSA and UMAP squeezed the inputs into
$\mathbb{R}^{2}$. LDLE produced mathematically correct embeddings by tearing
apart both inputs to embed them into $\mathbb{R}^{2}$ (see Figure 31).

###  6.5 High Dimensional Data

####  6.5.1 Synthetic sensor data

In Figure 21, motivated from [36], we embed a $42$ dimensional synthetic data
set representing the signal strength of $42$ transmitters at about $n=6000$
receiving locations on a toy floor plan. The transmitters and the receivers
are distributed uniformly across the floor. Let $(t_{r_{k}})_{k=1}^{42}$ be
the transmitter locations and $r_{i}$ be the $i$th receiver location. Then the
$i$th data point $x_{i}$ is given by
$(e^{-\left\|r_{i}-t_{r_{k}}\right\|_{2}^{2}})_{k=1}^{42}$. The resulting data
set is embedded using and other algorithms into $\mathbb{R}^{2}$. The
hyperparameters resulting in the most visually appealing embeddings were
identified for each algorithm and are provided in Table 2. The obtained
embeddings are shown in Figure 21. The shapes of the holes are best preserved
by LTSA, then LDLE followed by the other algorithms. The corners of the LDLE
embedding are more distorted. The reason for distorted corners is given in
Remark 4.

####  6.5.2 Face image data

In Figure 22, we show the embedding obtained by applying LDLE on the face
image data [44] which consists of a sequence of $698$ $64$-by-$64$ pixel
images of a face rendered under various pose and lighting conditions. These
images are converted to $4096$ dimensional vectors, then projected to $100$
dimensions through PCA while retaining about $98\%$ of the variance. These are
then embedded using LDLE and other algorithms into $\mathbb{R}^{2}$. The
hyperparameters resulting in the most visually appealing embeddings were
identified for each algorithm and are provided in Table 5. The resulting
embeddings are shown in Figure 23 colored by the pose and lighting of the
face. Note that values of the pose and lighting variables for all the images
are provided in the dataset itself. We have displayed face images
corresponding to few points of the LDLE embeddings as well. Embeddings due to
all the techniques except LTSA reasonably capture both the pose and lighting
conditions.

####  6.5.3 Rotating Yoda-Bulldog dataset

In Figure 23, we show the $2$d embeddings of the rotating figures dataset
presented in [28]. It consists of $8100$ snapshots taken by a camera of a
platform with two objects, Yoda and a bull dog, rotating at different
frequencies. Therefore, the underlying $2$d parameterization of the data
should render a torus. The original images have a dimension of $320\times
240\times 3$. In our experiment, we first resize the images to half the
original size and then project them to $100$ dimensions through PCA [24] while
retaining about $98\%$ variance. These are then embedded using LDLE and other
algorithms into $\mathbb{R}^{2}$. The hyperparameters resulting in the most
visually appealing embeddings were identified for each algorithm and are
provided in Table 5. The resulting embeddings are shown in Figure 23 colored
by the first dimension of the embedding itself. LTSA and UMAP resulted in a
squeezed torus. LDLE tore apart the underlying torus and automatically colored
the boundary of the embedding to suggest the gluing instructions. By tracing
the color on the boundary we have manually drawn the arrows. Putting these
arrows on a piece of paper and using cut and past argument one can establish
that the embedding represents a torus (see Figure 31). The images
corresponding to a few points on the boundary are shown. Pairs of images with
the same labels represent the two sides of the curve along which LDLE tore
apart the torus, and as is evident these pairs are similar.

##  7 Conclusion and Future Work

We have presented a new bottom-up approach (LDLE) for manifold learning which
constructs low-dimensional low distortion local views of the data using the
low frequency global eigenvectors of the graph Laplacian, and registers them
to obtain a global embedding. Through various examples we demonstrated that
LDLE competes with the other methods in terms of visualization quality. In
particular, the embeddings produced by LDLE preserved distances upto a
constant scale better than those produced by UMAP, t-SNE, Laplacian Eigenmaps
and for the most part LTSA too. We also demonstrated that LDLE is robust to
the noise in the data and produces fine embeddings even when the data is
sparse. We also showed that LDLE can embed closed as well as non-orientable
manifolds into their intrinsic dimension, a feature that is missing from the
existing techniques. Some of the future directions of our work are as follows.

  * •

It is only natural to expect real world datasets to have boundary and to have
many corners. As observed in the experimental results, when the boundary of
the manifold is unknown, then the LDLE embedding tends to have distorted
boundary. Even when the boundary is known, the embedding has distorted
corners. This is caused by high distortion views near the boundary (see
Figures 4 and 5). We aim to fix this issue in our future work. One possible
resolution could be based on [5] which presented a method to approximately
calculate the distance of the points from the boundary.

  * •

When the data represents a mixture of manifolds, for example, a pair of
possibly intersecting spheres or even manifolds of different intrinsic
dimensions, it is also natural to expect a manifold learning technique to
recover a separate parameterization for each manifold and provide gluing
instructions at the output. One way is to perform manifold factorization [48]
or multi-manifold clustering [46] on the data to recover sets of points
representing individual manifolds and then use manifold learning on these
separately. We aim to adapt LDLE to achieve this.

  * •

The spectrum of the Laplacian has been used in prior work for anomaly
detection [15, 33, 11, 10, 35]. Similar to our approach of using a subset of
Laplacian eigenvectors to construct low distortion local views in lower
dimension, in [34, 10], subsets of Laplacian eigenvectors were identified so
as to separate small clusters from a large background component. As shown in
Figures 4 and 5, LDLE produced high distortion local views near the boundary
and the corners, though these are not outliers. However, if we consider a
sphere with outliers (imagine a sphere with noise only at the north pole as in
Figure 24), then the distortion of the local views containing the outliers is
higher than the rest of the views. Therefore, the distortion of the local
views can help find anomalies in the data. We aim to further investigate this
direction to develop an anomaly detection technique.

  * •

Similar to the approach of denoising a signal by retaining low frequency
components, our approach uses low frequency Laplacian eigenvectors to estimate
local views. These eigenvectors implicitly capture the global structure of the
manifold. Therefore, to construct local views, unlike LTSA which directly
relies on the local configuration of data which may be noisy, LDLE relies on
the local elements of low frequency global eigenvectors of the Laplacian which
are supposed to be robust to the noise. Practical implication of this is shown
in Figure 17 to some extent while we aim to further investigate the
theoretical implications.

Rectangle ($4\times 0.25$) | barbell |  Square with two holes |  Swiss Roll with a hole | Noisy Swiss Roll
---|---|---|---|---
|  |  |  |
Figure 16: Violin plots [23, 2] for the distribution of $\mathcal{D}_{k}$ (See Eq. (64)). LDLE $\partial M$ means LDLE with boundary known apriori. The white point inside the violin represents the median. The straight line above the end of the violin represents the outliers. |  $\sigma=0.01$ |  $\sigma=0.015$ |  $\sigma=0.02$
---|---|---|---
Side view of Swiss Roll |  |  |
LDLE |  |  |
LTSA |  |  |
UMAP |  |  |
t-SNE |  |  |
Laplacian Eigenmaps |  |  |
Figure 17: Embeddings of the Swiss Roll with additive noise sampled from the Gaussian distribution of zero mean and a variance of $\sigma^{2}$ (see Section 6.2.2 for details). |  RES$=30$ ($n=990$) |  RES$=15$ ($n=280$) |  RES$=12$ ($n=184$) |  RES$=10$ ($n=133$)
---|---|---|---|---
Input |  |  |  |
LDLE |  |  |  |
LTSA |  |  |  |
UMAP |  |  |  |
t-SNE |  |  |  |
Laplacian Eigenmaps |  |  |  |
Figure 18: Embeddings of the Swiss Roll with decreasing resolution and increasing sparsity (see Section 6.2.3 for details). Note that when RES$=7$ ($n=70$) none of the above method produced a correct embedding. | Curved torus | Flat torus
---|---|---
Input |  |  |  See Eq. (59)
LDLE |  |  |  |
LTSA |  |  |  |
UMAP |  |  |  |
t-SNE |  |  |  |
Laplacian Eigenmaps |  |  |  |
Figure 19: Embeddings of $2$d manifolds without boundary into $\mathbb{R}^{2}$. For each manifold, the left and right columns contain the same plots colored by the two parameters of the manifold. A proof of the mathematical correctness of the LDLE embeddings is provided in Figure 31. | Möbius strip | Klein bottle
---|---|---
Input |  |  | See Eq. (61)
LDLE |  |  |  |
LTSA |  |  |  |
UMAP |  |  |  |
t-SNE |  |  |  |
Laplacian Eigenmaps |  |  |  |
Figure 20: Embeddings of $2$d non-orientable manifolds into $\mathbb{R}^{2}$.
For each manifold, the left and right columns contain the same plots colored
by the two parameters of the manifold. A proof of the mathematical correctness
of the LDLE embeddings is provided in Figure 31.

True floor plan |  LDLE |  LTSA |  UMAP |  t-SNE |  Laplacian eigenmaps
---|---|---|---|---|---
|  |  |  |  |
Figure 21: Embedding of the synthetic sensor data into $\mathbb{R}^{2}$ (see Section 6.5 for details). | LDLE
---|---
|
|  LDLE |  LTSA |  UMAP |  t-SNE
Pose |  |  |  |
Lighting |  |  |  |
Figure 22: Embedding of the face image data set [44] into $\mathbb{R}^{2}$
colored by the pose and lighting conditions (see Section 6.5 for details).
LDLE
---
LTSA |  UMAP |  t-SNE
|  |
Figure 23: Embeddings of snapshots of a platform with two objects, Yoda and a bull dog, each rotating at a different frequency, such that the underlying topology is a torus (see Section 6.5 for details). |  |
---|---|---
Figure 24: Local views containing outliers exhibit high distortion. (left)
Input data $(x_{k})_{k=1}^{n}$. (middle) $x_{k}$ colored by the distortion
$\zeta_{kk}$ of $\Phi_{k}$ on $U_{k}$. (right) $y_{k}$ colored by
$\zeta_{kk}$.

##  Appendix A First Proof of Theorem 2

Choose $\epsilon>0$ so that the exponential map
$\exp_{x}:T_{x}\mathcal{M}\rightarrow\mathcal{M}$ is a well defined
diffeomorphism on $\mathcal{B}_{2\epsilon}\subset T_{x}\mathcal{M}$ where
$T_{x}\mathcal{M}$ is the tangent space to $\mathcal{M}$ at $x$,
$\exp_{x}(0)=x$ and

$\displaystyle\mathcal{B}_{\epsilon}=\\{v\in T_{x}\mathcal{M}\ |\
\left\|v\right\|_{2}<\epsilon\\}.$ (65)

Then using [7, lem. 48, prop. 50, th. 51], for all $y\in B_{\epsilon}(x)$ such
that

$\displaystyle B_{\epsilon}(x)=\\{y\in\mathcal{M}\ |\ d_{g}(x,y)<\epsilon\\}$
(66)

we have,

$\displaystyle p(t,x,y)$
$\displaystyle=G(t,x,y)(u_{0}(x,y)+tu_{1}(x,y)+O(t^{2})),$ (67)

where

$\displaystyle G(t,x,y)$ $\displaystyle=\frac{e^{-d_{g}(x,y)^{2}/4t}}{(4\pi
t)^{d/2}},$ (68) $\displaystyle u_{0}(x,y)$
$\displaystyle=1+O(\left\|v\right\|^{2}),\ y=\exp_{x}(v),v\in
T_{x}\mathcal{M},$ (69)

and for $f\in C(\mathcal{M})$, the following hold

$\displaystyle f(x)$ $\displaystyle=\lim_{t\rightarrow
0}\int_{M}p(t,x,y)f(y)\omega_{g}(y)$ (70) $\displaystyle=\lim_{t\rightarrow
0}\int_{B_{\epsilon}(x)}p(t,x,y)f(y)\omega_{g}(y),$ (71) $\displaystyle f(x)$
$\displaystyle=\lim_{t\rightarrow
0}\int_{B_{\epsilon}(x)}G(t,x,y)f(y)\omega_{g}(y),$ (72) $\displaystyle
u_{1}(x,x)f(x)$ $\displaystyle=\lim_{t\rightarrow
0}\int_{B_{\epsilon}(x)}G(t,x,y)u_{1}(x,y)f(y)\omega_{g}(y).$ (73)

Using the above equations and the definition of $\Psi_{kij}(y)$ in Eq. (15)
and $A_{kij}$ in Eq. (16) we compute the limiting value of the scaled local
correlation (see Eq. (19)),

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t\rightarrow
0}\frac{A_{kij}}{2t}$ (74) $\displaystyle=\lim_{t\rightarrow
0}\frac{1}{2t}\int_{M}p(t,x_{k},y)\Psi_{kij}(y)\omega_{g}(y).$ (75)

which will turn out to be the inner product between the gradients of the
eigenfunctions $\bm{\phi}_{i}$ and $\bm{\phi}_{j}$ at $x_{k}$. We start by
choosing an $\epsilon_{k}>0$ so that $\exp_{x_{k}}$ is a well defined
diffeomorphism on $\mathcal{B}_{2\epsilon_{k}}\subset T_{x_{k}}\mathcal{M}$.
Using Eq. (71) we change the region of integration from $\mathcal{M}$ to
$B_{\epsilon_{k}}(x_{k})$,

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{B_{\epsilon_{k}}(x_{k})}p(t_{k},x_{k},y)\Psi_{kij}(y)\omega_{g}(y).$
(76)

Substitute $p(t_{k},x_{k},y)$ from Eq. (67) and simplify using Eq. (72, 73)
and the fact that $\Psi_{kij}(x_{k})=0$ to get

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{B_{\epsilon_{k}}(x_{k})}G(t_{k},x_{k},y)(u_{0}(x_{k},y)+t_{k}u_{1}(x_{k},y)+O(t_{k}^{2}))\Psi_{kij}(y)\omega_{g}(y).$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\left(\frac{1}{2t_{k}}\int_{B_{\epsilon_{k}}(x_{k})}G(t_{k},x_{k},y)u_{0}(x_{k},y)\Psi_{kij}(y)\omega_{g}(y)+\right.$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\left.\frac{t_{k}u_{1}(x_{k},x_{k})\Psi_{kij}(x_{k})+O(t_{k}^{2})\Psi_{kij}(x_{k})}{2t_{k}}\right)$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{B_{\epsilon_{k}}(x_{k})}G(t_{k},x_{k},y)u_{0}(x_{k},y)\Psi_{kij}(y)\omega_{g}(y).$
(77)

Replace $y\in B_{\epsilon_{k}}(x_{k})$ by $\exp_{x_{k}}(v)$ where
$v\in\mathcal{B}_{\epsilon_{k}}\subset T_{x_{k}}\mathcal{M}$ and
$\left\|v\right\|=d_{g}(x_{k},y)$. Denote the Jacobian for the change of
variable by $J(v)$ i.e. $J(v)=\frac{d}{dv}\exp_{x_{k}}(v)$. Note that
$\exp_{x_{k}}(0)=x_{k}$ and $J(0)=I$. Using the Taylor expansion of
$\bm{\phi}_{i}$ and $\bm{\phi}_{j}$ about $0$ we obtain

$\displaystyle\phi_{s}(y)=\phi_{s}(\exp_{x_{k}}(v))$
$\displaystyle=\phi_{s}(\exp_{x_{k}}(0))+\nabla\phi_{s}(\exp_{x_{k}}(0))^{T}J(0)v+O(\left\|v\right\|^{2})$
$\displaystyle=\phi_{s}(x_{k})+\nabla\phi_{s}(x_{k})^{T}v+O(\left\|v\right\|^{2}),\
s=i,j.$ (78)

Substituting the above equation in the definition of $\Psi_{kij}(y)$ (see Eq.
(15)) we get

$\displaystyle\Psi_{kij}(y)$ $\displaystyle=\Psi_{kij}(\exp_{x_{k}}(v))$
$\displaystyle=v^{T}\nabla\phi_{i}\nabla\phi_{j}^{T}v+(\nabla\phi_{i}^{T}v+\nabla\phi_{j}^{T}v)O(\left\|v\right\|^{2})+O(\left\|v\right\|^{4}),$
(79)

where $\nabla\phi_{s}\equiv\nabla\phi_{s}(x_{k}),s=i,j$. Now we substitute Eq.
(79, 68, 69) in Eq. (77) while replacing variable $y$ with $\exp_{x_{k}}(v)$
where $J(v)$ is the Jacobian for the change of variable as before, to get

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{\mathcal{B}_{\epsilon_{k}}}\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}(1+O(\left\|v\right\|^{2}))\Psi_{kij}(\exp_{x_{k}}(v))J(v)dv$
$\displaystyle=L_{1}+L_{2},$ (80)

where $L_{1}$ and $L_{2}$ are the terms obtained by expanding
$1+O(\left\|v\right\|^{2})$ in the integrand. We will show that $L_{2}=0$ and
$\widetilde{A}_{kij}=L_{1}=\nabla\phi_{i}^{T}\nabla\phi_{j}$.

$\displaystyle L_{2}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{\mathcal{B}_{\epsilon_{k}}}\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}O(\left\|v\right\|^{2})(\operatorname{tr}(\nabla\phi_{i}\nabla\phi_{j}^{T}vv^{T})+$
$\displaystyle\qquad\qquad\qquad\qquad\qquad(\nabla\phi_{i}^{T}v+\nabla\phi_{j}^{T}v)O(\left\|v\right\|^{2})+O(\left\|v\right\|^{4}))J(v)dv$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}(O(t_{k}^{2})+0+0+O(t_{k}^{4}))$ $\displaystyle=0.$ (81)

Therefore,

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=L_{1}$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{\mathcal{B}_{\epsilon_{k}}}\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}\Psi_{kij}(\exp_{x_{k}}(v))J(v)dv$ (82)
$\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{\mathcal{B}_{\epsilon_{k}}}\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}(v^{T}\nabla\phi_{i}\nabla\phi_{j}^{T}v+$
$\displaystyle\qquad\qquad\qquad\qquad(\nabla\phi_{i}^{T}v+\nabla\phi_{j}^{T}v)O(\left\|v\right\|^{2})+O(\left\|v\right\|^{4}))J(v)dv$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{\mathcal{B}_{\epsilon_{k}}}\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}v^{T}\nabla\phi_{i}\nabla\phi_{j}^{T}vJ(v)dv+\frac{0+0+O(t_{k}^{2})}{2t_{k}}$
$\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{\mathcal{B}_{\epsilon_{k}}}\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}v^{T}\nabla\phi_{i}\nabla\phi_{j}^{T}vJ(v)dv.$ (83)

Substitution of $t_{k}=0$ leads to the indeterminate form $\frac{0}{0}$.
Therefore, we apply L’Hospital’s rule and then Leibniz integral rule to get,

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2}\int_{\mathcal{B}_{\epsilon_{k}}}\left(\frac{\left\|v\right\|^{2}}{4t_{k}^{2}}-\frac{d}{2t_{k}}\right)\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}v^{T}\nabla\phi_{i}\nabla\phi_{j}^{T}vJ(v)dv$
$\displaystyle=\operatorname{tr}\left(\frac{1}{2}\nabla\phi_{i}\nabla\phi_{j}^{T}\lim_{t_{k}\rightarrow
0}\int_{\mathcal{B}_{\epsilon_{k}}}\left(\frac{\left\|v\right\|^{2}}{4t_{k}^{2}}-\frac{d}{2t_{k}}\right)\frac{e^{-\left\|v\right\|^{2}/4t_{k}}}{(4\pi
t_{k})^{d/2}}vv^{T}J(v)dv\right)$
$\displaystyle=\operatorname{tr}\left(\frac{1}{2}\nabla\phi_{i}\nabla\phi_{j}^{T}\left(\lim_{t_{k}\rightarrow
0}\left(\frac{(12+4(d-1))t_{k}^{2}}{4t_{k}^{2}}-\frac{2t_{k}d}{2t_{k}}\right)I+O(t_{k})I\right)\right)$
$\displaystyle=\nabla\phi_{i}^{T}\nabla\phi_{j}.$ (84)

Finally, note that the Eq. (82) is same as the following equation with $y$
replaced by $\exp_{x_{k}}(v)$,

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{1}{2t_{k}}\int_{B_{\epsilon_{k}}(x_{k})}G(t_{k},x_{k},y)\Psi_{kij}(y)\omega_{g}(y).$
(85)

We used the above equation to estimate $\widetilde{A}_{kij}$ in Section 3.1. ∎

##  Appendix B Second Proof of Theorem 2

Yet another proof is based on the Feynman-Kac formula [41, 42],

$\displaystyle A_{kij}$
$\displaystyle=[e^{-t_{k}\Delta_{g}}((\phi_{i}-\phi_{i}(x_{k}))(\phi_{j}-\phi_{j}(x_{k}))](x_{k}).$
(86)

where

$\displaystyle[e^{-t\Delta_{g}}f](x)=\sum_{i}e^{-\lambda_{i}t}\langle\phi_{i},f\rangle\phi_{i}(x)$
(87)

and therefore,

$\displaystyle\widetilde{A}_{kij}$ $\displaystyle=\lim_{t_{k}\rightarrow
0}\frac{A_{kij}}{2t_{k}}=\left.\frac{1}{2}\frac{\partial A_{kij}}{\partial
t_{k}}\right|_{t_{k}=0}$ (88)
$\displaystyle=\frac{-1}{2}\left\\{\Delta_{g}[(\phi_{i}-\phi_{i}(x_{k}))(\phi_{j}-\phi_{j}(x_{k}))](x_{k})\right\\}$
(89)
$\displaystyle=\frac{-1}{2}\left\\{0+0-2\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k})\right\\}$
(90) $\displaystyle=\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k})$ (91)

where we used the fact
$\Delta_{g}(f_{i}f_{j})=f_{j}\Delta_{g}f_{i}+f_{i}\Delta_{g}f_{j}-2\langle\nabla_{g}f_{i}(x),\nabla_{g}f_{j}(x)\rangle_{g}$.
Note that as per our convention $\nabla\phi_{i}(x_{k})=\nabla(\phi_{i}\ \circ\
\text{exp}_{x_{k}})(0)$ and therefore
$\langle\nabla_{g}\phi_{i}(x),\nabla_{g}\phi_{j}(x)\rangle_{g}=\nabla\phi_{i}(x_{k})^{T}\nabla\phi_{j}(x_{k})$.

##  Appendix C Rationale Behind the Choice of $t_{k}$ in Eq. (25)

Since $|\mathcal{M}|\leq 1$, we note that

$\displaystyle\epsilon_{k}\leq\Gamma(d/2+1)^{1/d}/\sqrt{\pi}$ (92)

where the maximum can be achieved when $\mathcal{M}$ is a $d$-dimensional ball
of unit volume. Then we take the limiting value of $t_{k}$ as in Eq. (25)
where chi2inv is the inverse cdf of the chi-squared distribution with $d$
degrees of freedom evaluated at $p$. Since the covariance matrix of
$G(t_{k},x,y)$ is $\sqrt{2t_{k}}I$ (see Eq. (21)), the above value of $t_{k}$
ensures $p$ probability mass to lie in $B_{\epsilon_{k}}(x_{k})$. We take $p$
to be $0.99$ in our experiments. Also, using Eq. (92) and Eq. (25) we have

$\displaystyle
t_{k}\leq\frac{1}{2\pi}\frac{\Gamma(d/2+1)^{2/d}}{\text{chi2inv}(p,d)}<<1,\
\text{when }p=0.99.$ (93)

Using the above inequality with $p=0.99$, for $d=2,10,100$ and $1000$, the
upper bound on $t_{k}=0.0172,0.018,0.0228$ and $0.0268$ respectively. Thus,
$t_{k}$ is indeed a small value close to $0$.

##  Appendix D Computation of $(s_{m},p_{s_{m}})_{m=1}^{M}$ in Algo. 5

Algo. 5 aligns the intermediate views in a sequence. The computation of the
sequences $(s_{m},p_{s_{m}})_{m=1}^{M}$ is motivated by the necessary and
sufficient conditions for a unique solution to the standard orthogonal
Procrustes problem [39]. We start by a brief review of a variant of the
orthogonal Procrustes problem and then explain how these sequences are
computed.

###  D.1 A Variant of Orthogonal Procrustes Problem

Given two matrices $A$ and $B$ of same size with $d$ columns, one asks for an
orthogonal matrix $T$ of size $d\times d$ and a $d$-dimensional columns vector
$v$ which most closely aligns $A$ to $B$, that is,

$\displaystyle
T,v=\mathop{\mathrm{argmin}}\limits_{\Omega,\omega}\left\|A\Omega+\mathbf{1}_{n}\omega^{T}-B\right\|^{2}_{F}\text{
such that }\Omega^{T}\Omega=I.$ (94)

Here $\mathbf{1}_{n}$ is the $n$-dimensional column vector containing ones.
Equating the derivative of the objective with respect to $\omega$ to zero, we
obtain the following condition for $\omega$,

$\displaystyle\omega=\frac{\mathbf{1}_{n}}{n}^{T}(A\Omega-B).$ (95)

Substituting this back in Eq. (94), we reduce the above problem to the
standard orthogonal Procrustes problem,

$\displaystyle
T=\mathop{\mathrm{argmin}}\limits_{\Omega}\left\|\overline{A}\Omega-\overline{B}\right\|_{F}^{2}$
(96)

where

$\displaystyle\overline{X}=\left(I-\frac{1}{n}\mathbf{1}_{n}\mathbf{1}_{n}^{T}\right)X$
(97)

for any matrix $X$. This is equivalent to subtracting the mean of the rows in
$X$ from each row of $X$.

As proved in [39], the above problem, and therefore the variant, has a unique
solution if and only if the square matrix $\overline{A}^{T}\overline{B}$ has
full rank $d$. Denote by $\sigma_{d}(X)$ the $d$th smallest singular value of
$X$. Then $\overline{A}^{T}\overline{B}$ has full rank if
$\sigma_{d}(\overline{A}^{T}\overline{B})$ is non-zero, otherwise there exists
multiple $T$ which minimize Eq. (94).

###  D.2 Computation of $(s_{m},p_{s_{m}})_{m=1}^{M}$

Here, $s_{m}$ corresponds to the $s_{m}$th intermediate view and $p_{s_{m}}$
corresponds to its parent view. The first view in the sequence corresponds to
the largest cluster and it has no parent, that is,

$\displaystyle
s_{1}=\mathop{\mathrm{argmax}}\limits_{m=1}^{M}|\mathcal{C}_{m}|\text{ and
}p_{s_{1}}=\text{none}.$ (98)

For convenience, denote $s_{m}$ by $s$, $p_{s_{m}}$ by $p$ and
$V_{mm^{\prime}}$ by $\widetilde{\Phi}^{g}_{m}(\widetilde{U}_{mm^{\prime}})$.
We choose $s$ and $p$ so that the view $V_{sp}$ can be aligned with the view
$V_{ps}$ without any ambiguity. In other words, $s$ and $p$ are chosen so that
there is a unique solution to the above variant of orthogonal Procrsutes
problem (see Eq. (94)) with $A$ and $B$ replaced by $V_{sp}$ and $V_{ps}$,
respectively. Therefore, an ambiguity (non-uniqueness) would arise when
$\sigma_{d}(\overline{V}_{sp}^{T}\overline{V}_{ps})$ is zero. We quantify the
ambiguity in aligning arbitrary $m$th and the $m^{\prime}$th intermediate
views on their overlap, that is, $V_{mm^{\prime}}$ and $V_{m^{\prime}m}$, by

$\displaystyle
W_{mm^{\prime}}=\sigma_{d}(\overline{V}_{mm^{\prime}}^{T}\overline{V}_{m^{\prime}m}).$
(99)

Note that $W_{mm^{\prime}}=W_{m^{\prime}m}$. A value of $W_{mm^{\prime}}$
close to zero means high ambiguity in the alignment of $m$th and
$m^{\prime}$th views. By default, if there is no overlap between $m$th and
$m^{\prime}$th view then $W_{mm^{\prime}}=W_{m^{\prime}m}=0$.

Finally, we compute the sequences $(s_{m},p_{s_{m}})_{m=2}^{M}$ so that
$\sum_{m=2}^{M}W_{s_{m}p_{s_{m}}}$ is maximized and therefore the net
ambiguity is minimized. This is equivalent to obtaining a maximum spanning
tree $T$ rooted at $s_{1}$, of the graph with $M$ nodes and $W$ as the
adjacency matrix. Then $(s_{m})_{m=2}^{M}$ is the sequence in which a breadth
first search starting from $s_{1}$ visits the nodes in $T$. And $p_{s_{m}}$ is
the parent of the $s_{m}$th node in $T$. Thus,

$\displaystyle(s_{m})_{m=2}^{M}=\text{Breadth-First-Search}(T,s_{1})\text{ and
}p_{s_{m}}=\text{parent of }s_{m}\text{ in }T.$ (100)

##  Appendix E Computation of $\widetilde{U}^{g}_{mm^{\prime}}$ in Eq. (58)

Recall that $\widetilde{U}^{g}_{mm^{\prime}}$ is the overlap between the $m$th
and $m^{\prime}$th intermediate views in the embedding space. The idea behind
its computation is as follows. We first compute the discrete balls $U^{g}_{k}$
around each point $y_{k}$ in the embedding space. These are the analog of
$U_{k}$ around $x_{k}$ (see Eq. 26) but in the embedding space, and are given
by

$\displaystyle U^{g}_{k}=\\{y_{k^{\prime}}\ |\
d_{e}(y_{k},y_{k^{\prime}})<\epsilon^{g}_{k}\\}.$ (101)

An important point to note here is that while in the ambient space, we used
$\epsilon_{k}$, the distance to the $k_{\text{lv}}$th nearest neighbor, to
define a discrete ball around $x_{k}$, in the embedding space, we must relax
$\epsilon_{k}$ to account for a possibly increased separation between the
embedded points. This increase in separation is caused due to the distorted
parameterizations. Therefore, to compute discrete balls in the embedding
space, we used $\epsilon^{g}_{k}$ in Eq. (101), which is the distance to the
$\nu k_{\text{lv}}$th nearest neighbor of $y_{k}$. In all of our experiments,
we take $\nu$ to be $3$.

Recall that $c_{k}$ is the cluster label for the point $x_{k}$. Using the same
label $c_{k}$ for the point $y_{k}$, we construct secondary intermediate views
$\widetilde{U}^{g}_{m}$ in the embedding space,

$\displaystyle\widetilde{U}^{g}_{m}=\cup_{c_{k}=m}U^{g}_{k}.$ (102)

Finally, same as the computation of $\widetilde{U}_{mm^{\prime}}$ in Eq. (53),
we compute $\widetilde{U}^{g}_{mm^{\prime}}$ as the intersection of
$\widetilde{U}^{g}_{m}$ and $\widetilde{U}^{g}_{m^{\prime}}$,

$\displaystyle\widetilde{U}^{g}_{mm^{\prime}}=\widetilde{U}^{g}_{m}\cap\widetilde{U}^{g}_{m^{\prime}}.$
(103)

##  Appendix F Comparison with the Alignment Procedure in LTSA

In the following we use the notation developed in this work. LTSA [49]
computes the global embedding $Y_{m}$ of the $m$th intermediate view
$\widetilde{U}_{m}$ so that it respects the local geometry determined by
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$. That is,

$\displaystyle
Y_{m}=\widetilde{\Phi}_{m}(\widetilde{U}_{m})L_{m}+e_{m}v_{m}^{T}+E_{m}.$
(104)

Here, $Y=[y_{1},y_{2},\ldots,y_{n}]^{T}$ where $y_{i}$ is a column vector of
length $d$ representing the global embedding of $x_{i}$, $Y_{m}$ is a
submatrix of $Y$ of size $|\widetilde{U}_{m}|\times d$ representing the global
embeddings of the points in $\widetilde{U}_{m}$, and
$\widetilde{\Phi}_{m}(\widetilde{U}_{m})$ is a matrix of size