
\section{Introduction}
Manifold learning or nonlinear dimensionality reduction (NLDR) methods play an important role in understanding high-dimensional data across a wide range of disciplines such as computer science \cite{tang2015line}, medicine and biology \cite{becht2019dimensionality, esteva2017dermatologist}, and chemistry \cite{dsilva2013nonlinear}. Understanding the manifold structure underlying a model or the raw data can have significant implications for explanatory aspects of machine learning. For example, it can engender trust in models diagnosing cancer \cite{esteva2017dermatologist} or aid scientific discovery by helping biologists identify novel cell types not only via clustering of RNA sequence data, but also by biologically meaningful properties \cite{rostom2017interpretscRNA, haque2017practical}.  Unfortunately, despite two decades of effort, manifold learning is still not at the point that it can be used systematically in applications.  For many of the most popular methods, there is no theory that provides general conditions under which they yield suitable embeddings, and indeed there are known counterexamples.  Of particular concern, extant theory often requires a strong assumption that there is no noise, and the methods indeed often fail in the presence of noise.  These methods can fail dramatically with even small perturbations to the data, such as a modest rescaling or a modest amount of additive noise~\cite{goldberg2008manifold}.


In the current paper, we provide a new framework, and accompanying theory, that mitigates several of these problems.  The framework comprises two algorithmic components.  The primary component is \emph{Manifold Deflation}, a technique that iteratively reduces the dimension of the underlying manifold learning problem as coordinates are added. Each additional coordinate adds a constraint which ensures the next coordinate yields appropriate local geometry and identifies meaningful directions in the tangent space.  The second component is \emph{Vector Field Inversion}, which removes boundary bias by casting a coordinate as a solution to a regression problem that preserves local geometry rather than an eigenproblem that does not preserve geometric properties.  Working together, these two components provide a way to recover novel embeddings and solve or mitigate a number of failures of manifold learning due to scale, noise, holes, and boundaries.

Our framework is best viewed as a modification rather than a wholesale replacement of existing methods.  In particular, the framework requires the specification of a base NLDR method, which can be any of a wide class of local, spectral methods.  This includes some of the most popular NLDR methods, such as Locally Linear Embedding (LLE) \cite{roweis00LLE}, Laplacian Eigenmaps (LE) \cite{belkin2003laplacian}, Diffusion maps \cite{coifman2006diffusion}, Local tangent space alignment (LTSA) \cite{zhang2004principal}, and Hessian LLE (HLLE) \cite{HessEig}.  Used alone, these methods can and do fail, but, as we show, many of the failures can be attributed to a failure to ensure that output embedding coordinates jointly yield suitable local geometries, despite using statistics about local neighborhoods as inputs to generate the embedding.  Our framework provably avoids these failures.










Both Manifold Deflation and Vector Field Inversion are based on a recent dual characterization of local spectral methods, first in terms of a differential operator and a set of boundary conditions and second as a linear smoother \cite{tingnldr}. The original differential operator can be expressed in local coordinate systems as derivatives on $m$ variables where $m$ is the intrinsic dimension of the manifold. Manifold deflation removes one variable as each coordinate is added, thus ``deflating'' the problem into an $(m-1)$-dimensional problem.

Performing this deflation operation requires estimating fundamental objects for a manifold: the tangent vector field induced by a coordinate function and the tangent space at every point.  These are  inferential problems, and we accordingly provide statistical methodology for solving them, by showing how the tangent vector field can be represented as a linear smoother and estimated using local linear regression. We also provide statistical theory, establishing the consistency of this estimator.







\subsection{Related work}
Our methods augment any local, spectral NLDR method \cite{tingnldr}. As we have noted, many of the most popular NLDR methods fall in this class \cite{roweis00LLE, coifman2006diffusion, belkin2003laplacian, zhang2004principal, HessEig}, but there are other methods that do not.  These include Isomap \cite{bernstein2000graph}, t-SNE \cite{maaten2008tsne}, Maximum Variance Unfolding (MVU) \cite{weinberger2006unsupervised}, and UMAP \cite{mcinnes2018umap}.
 
We draw inspiration from work by \cite{goldberg2008manifold} and \cite{gerber2007robust}, who isolated a key structural problem underlying manifold learning algorithms: the problem of scale and repeated eigendirections. Recent work has attempted to address these problems by skipping eigenvectors or by ensuring the next coordinate cannot be predicted from other coordinates~\cite{chen2019selecting, blau2017non-redundant, dsilva2018parsimonious}.   Alternatively, \cite{gerber2007robust} proposed mapping an $m$-dimensional manifold to a lower dimensional manifold and recomputing the Laplacian.

There has been little work addressing boundary bias problems and non-geodesically convex manifolds with holes. Theoretical analyses of methods such as Isomap and MVU \cite{arias2013convergence} assume convexity. While HLLE \cite{HessEig} and the asymptotically equivalent LTSA \cite{zhang2017equivalence, tingnldr} avoid the convexity assumption and boundary bias when there is no noise, their behavior degenerates when there is noise (see Figure \ref{fig:scurve}). Work on handling noise off the manifold in NLDR methods focus on denoising the manifold in the ambient space \cite{hein2007manifold, fefferman2018fitting} or treating the problem as one of outlier detection \cite{chen2006robust}.












\section{Preliminaries}
\label{sec:preliminaries}
Consider an $m$-dimensional compact Riemannian manifold $\ensuremath{\mathcal M}$  with smooth boundary smoothly embedded in a high-dimensional ambient space $\mathbb{R}^d$ with $d \geq m$. All of the manifolds we consider in this paper will be assumed to satisfy these properties.
Given $n$ samples $y^{(1)}, \ldots, y^{(n)}$ from $\ensuremath{\mathcal M}$, potentially with noise in $\mathbb{R}^d$, we wish to find a good low-dimensional representation of $\ensuremath{\mathcal M}$. 
This representation is a mapping $\phi$ from $\ensuremath{\mathcal M}$ into $\mathbb{R}^{d'}$ that defines a set of global coordinate functions, $\phi_1(p), \ldots \phi_{d'}(p)$, where $p \in \ensuremath{\mathcal M}$. 

The main relevant concepts from differential geometry are the relationships among tangent vectors, differential operators, coordinate charts, and geodesic curves passing through a point. We briefly review these relationships.
Associated with each point $p\in \ensuremath{\mathcal M}$ is a tangent vector space $T_p(\ensuremath{\mathcal M})$. Each unit tangent vector $v \in T_p(\ensuremath{\mathcal M})$ is associated with a unique one-dimensional geodesic curve with unit geodesic length passing through $p$ in the direction $v$. In a local neighborhood, a geodesic curve is the unique shortest path connecting two points $p, q$. The point $q$ can be mapped to the vector $v$ tangent to the geodesic curve at $p$ with magnitude $\|v\|$ equal to the geodesic length. Expressing $v$ in a basis of $T_p(\ensuremath{\mathcal M})$ gives a normal coordinate chart for points in a neighborhood of $p$.
 Finally, a tangent vector $v$ naturally defines a differential operator by the partial derivative $f(p) \mapsto \partial f(p) / \partial v$ in the normal coordinate chart. This defines a bijective relationship between tangent vectors and partial derivatives.

Given tangent vectors $v_p$ for all $p \in \ensuremath{\mathcal M}$, a vector field $V$ can thus be defined as the mapping $V : C^\infty(\ensuremath{\mathcal M}) \to C^\infty(\ensuremath{\mathcal M})$ such that $V f(p) = \partial f(p) / \partial v_p$. Since derivatives are linear operators, finite approximations for $V$ can be represented as a matrix. Our statistical methods are based on estimating and exploiting this matrix. Figure \ref{fig:diff geom} illustrates these relationships. 


\begin{figure}
    \hspace{-0.6cm} 
    \includegraphics[width=6in]{figures/Manifold_deflation_diagram.png}
    \caption{Relationships among differential geometry concepts. Pink denotes the areas our methods uniquely contribute to. These close the loop and allow for iterative refinement of coordinates.}
    \label{fig:diff geom}
\end{figure}

\subsection{Problems}
Given samples from an $m$-dimensional manifold $\ensuremath{\mathcal M}$, potentially corrupted with high-dimensional noise to form $\ensuremath{\mathcal M} + \mathcal{E}$, our goal is to generate an embedding with
coordinate functions $\phi_i : \ensuremath{\mathcal M} + \mathcal{E} \to \ensuremath{\mathcal M}' \subset \mathbb{R}^m$ that preserves tangent spaces in the following ways. 
The tangent vectors $\nabla \phi_i(p)$ form an orthogonal basis of $T_p(\ensuremath{\mathcal M})$,
and furthermore the basis is orthonormal. If these conditions are perfectly satisfied at all points $p \in \ensuremath{\mathcal M}$, then the embedding is an \emph{isometric embedding}. If only orthogonality is satisfied for all points, then one obtains an \emph{immersion}. As we will see, the latter weakening addresses the practical problem of how to generate good embeddings when local, spectral methods fail and permits discovery of novel, interesting embeddings that reveal the manifold's structure even when an isometric embedding is impossible to attain.


Situations under which local, spectral methods fail include simple problems of scale, non-convexity, and noise. These are illustrated in Figure \ref{fig:scurve}. 
If a manifold has a large aspect ratio---that is, it has underlying coordinates where the range of one coordinate is much longer than the others---multiple eigenvectors may encode the same coordinate direction. This problem is called the \emph{repeated eigendirections problem}~\cite{goldberg2008manifold, gerber2007robust}. Even on a simple two-dimensional s-curve, this problem leads to failures for Laplacian Eigenmaps.

Embeddings can also be changed dramatically by boundaries and holes. 
Some methods, such as LTSA, require estimating local tangent spaces and knowing the intrinsic dimension of the manifold. When a small amount of off-manifold additive noise is added, we prove in the supplementary material that it is impossible to estimate the tangent space when neighborhood sizes shrink to zero. Points are nearly uniformly distributed in a small ball so that meaningful tangent space directions cannot be differentiated from noise directions. Figure \ref{fig:scurve} shows even a small amount of noise causes LTSA to behave like the regular Laplacian. 






\begin{figure}
    \hspace{-1.5cm}
    \begin{tabular}{c|c}
    \raisebox{1.2cm}[0pt][0pt]{\rotatebox{90}{Scale}} &  \hspace{-0.4cm}\includegraphics[width=6in,height=1.0in]{figures/scurve.png} \\
    \raisebox{0.3cm}[0pt][0pt]{\rotatebox{90}{Non-convexity}} & \hspace{-0.4cm}\includegraphics[width=6in,height=1.0in]{figures/scurve_w_offset_hole.png}  \\ 
    \raisebox{1.2cm}[0pt][0pt]{\rotatebox{90}{Noise}}  & \hspace{-0.4cm}\includegraphics[width=6in,height=1.0in]{figures/scurve_w_offset_hole_noise.png} 
    \end{tabular}
    \caption{Scurve}
    \label{fig:scurve}
\end{figure}


\section{Estimators and Embeddings using Differential Operators}
We first briefly review local, spectral methods and their connections to partial differential equations (PDEs) on $m$ variables. We then show manifold deflation corresponds to eliminating one variable at a time from the PDE and describe the geometric connection to orthogonal tangent vectors.

Local, spectral methods consist of three main steps: 1) constructing a neighborhood graph, 2) constructing a sparse matrix $L^{(n)}$ from the local neighborhoods, and 3) solving an eigenproblem using $L^{(n)}$, where the eigenvectors correspond to the embedding. 
Recent work shows that for any local, spectral method, if $L^{(n)}$ converges as neighborhoods shrink in radius, then it must converge to a differential operator and the domain of the operator is specified by a set of boundary conditions~\cite{tingnldr}. For example,  the negative weighted Laplace-Beltrami operator is the limit differential operator for Laplacian Eigenmaps \cite{belkin2003laplacian} and Diffusion Maps \cite{coifman2006diffusion}. Concretely, when points are sampled uniformly from the manifold, the following unweighted Laplace-Beltrami operator is the limit operator, and the domain imposes a Neumann boundary condition:
\begin{align}
 \label{eqn:laplace-beltrami}
    -L^{(n)} f \to \Delta f &= \sum_{i=1}^{m} \frac{\partial^2 f}{\partial x_i^2},\quad \mbox{for}\quad f \in dom(\Delta) = \left\{g : \frac{\partial g}{\partial \nu_p} = 0, \quad\forall p \in \partial \ensuremath{\mathcal M} \right\}.
 \end{align} 
 Here, $\nu_p$ is a tangent vector normal to the boundary at $p \in \ensuremath{\mathcal M}$. The partial derivatives are taken with respect to local, normal coordinate charts, and convergence is pointwise under the uniform norm. 
Thus, the limit eigenproblem yields estimated coordinates $\hat{\phi} \in dom(\Delta)$ that satisfy the PDE on $m$ variables, $\sum_{i=1}^{m} \frac{\partial^2 \hat{\phi}}{\partial x_i^2} = \lambda \phi$, for some $\lambda$ in a neighborhood with normal coordinates denoted by $\mathbf{x}$. %

\subsection{Manifold deflation}
\label{sec:manifold deflation}
The main idea behind manifold deflation is to eliminate one coordinate from the differential operator. The solutions of the modified differential equations are then functions of the remaining coordinates, thus ``deflating'' the problem from one on $m$ dimensions to one on $m-1$ dimensions. 
 









 We do this by restricting the domain of the operator. Given a coordinate function $\phi_k(p)$, then, with a slight abuse of notation, we impose the constraint that $\partial f(p) / \partial \phi_k(p) = 0$. In Euclidean space this encodes the fact that in Cartesian coordinates the coordinate axes are perpendicular to each other. Correspondingly, the gradients of the coordinate functions are orthogonal, and the directional derivatives satisfy $\partial \phi_j(p) / \partial \phi_k = 0$ whenever $j \neq k$.
 Imposing this constraint forces any term in the differential operator containing the coordinate $\phi_k$ to be zero, so it can be eliminated. 
 
 To formally define the constraint, note that the gradients $\nabla \phi_k(p)$ generate a collection of tangent vectors on $\ensuremath{\mathcal M}$. From Section \ref{sec:preliminaries}, each tangent vector corresponds to a directional derivative in a local normal coordinate system. This defines a vector field $\mathbf{V}_k$ mapping a function to its directional derivative at each point.  Thus, the constraint can be expressed as $\mathbf{V}_k f = 0$. 
 We can now formalize our notation, by defining
 $\partial f(p) / \partial \phi_k := \|\nabla \phi_k(p)\|^{-2} V_k(p) f$
 as a rescaling of $V_k$.
 
Given a local, spectral method that converges to a positive semidefinite operator $L^\infty$, manifold deflation iteratively solves the following:
 \begin{align}
     \phi_{k+1} &= \argmin_{\substack{f \in dom(L^\infty) \cap \Omega_k \\ 
                                    \langle f, 1  \rangle = 0, \, 
                                    \| f \| = 1}}
                        f^T L^\infty f \quad \mbox{where }\quad 
    \Omega_i = \left\{f \in C^{\infty}(\ensuremath{\mathcal M}): V_j f = 0\quad \forall j \leq k\right\}.
\end{align}
 
We note the similarity to matrix deflation methods in numerical linear algebra. The only difference is that matrix deflation imposes an orthogonality constraint $\Omega_k = (\mathrm{span} \{\phi_j : j \leq k\})^\perp$  instead of the vector field constraint. The former are constraints in a function space while the latter are constraints in local neighborhoods and have a geometric interpretation.
The hard constraints on the derivative can be expressed as soft constraints as follows: 
\begin{align}
\label{eqn:soft constraints}
        \phi_{k+1}^{(\lambda)} &= \argmin_{\substack{f \in dom(L^\infty) \cap \Omega_k \\ 
                                    \langle f, 1  \rangle = 0, \, 
                                    \| f \| = 1}}
                        \langle f, L^\infty f\rangle + \lambda \sum_{j \leq i}\| V_j f \|_{L_2(\nu)}^2,
\end{align}
for any measure $\nu$ on $\ensuremath{\mathcal M}$ with density bounded away from 0 and $\infty$.
This yields the algorithm:
\begin{enumerate}[noitemsep]
    \item Take as input $L$ from any local, spectral method and a regularization parameter $\lambda$. 
    \item Compute the bottom, non-constant eigenvector $\phi_k$ of $L + \lambda \sum_{j=1}^{k-1} \hat{V}_j^T \hat{V}_j$.
    \item Add $\phi_k$ to the list of coordinates and estimate the matrix valued vector field $\hat{V}_k$ from $\phi_k$.
    \item Increment $k$ and repeat from step 2 until $k = dim(\ensuremath{\mathcal M})$.
\end{enumerate}
    
\subsection{Vector field and tangent space estimation}
The main challenge in performing manifold deflation is estimating the vector field $\hat{V}_k$.
We show that this estimation problem can be reduced to a local linear regression, yielding our proposed Tangent Derivation Estimator.



 
 

First, we examine derivative estimation when a coordinate system is known, which we can solve by extending techniques from local polynomial regression \cite{fan1996localpolyregression}. Consider a Taylor expansion in a local, normal coordinate system for a ball $B(p_0,h)$. A point $q$ in the neighborhood has coordinates $x_1(q), \ldots, x_m(q)$, so that
 \begin{align}
     f(q) &= f(p_0) + \sum_{i=j}^{m} \frac{\partial f(p_0)}{\partial x_j}x_j(q) + \epsilon_q,
 \end{align}
 where $\epsilon_q = O(\|h\|^2)$.
Given points $\ensuremath{\mathcal N} = \{q_1, \ldots q_\ell\}$ in a small ball around $p_0$ on the manifold and a function $f$, the partial derivative with respect to $x_i$ can be estimated by a local linear regression where the $\ell \times (m+1)$ matrix of regressors is given by the coordinates $M_{ai} = x_i(q_a)$ and the intercept $M_{a,m+1} = 1$. Writing $\beta_k = \frac{\partial f(p_0)}{\partial x_k}$, the derivative estimate is the regression coefficient $\hat{\beta}_k(f) = e_{k}^T (M^T M)^{-1} M^T S f_\ensuremath{\mathcal N} \approx \frac{\partial f(p_0)}{\partial x_j}$, where $f_\ensuremath{\mathcal N} = (f(q_1), \ldots, f(q_\ensuremath{\mathcal N}))$.

Extending this derivative estimator to the case of manifold deflation requires computing an estimate when only a subset of $k$ coordinates out of $m$ are known. We show in the supplementary material that the unknown coordinates can be ignored
when the sampled neighborhoods are sampled uniformly from symmetric balls contained in the interior of $\ensuremath{\mathcal M}$. In this case, the normal coordinate functions in the neighborhood asymptotically have zero mean and are orthogonal. Since the unknown coordinates' regressors are orthogonal to the $k^{th}$ coordinate's regressor, they can be ignored. This also implies that other known coordinates can be ignored and the derivative can be estimated using a simple linear regression that depends only on the $k^{th}$ coordinate. 

We now define the Tangent Derivation Estimator.
Let $\mathcal{A}$ be the indices for the $a^{th}$ point's neighbors and let $\phi_{k, \mathcal{A}}$ denote the values of the coordinate of interest, $\phi_k$, evaluated at the $\mathcal{A}$ neighbors. Let $\overline{\phi}_{k, \mathcal{A}}$ be the mean of those values.
A simple linear regression in each neighborhood gives that vector field can be estimated with an $n \times n$ matrix $\hat{V}_k$, where
\begin{align}
    \hat{V}_{k, a, \ensuremath{\mathcal A}} &= (\phi_{k,\ensuremath{\mathcal A}} - \overline{\phi}_{k,\ensuremath{\mathcal A}}) / \|\phi_{k,\ensuremath{\mathcal A}} - \overline{\phi}_{k,\ensuremath{\mathcal A}}\|^2. \label{eqn:vf approx}
\end{align}
Here $\hat{V}_{k, i, J}$ denotes the entries indexed by $J$ in the $i^{th}$ row. This estimator has the consistency property given in Theorem \ref{thm:derivative consistency}. The derivation of the estimator given above sketches the proof that it converges to the correct limit.





























\begin{theorem}[Consistency of vector field estimate]
\label{thm:derivative consistency}
    Let $\ensuremath{\mathcal Y}^{(n)}$ be a uniform sample from a smooth $m$-dimensional manifold $\ensuremath{\mathcal M} \subset \mathbb{R}^d$ and let the sequence
    $(\epsilon_n)$ decrease to zero. 
    Let  $\phi : \ensuremath{\mathcal M} \to \mathbb{R}$ be a smooth function such that $\inf_{p\in \ensuremath{\mathcal M}} \|\nabla \phi(p)\| > 0$. Given a smooth function $f$ with bounded gradient.
    Let $\hat{D}_{\phi,\epsilon_n}$ be the estimate of $\partial f / \partial \phi$ given by Eq.~\ref{eqn:vf approx} on an $\epsilon_n$-neighborhood graph.
    Then, for some constant $\alpha$ and for any $h > 0$, if $\epsilon_n \to 0$ sufficiently slowly as $n \to \infty$,
    \begin{align}
    \label{eqn:consistency}
        \sup_{\substack{p \in \ensuremath{\mathcal M} \\ d(p, \partial \ensuremath{\mathcal M}) > h}} \left|\frac{\partial f}{\partial \phi}(p) - \alpha \hat{D}_{\phi,\epsilon_n} (p)\right| \ensuremath{\stackrel{p}{\to}} 0.
    \end{align}
\end{theorem}

Since for any $\epsilon > 0$, the estimator $\hat{D}_{\phi,\epsilon}$ is a continuous function of $\phi$ under the uniform norm, the continuous mapping theorem extends the theorem to estimated coordinates $\phi^{(n)}$ that converge to  $\phi$.
\begin{corollary}
 Let $\phi^{(n)}$ be a sequence of smooth, real valued functions such that 
 $\phi^{(n)} \to \phi$ under the uniform norm. 
 Then, under the conditions of theorem \ref{thm:derivative consistency}, Eq. \ref{eqn:consistency} holds with 
 $\hat{D}_{\phi,\epsilon_n}$ replaced with $\hat{D}_{\phi_n,\epsilon_n}$.
\end{corollary}
 
 
 

 
 \subsection{Example: Laplacian Eigenmaps}
 To illustrate manifold deflation, consider Laplacian Eigenmaps on a uniform sample drawn from a 3D strip,
$[0,9\pi] \times [0, 3\pi] \times [0, \pi] \subset \mathbb{R}^3$. Ignoring the non-smoothness at the corners, this method asymptotically generates the Laplace-Beltrami operator given in Eq.~\ref{eqn:laplace-beltrami}, with boundary conditions
 \begin{align}
 \label{eqn:laplace-beltrami boundary}
    \frac{\partial f(p)}{\partial x_i} = 0 \quad \mbox{if } p_i \in \{0,b_i \pi\} \quad 
    \mbox{where $b_1 = 9, b_2 = 3, b_3=1$}.
 \end{align} 
 

The resulting PDE is separable, such that each coordinate can be treated as an ordinary differential equation (ODE), $-\frac{\partial^2 f}{\partial x_i^2} = \lambda f$, with $\frac{\partial f(0)} {\partial x_i} = \frac{\partial f(b_i)}{\partial x_i} = 0$. The eigenfunctions for this ODE are given by $\cos\left( j x_1 / b_i\right)$ 
for $j \in \mathbb{N}^0$. %

The resulting eigenfunctions for the PDE include these eigenfunctions and their products. The first three nonzero eigenvalues and their corresponding eigenfunctions are given by

\begin{tabular}{lllll}
Eigenvalue & $1/9^2$ & $(2/9)^2$ & $1/3^2$ & $1/3^2 + 1/9^2$  \\ 
Eigenfunctions & $\cos\left( \frac{x_1}{9} \right)$ & $\cos\left( \frac{2x_1}{9} \right)$ & $\cos\left( \frac{x_1}{3} \right), \, \cos\left( \frac{x_2}{3} \right)$ &
$\cos\left( \frac{x_1}{9} \right) \cos\left( \frac{x_2}{3} \right)$ 
\end{tabular}

Thus, the Laplacian Eigenmaps method chooses a function of only $x_1$ for the first three coordinates. A further computation shows that the first 17 eigenfunctions will be solely functions of $x_1$ and $x_2$, so that $x_3$ is not picked up.
Note that the second time $x_2$ appears is in the cross-term
$\cos(x_1/9) \cos(x_2/3)$. This has an eigenvalue that is close to the eigenvalue of $\cos(x_2/3)$ which can cause commonly observed mixing of the eigenvectors when points are sampled randomly. In this case, methods using a predictability heuristic~\cite{blau2017non-redundant, dsilva2018parsimonious} or eigenvector skipping \cite{chen2019selecting} may fail to appropriately distinguish the eigenvectors. We are unaware of work that addresses this form of repeated eigendirections problem. 
On the other hand, manifold deflation eliminates all other eigenfunctions containing $x_1$, including the cross-term $\cos(x_1/9) \cos(x_2/3)$. If the  coordinate $x_3$ with the shortest range is considered a noise dimension, manifold deflation yields a two-dimensional embedding that ignores this noise.













\subsection{Boundary bias and the Vector Field Inversion embedding}
The estimated vector fields can also be used outside of the manifold deflation procedure. We demonstrate one use in removing boundary bias.
Methods that impose a Neumann boundary condition, such as Laplacian Eigenmaps, cannot recover an isometric embedding due to boundary bias~\cite{tingnldr}. 
Trivially, every coordinate function in an isometric embedding must have a gradient with unit norm everywhere. However, the Neumann boundary condition implies the gradient must be zero at some point on the boundary. 

Suppose there exist global coordinate functions $\phi_1, \ldots, \phi_m$ that generate an isometric embedding of the manifold $\ensuremath{\mathcal M}$.
The vector field $V_k := \partial / \partial \phi_k$ satisfies $V_k \phi_k = 1$ where $1$ denotes the constant function equal to one. 
To obtain a bias-corrected coordinate, assume $\phi_k$ lies in some reproducing kernel Hilbert space (RKHS). Applying kernel ridge regression to the problem $1 = V_k \phi_k + \epsilon$ yields:
\begin{align}
\label{eqn:debiasing}
    \tilde{\phi}_k &= K (K V_k^T V_k K + \alpha I)^{-1} K V_k^T 1,
\end{align}
for some kernel $K$ and regularizer $\alpha$.
We note that naive matrix inversion of the finite matrix approximation to $V_k$ will not yield meaningful results. We posit this is because the continuous, limit operator $V_k$ is not invertible, since the inverse is not bounded.


\vspace{-0.2cm}
\subsection{Space and time complexity}
Manifold deflation consists of two phases: 1) estimating a vector field $V_k$ and penalty matrix $V_k^T V_k$, and 2) estimating the next coordinate. 
The vector field estimator given in Eq.~\ref{eqn:vf approx} takes time $O(|E|)$ as it simply copies and rescales values from the coordinate into a sparse matrix with $|E|$ nonzero entries where $E$ are the edges in the neighborhood graph.
The number of nonzero entries in the penalty is equal to the number of node pairs two hops or less away in the neighborhood graph. Note there is no dependence on the dimension of the original data or on the intrinsic dimension for the steps unique to manifold deflation. Computing the second smallest eigenvector of a positive semi-definite sparse matrix can be done efficiently with iterative algorithms such as LOBPCG \cite{Knyazev2001TowardTO} or the Lanczos method with inverse iterations.





\vspace{-0.1cm}
\section{Experiments}
We demonstrate the properties of manifold deflation and vector field inversion with experiments on three datasets: an S-curve, a sphere, and the Fashion MNIST dataset. 
In each experiment, we use Laplacian Eigenmaps with an unnormalized Laplacian as the base method for manifold deflation and k-nearest neighbors with 15 neighbors to define the adjacency graph. Additional details about the experiment as well as additional results are included in the supplementary material. 

The S-curve is an isometric embedding of a $3 \times 1$ rectangle in $\mathbb{R}^2$ into an S-shaped curve in $\mathbb{R}^3$. Figures \ref{fig:scurve} shows the behavior of NLDR methods when various challenges are incorporated.  When there is a hole or additive noise, only manifold deflation preserves the two-dimensional structure on the S-curve. Other methods fail despite the small magnitude of the noise, which is additive noise drawn uniformly on a cube with half-width $0.1$.  Non-redundant spectral embedding or methods that skip eigenvectors to deal with the repeated eigendirections problem yield worse embeddings and cannot deal with holes. The supplementary material shows that other combinations of eigenvectors still yield poor embeddings so any method that selects a subset of eigenvectors does poorly.

Due to its use of Laplacian Eigenmaps, manifold deflation exhibits boundary bias in figure \ref{fig:debias} when embedding an S-curve with additive noise. Vector field inversion debiases the coordinates near the boundary to yield strips with uniform width. 

Manifold deflation can also recover embeddings that are useful and qualitatively different from other methods. Figure \ref{fig:sphere} shows that it recovers a rescaling of polar coordinates on a sphere. Points on the sphere are selected from a spherical Fibonacci lattice \cite{hannay2004fibonacci} which deterministically generates uniformly spaced points on a sphere. The north-south and east-west axes are  slightly stretched to make the order of eigenvectors deterministic. Thus, if applied to the globe, manifold deflation approximately recovers a rectangular map. By comparison, all other local, spectral methods do not yield meaningful dimensionality reduction and simply return  a linear projection.

Figure \ref{fig:fmnist} shows embeddings of the Fashion MNIST dataset. Laplacian Eigenmaps essentially collapses the global structure into a collection of one dimensional manifolds. Manifold deflation generates an embedding with interesting manifold structure. The figure shows the embedding on different types of shoes. 
For Laplacian Eigenmaps, it is unclear what can learned from the visualization. Manifold deflation, on the other hand, clearly places darker shades to the left and sandals to the right. Shoes with high tops are close to the bottom while those with low profiles are near the top. The supplementary material provide additional plots comparing manifold deflation to Laplacian Eigenmaps.
We note that while methods such as t-SNE and UMAP provide visualizations that are adept at uncovering the cluster structure of the data, manifold deflation yields insight into the continuous manifold structure of the data.



\begin{figure}
\vspace{-1cm}
\hspace{-1cm}
\begin{tabular}{cc}
    \hspace{-0.2cm}\includegraphics[width=3.0in, height=2.9in]{figures/FMNIST_LE_shoes.png} &
    \hspace{-0.3cm}\includegraphics[width=3.3in, height=2.9in]{figures/FMNIST_MD_shoes.png} \\
    \hspace{-0.3cm}\includegraphics[width=3.0in,height=1.5in]{figures/FMNIST_LE_raw.png} &
    \hspace{-0.3cm}\includegraphics[width=3.3in,height=1.5in]{figures/FMNIST_MD_raw.png}
    \end{tabular}
    \caption{Embedding of the Fashion MNIST dataset colored by category (bottom) and restricted to footwear (top). Manifold deflation avoids collapsing the manifold and reveals interesting variation.}
    \label{fig:fmnist}
\end{figure}


\begin{figure}
    \hspace{-1cm}\includegraphics[width=2.1in, height=1.4in]{figures/debiased_coordinate.png}
    \includegraphics[width=4.2in, height=1.4in]{figures/debiased_plot.png}
        \caption{The estimated tangent vector field can be used to debias the coordinate estimates from manifold deflation on an S-curve with noise. The left plots show the first estimated coordinate function versus the true underlying coordinate. Vector field inversion ensures that this relationship is linear; accordingly, the strips in the right plots are of even width.}
    \label{fig:debias}
\end{figure}


\begin{figure}
    \hspace{-0.5cm}\includegraphics[width=6in, height=1.2in]{figures/sphere.png}
        \caption{When embedding the surface of a sphere in $\mathbb{R}^2$, manifold deflation recovers polar coordinates (latitude, longitude) for each hemisphere. Laplacian Eigenmaps and other local, spectral methods linearly project into two dimensions. The left plots show slices of longitude values map roughly to vertical slices in the manifold deflation embedding. The right plots show that each recovered coordinate is a monotone transformation of the corresponding polar coordinate.} 
    \label{fig:sphere}
\end{figure}

\vspace{-0.1cm}
\section{Conclusions}
We have proposed a novel framework for manifold learning and nonlinear dimensionality reduction based on Manifold Deflation and Vector Field Inversion. We demonstrate the effectiveness of this framework by successfully recovering a manifold's structure in cases where a local, spectral method will provably fail. 
The common thread among the components of our framework is that they encode information or enforce geometric constraints on coordinates in the coordinate-free language of differential operators that are the foundation of local, spectral NLDR methods. By doing so, they iteratively refine the differential operator of a local, spectral NLDR method as coordinates are sequentially incorporated into an embedding. We crucially develop a key component that enables this expression as differential operators by estimating vector fields and components of the tangent space using the Tangent Derivation Estimator. We have exhibited specific algorithmic applications of this paradigm, but we note that the paradigm may apply much more widely. 




\section*{Broader Impact}
We see this work's impact to be a step towards improved visualizations and understanding of high dimensional data and complex machine learning models.
Both of these areas can have significant positive societal impact, but we are unable to anticipate if there are any negative consequences. 
We also believe tangent space estimation may have impact beyond machine learning, such as in mathematics research.


\bibliographystyle{plainnat}
