# Trustworthy Dimensionality Reduction
A dissertation report presented for the completion of
the degree of Master of Statistics (M. Stat.) from Indian Statistical
Institute

Subhrajyoty Roy

(June 25, 2021)

###### Abstract

Different unsupervised models for dimensionality reduction like PCA, LLE,
Shannon’s mapping, tSNE, UMAP, etc. work on different principles, hence, they
are difficult to compare on the same ground. Although they are usually good
for visualisation purposes, they can produce spurious patterns that are not
present in the original data, losing its trustability (or credibility). On the
other hand, information about some response variable (or knowledge of class
labels) allows us to do supervised dimensionality reduction such as SIR, SAVE,
etc. which work to reduce the data dimension without hampering its ability to
explain the particular response at hand. Therefore, the reduced dataset cannot
be used to further analyze its relationship with some other kind of responses,
i.e., it loses its generalizability.

To make a better dimensionality reduction algorithm with a better balance
between these two, we shall formally describe the mathematical model used by
dimensionality reduction algorithms and provide two indices to measure these
intuitive concepts such as trustability and generalizability. Then, we propose
Localized Skeletonization and Dimensionality Reduction (LSDR) algorithm which
approximately achieves optimality in both these indices to some extent. The
proposed algorithm has been compared with the state of the art algorithms such
as tSNE and UMAP and is found to be better overall in preserving global
structure while retaining useful local information as well. We also propose
some of the possible extensions of LSDR which could make this algorithm
universally applicable for various types of data similar to tSNE and UMAP.

## Notations

$x$ | A fixed scalar variable
---|---
$\boldsymbol{x}$ | A fixed scalar vector
$\boldsymbol{X}$ | A fixed scalar matrix
$\mathsf{x}$ | A scalar random variable
$\boldsymbol{\mathsf{x}}$ | A random vector
$\boldsymbol{\mathsf{X}}$ | A random matrix
$\mathbb{R}$ | Set of real numbers
$\boldsymbol{\mathsf{1}}_{n}$ | The $(n\times 1)$ column vector of all ones
$\mathbb{E}(\cdot)$ | Expectation operator
$\text{Var}(\cdot)$ | Variance operator
$\text{Cov}(\cdot)$ | Covariance operator

###### Contents

  1. 1 Introduction
    1. 1.1 Mathematical Setup of Dimensionality Reduction
    2. 1.2 Existing Dimensionality Reduction Algorithms
      1. 1.2.1 Metric Preserving Algorithms
      2. 1.2.2 Graphical Algorithms
      3. 1.2.3 Supervised Dimensional Reduction Algorithms
    3. 1.3 Existing Metrics
    4. 1.4 Contributions of this study
  2. 2 A Generalized Framework
    1. 2.1 The Framework for Dimensionality Reduction
    2. 2.2 Dimensionality reduction algorithms
    3. 2.3 Nonparametric estimator of out-of-the-sample reduction
    4. 2.4 Building the reconstruction function
    5. 2.5 Trustability and Consistency index
  3. 3 Localized Skeletonization and Dimensionality Reduction (LSDR)
    1. 3.1 Manifold Approximation Graph
    2. 3.2 Skeletonization
    3. 3.3 Multidimensional Scaling and Kernel embedding
  4. 4 Applications
    1. 4.1 Simulation Studies
      1. 4.1.1 Performance indices
      2. 4.1.2 Performance of LSDR
  5. 5 Conclusions and Future Aspects
  6. A Prerequisites
    1. A.1 Reproducible Kernel Hilbert Spaces (RKHS)
    2. A.2 Vector Valued RKHS
    3. A.3 Procrustes Analysis
  7. B Proofs
    1. B.1 Proof of Theorem 1
    2. B.2 Proof of Theorem 2
    3. B.3 Proof of Theorem 3
    4. B.4 Proof of Theorem 4
    5. B.5 Proof of Theorem 5
    6. B.6 Proof of Theorem 6
    7. B.7 Proof of Theorem 8
    8. B.8 Proof of Theorem 9
  8. C More Simulated Examples

##  1 Introduction

During the current age of big data, data sources have become easily
accessible, leading to an abundance of available real-world data to solve any
practical problems. This abundance is mostly related to the sheer increase in
the number of features or variables about a particular sample, rather than the
increase in the number of sample points. However, analysing these high-
dimensional datasets in their entirety can pose various challenges. Because of
the curse of dimensionality, the samples become sparse in high dimensions
leading to the overfitting of many models. Also, most of these high-
dimensional features are not linked to the response variable under study, and,
one cannot use visualisation techniques to see the common patterns to obtain
insight into the variables. Therefore, it is natural to have a reduced
representation of these high dimensional datasets which should not compromise
the subsequent analysis. Dimensionality Reduction is such a transformation of
high-dimensional data into a meaningful representation of reduced
dimensionality [51]. Ideally, this reduced dimension should be equal to the
intrinsic dimension of the data, which is described as the minimum number of
parameters required to account for the variability presented in the observed
data [15]. For example, it is well known that only $(k-1)$ variables are
enough to describe the segregation between $k$ well-separated classes of any
dimension [15].

There are various applications of dimensionality reduction in different fields
of science. Compression of digital signals from audio, video and images [22],
analysis of fMRI scans, gene expression, protein-coding datasets in
biostatistics [47], finding word embeddings in Natural Language Processing
[21], stock market analysis in Finance [6, 55] are some of the real-life use-
cases for dimensionality reduction techniques.

Dimensionality reduction algorithms can broadly be classified into three
groups based on their focus and approach.

  1. 1.

Metric preserving Algorithms focus on reducing any high dimensional data by
preserving local distances or increasing global variability. PCA, Kernel PCA,
Maximum variance unfolding, Hessian Map, Multidimensional scaling, tSNE etc.
are examples of such algorithms. Neural network models such as autoencoders
try to minimize the reconstruction error which is also a metric in high
dimensional space.

  2. 2.

Graphical Algorithms focus on reducing the dimensionality by approximating the
geodesic distance based on the creation of neighbourhood graphs. Isomap, UMAP,
Locally linear embedding, Laplacian eigenmaps etc. are examples of such
algorithms.

  3. 3.

Supervised Algorithms focus on reducing the dimensionality of predictor
variables so that the ability to explain the response variable is not reduced.
Sliced Inverse Regression (SIR), Sliced Average Variance Estimation (SAVE),
Mean Average Variance Estimation (MAVE) are some of the examples of this
nature.

###  1.1 Mathematical Setup of Dimensionality Reduction

The mathematical setup of dimensionality reduction starts with an $n\times p$
size dataset $\boldsymbol{X}$, where each row represents a sample
$\boldsymbol{x}_{i}\in\mathbb{R}^{p}$, for $i=1,2,\dots n$. For the sake of
simplicity, one assumes that these datapoints
$\boldsymbol{x}_{1},\dots\boldsymbol{x}_{n}$ are independent and identically
distributed observations from a distribution $F$. To make dimensionality
reduction meaningful, it is also assumed that the support of this distribution
$F$ has an intrinsic dimension $d$ (with $d<<p$), i.e. the datapoints
$\boldsymbol{x}_{i}$ lie on or near a $d$-dimensional manifold in the
$\mathbb{R}^{p}$ space.

It is often easier to work if we restrict the attention of the Riemannian
manifold only, but the general problem makes no assumption about the manifold
and allows discontinuous clusters to be present. The goal is to obtain a low
dimensional representation $\boldsymbol{y}_{i}\in\mathbb{R}^{d}$ corresponding
to each $\boldsymbol{x}_{i}$ such that the geometry of the manifold is
preserved in the low dimensional space as much as possible. Since the problem
does not associate a metric to measure the distortion, it is an ill-posed
problem. Almost all of the existing DR algorithms define objective functions
(or a metric) which is a characterization of the dissimilarity between
geometric properties of the high and low dimensional manifolds, and then it
tries to minimize that metric to obtain the “best” lower dimensional
representation. Since there is no universal characterization, most of these
objective function building process is ad-hoc in nature. Later, in Section
2.1, we shall provide a unified framework to define the problem more
concretely, and also formally define what is meant by a dimensionality
reduction algorithm.

###  1.2 Existing Dimensionality Reduction Algorithms

####  1.2.1 Metric Preserving Algorithms

Principal Component Analysis (PCA) is the most popular metric-preserving DR
algorithm. It aims to find the linear projections
$\boldsymbol{\mathsf{y}}=\boldsymbol{M}\boldsymbol{\mathsf{x}}$ of the data
$\boldsymbol{\mathsf{x}}$ to preserve the highest variability. The objective
becomes same as maximizing
$\text{trace}\left(\boldsymbol{M}^{\intercal}\text{Cov}(\boldsymbol{\mathsf{x}})\boldsymbol{M}\right)$.
An analytical solution of this problem can be obtained by finding the
normalized eigenvectors associated with the $d$ largest eigenvalues of the
empirical variance-covariance matrix $\boldsymbol{\Sigma}$ of the zero mean
data (i.e. the centered data). Although it has been used in various
applications [48, 37, 20, 23], it suffers from many problems.

  1. 1.

PCA restricts its attention to only linear functions of the variables.
However, if the underlying manifold is nonlinear, PCA cannot recover such
structure. Localized PCA [24, 38] and Kernel PCA [42, 41] are two extensions
that considers functions in an RKHS (see Appendix A.1) instead of only linear
functions to counter this problem.

  2. 2.

It is not necessarily true that the information contained in the dataset must
be explained by variance (not by other moments). An instance of this problem
is described in detail in Example 1.

######  Example 1.

Since PCA, Kernel PCA or probabilistic PCA uses the variance in the data as a
measure of information content, it can fail when most of the information
presented in the data is along the direction with minimum variance. For
example, if the dataset contains two long linear clusters as shown in Figure
1, PCA fails to detect the shape or the structure of the manifold. A desirable
dimensionality reduction algorithm should be able to detect these clusters,
and create a one-dimensional mapping where the clusters are mapped to some
disjoint intervals in $\mathbb{R}$.

(a) A 2-dimensional dataset with two linear clusters

(b) Output of 1 dimensional PCA with random noise in $y$-axis

Figure 1: PCA fails to capture the class difference by first principal
component

  3. 3.

The third problem is the trustability issues. Since PCA takes the approach to
maximize the variances, and the solution is given by the principal
eigenvectors of the covariance matrix. These principal eigenvectors have an
asymptotic normal distribution irrespective of the distribution of the
original $\boldsymbol{\mathsf{x}}$ [30], hence the top few principal
components tend to look like Gaussian distributed random variables. A
simulated experiment has been performed in Example 2 to visualize this issue.
Although all the principal components taken together returns only a rotated
version of the original dataset, which makes PCA very desirable, its lower
dimensional embeddings tend to distort the distribution of the data, which is
undesirable. Any statistical inference based on the PCA reduced form could
then be misleading in nature.

######  Example 2.

As noted earlier, PCA has trustability issues as the principal eigenvectors
are asymptotically normally distributed [30], hence the low dimensional
representation usually tend to look like normally distributed. A simulation
was performed to demonstrate this effect.

  1. 1.

In first simulation, $1000$ datapoints are generated from
$N(\boldsymbol{0},\boldsymbol{I}_{5})$. From the data, the first two principal
component directions were estimated, and then the projection onto them is
obtained. As shown in Figure 2, the algorithm outputs a 2-dimensional Gaussian
data as best embedding of the original $5$-dimensional data.

  2. 2.

In the second simulation, $1000$ datapoints are generated from the unit
$5$-dimensional hypercube, where each of the variables are i.i.d. uniform
$(0,1)$. Again, the first two principal components are extracted and plotted
in Figure 3. The output again looks similar to a Gaussian distribution.

Figure 2: First two principal components as extracted by normally distributed
noise data

Figure 3: First two principal components as extracted by uniformly distributed
noise data

Therefore, despite the original data being very different, the reduction
outputted by PCA is similar, resulting in similar inference about the data
(whenever the analysis first performs PCA as dimensionality reduction and
analyzes the reduced data).

####  1.2.2 Graphical Algorithms

Another approach to dimensionality reduction is by approximating the geodesic
distance through the curved manifold. These algorithms start by creating a
neighbourhood graph $G$ where each datapoint $\boldsymbol{x}_{i}$ is connected
to its $k$-nearest neighbours. The length of the shortest path between two
points $\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$ through the graph $G$
serves as an approximation of the geodesic distance between
$\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$. Isomap, a popular graph-based
DR algorithm [29, 35] uses multidimensional scaling to find low dimensional
points to approximate these distances through $G$. Uniform Manifold
Approximation and Projection (UMAP) [31], one of the most recent DR
algorithms, creates the $k$-neighbourhood graph $G$ and computes the
simplicial sets [13] generated by that neighbourhood graph. Then the high and
low dimensional simplicial sets are matched against each other.

These graph-based algorithms mainly suffer from two problems. The first one
being the choice of the parameter $k$, the optimal number of neighbours, which
if not properly chosen, may lead to artificial holes. For the purpose of
demonstration, we generated $100$ datapoints uniformly from the unit square
and obtained the neighbourhood graphs for $k=1,3$ and $10$. As shown in Figure
4, for $k=1$, the neighbourhood graph is highly disconnected, and even for
$k=10$, there remains an artificial hole in the dataset. There is no consensus
about the choice of $k$ in a data-dependant way, since there is no established
theory about what metric should one look at.

Figure 4: Neighbourhood graphs with $k=1,3,10$, for $100$ points generated
uniformly in the unit square

On the other hand, another problem occurs when there are tightly clustered
points having small intra-cluster distances and having high inter-cluster
distances. In such a case, the neighbourhood graph often becomes disconnected.
Therefore, the preservation of distances based on the neighbourhood graph
leads to distortion of the global structure of the dataset. UMAP, the most
recent graph-based DR algorithm, suffers from a similar problem [50, 49, 36].
The algorithm assumes that the data is uniformly distributed on the manifold
$\mathcal{M}$, and hence, away from any boundaries, any ball of fixed volume
should contain approximately the same number of points of $\boldsymbol{X}$
regardless of where on the manifold it is centered. However, if there are more
than one clusters present in the data, the neighbourhood graph $G$ becomes
disconnected, i.e. for two clusters present in the data, the neighbourhood
graph can be split into two disjoint neighbourhood graph $G=G_{1}\cup G_{2}$,
where $G_{1},G_{2}$ are neighbourhood graphs for each of the clusters. In that
case, the ordering or the position of the clusters have very little effect,
since the disconnected graphs can be rearranged without any constraint. For
this reason, UMAP fails to convey the ordering of the clusters (present in the
original data), into the reduced dimensionality, as shown in Figure 5.

(a) The Neighborhood graph

(b) 2 dimensional projection by UMAP

Figure 5: UMAP fails to preserve global structure because of disconnected
graphs

A few other graph-based algorithms such as Locally Linear Embedding (LLE)
[39], Laplacian eigenmaps [1] etc., all suffer from similar problems described
above.

####  1.2.3 Supervised Dimensional Reduction Algorithms

Sufficient Dimensionality Reduction Framework describes the supervised
framework for dimensionality reduction. Any kind of dimensionality reduction
analysis is generally followed by clustering, classification, regression
analysis using the reduced set of variables. If such classification labels or
regression variables are known a-priori, one may use this information on
response variable in formulation of the dimensionality reduction problem, so
that the reduced set of variables is more meaningful. It starts with a
framework given as
$\boldsymbol{\mathsf{Y}}\perp\\!\\!\\!\perp\boldsymbol{\mathsf{X}}\mid\sigma\left(\left\\{\phi_{1}(\boldsymbol{\mathsf{X}}),\dots\phi_{d}(X)\right\\}\right)$,
where $\phi_{1},\dots\phi_{d}$ are unknown functions to be estimated.
Specifically for linear embedding functions, one needs to estimate
$\beta_{1},\dots\beta_{d}$ such that
$\boldsymbol{\mathsf{Y}}\perp\\!\\!\\!\perp\boldsymbol{\mathsf{X}}\mid(\beta_{1}^{\intercal}\boldsymbol{\mathsf{X}},\beta_{2}^{\intercal}\boldsymbol{\mathsf{X}},\dots\beta_{d}^{\intercal}\boldsymbol{\mathsf{X}})$.
Bing Li [27] provides a useful discussion on the numerous methods that has
been developed to solve this problem. Popular methods under sufficient
dimensionality reduction framework includes,

  1. 1.

Sliced Inverse Regression [28, 53]. Here, the idea is that the difference
between the conditional expectations
$\mathbb{E}(\boldsymbol{\mathsf{X}}\mid\boldsymbol{\mathsf{Y}}=y_{1})-\mathbb{E}(\boldsymbol{\mathsf{X}}\mid\boldsymbol{\mathsf{Y}}=y_{2})$
contains only the effect due the principal SDR directions.

  2. 2.

Sliced Average Variance Estimation (SAVE) [11]. The basic idea is that
$\text{Var}(\boldsymbol{\mathsf{X}}\mid\boldsymbol{\mathsf{Y}}=y)$ contains
only the variation due to the directions orthogonal to SDR.

  3. 3.

Generalized SIR and Generalized SAVE attempts to solve the general problem by
restricting the $\sigma$-field
$\sigma\left(\left\\{\phi_{1}(\boldsymbol{\mathsf{X}}),\dots\phi_{d}(X)\right\\}\right)$
to a RKHS instead.

The main problem with the SDR approach is that it performs poorly if the
response variable is categorical. For example, if there are $c$ levels in the
categorical variable $\boldsymbol{\mathsf{Y}}$, then only $(c-1)$ SDR
directions can be obtained. Also, since the reduction is performed only with a
particular response in mind, the reduced set of variables cannot be reused in
other problems. An example of how SIR may result in loss of information has
been given as follows.

######  Example 3.

We apply dimensionality reduction for Wisconsin breast cancer data based on
two algorithms, tSNE and SIR. As shown in Figure 6, they output a very
different shape of the reduced data. While SIR only tries to best preserves
the classification of benign and malignant cells, tSNE does it unsupervisedly,
and outputs a one dimensional curvilinear manifold of the data. It was seen
that the horizontal axis outputted by tSNE roughly aligns to the size of the
tumour, with a strongly negative correlation of $-0.7013$. The tSNE output
also says that it is much easier to classify the malignant and benign cells if
the tumour is large, while it is much harder to segregate if it is small. SDR
techniques will generally fail to consider this as it is optimized only to
discriminate between the two classes.

(a) Dimensionality reduction via tSNE

(b) Dimensionality reduction via SIR

Figure 6: Dimensionality reduction via tSNE and SIR for Winconsin breast
cancer data, Red means malignant and Turquoise means benign cancer cells.

###  1.3 Existing Metrics

As described in Section 1.2, despite the existence of a large number of
dimensionality reduction algorithms, most of these suffer from some kind of
trustability issues. Also, there is little consensus about how to evaluate the
performance of a dimensionality reduction algorithm since there exists a
diverse range of frameworks under which the existing algorithms are derived
and they use different metrics for defining their corresponding loss
functions. Let us define $A_{k}(\boldsymbol{x}_{i})$ (and
$B_{k}(\boldsymbol{y}_{i})$) as the set of indices $j$ such that the $j$-th
datapoint is a $k$ nearest neighbour of the $i$-th datapoint in high (and low)
and $r_{k,i}(\boldsymbol{x}_{j})$ (and
$\widetilde{r}_{k,i}(\boldsymbol{y}_{j})$) denote the corresponding ranks
according to the distances from $i$-th datapoint in high (and low) dimensional
space, with the nearest neighbour having rank $1$. Almost all of the existing
metrics to evaluate dimensionality reduction work on the idea of preserving
local structures, as follows

  1. 1.

Separability indices [18, 34] such as
$\text{TSI}(\boldsymbol{X},\boldsymbol{Y})=1-\dfrac{1}{n}\sum_{i=1}^{n}\dfrac{1}{k}\sum_{j\in
A_{k}(\boldsymbol{x}_{i})}\boldsymbol{1}_{j\not\in
B_{k}(\boldsymbol{y}_{i})},$

  2. 2.

“Trustworthiness” and “continuity” measures proposed by Kaski et al. [25],
given by

Trustworthiness
$\displaystyle=1-\dfrac{2}{nk(2n-3k-1)}\sum_{i=1}^{n}\sum_{j\in
B_{k}(\boldsymbol{y}_{i})\backslash
A_{k}(\boldsymbol{x}_{i})}(r_{k,i}(\boldsymbol{y}_{j})-k)$ Continuity
$\displaystyle=1-\dfrac{2}{nk(2n-3k-1)}\sum_{i=1}^{n}\sum_{j\in
A_{k}(\boldsymbol{x}_{i})\backslash
B_{k}(\boldsymbol{y}_{i})}(\widetilde{r}_{k,i}(\boldsymbol{x}_{j})-k)$

These two quantities are created based on the intuitive idea that to measure
trustworthiness of the output of an algorithm, if the low dimensional
visualization shows that two points are close together, then they should also
remain close together in the original dataset. The continuity is the other
side of the coin, points in a vicinity in the original space should be
projected close together, so that the visualization of a new test datapoint
can be obtained based on the proximity calculations.

###  1.4 Contributions of this study

Despite having a large number of different metrics for evaluating different
dimensionality reduction algorithms, there is no consensus about a single
metric. Also, most of these metrics, using some form of $k$-nearest
neighbours, only focus on retaining the local structure, then they are only
meaningful when the reduced data is used for clustering and classification
using nearest neighbour-based algorithms. Our contribution would be to
formally define a dimensionality reduction framework, under which we will
provide a precise definition for an algorithm to qualify as a dimensionality
reduction algorithm such that almost all of the existing dimensionality
reduction algorithms will obey the definition. Then, we will consider two
metrics under that general setup based on two important properties of the
ideal dimensionality reduction algorithm, named as “Trustability index” and
“Consistency index”.

While these two metrics will allow us to compare performances of different
algorithms, we shall try to explore how these metrics can be reduced and what
are the essential components of an algorithm that naturally correlates with
these metrics. These, we shall try to build a dimensionality reduction
algorithm by integrating those components which should guarantee better
performances than the existing algorithms.

For brevity of presentation, the proofs of all the theorems are given in the
Appendix B.

##  2 A Generalized Framework

###  2.1 The Framework for Dimensionality Reduction

We start with a generic setup where $\boldsymbol{\mathsf{x}}$ is a
$\mathbb{R}^{p}$-valued random variable where $p$ is large. The support of the
distribution of $\boldsymbol{\mathsf{x}}$ is approximately a lower dimensional
manifold. Denoting this intrinsic dimension as $d$, there should exist
$\phi_{1},\dots\phi_{d}$ functions from $\mathbb{R}^{p}\rightarrow\mathbb{R}$
such that each of them represents a coordinate system along the manifold. When
one takes account of all these coordinates, the part left out must be a random
variation that the coordinate system cannot take an account for. With this
interpretation, we define the following.

######  Definition 1.

A $d$-dimensional reduction for the $\mathbb{R}^{p}$ valued random variable
$\boldsymbol{\mathsf{x}}$ is a specification of the functions
$\phi_{1},\phi_{2},\dots\phi_{d}$ with each
$\phi_{i}:\mathbb{R}^{p}\rightarrow\mathbb{R}$ for $i=1,2,\dots d$ and
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ such that,

  1. 1.

$\boldsymbol{\mathsf{x}}=f\left(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}})\right)+\epsilon(\boldsymbol{\mathsf{x}})$
almost surely.

  2. 2.

$\phi_{1},\phi_{2},\dots\phi_{d}$ are measurable (may be continuous and
smooth) functions. Each of these embeds the original data into the lower
dimensional space, hence are called embedding functions.

  3. 3.

$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ is a measurable reconstruction
function that approximately retrieves $\boldsymbol{\mathsf{x}}$.

  4. 4.

$\epsilon(\boldsymbol{\mathsf{x}})$ is the error in reconstruction, which is
independent of the $\sigma$-field generated by
$\\{\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}})\\}$.

  5. 5.

$\mathbb{E}(\epsilon(\boldsymbol{\mathsf{x}}))=0$.

Following are some of the examples to show that the above definition is a very
general one.

######  Example 4.

The first and trivial example of a dimensional reduction is the set of
coordinate functions (or any permutation of that). Formally, a $p$-dimensional
reduction for any $\mathbb{R}^{p}$-valued random variable
$\boldsymbol{\mathsf{x}}$ is given by $\phi_{i}(\boldsymbol{x})=x_{\pi(i)}$,
where $x_{i}$ is the $i$-th coordinate of $\boldsymbol{x}$ and
$\pi:\\{1,2,\dots p\\}\rightarrow\\{1,2,\dots p\\}$ is a permutation function.
Clearly, $\epsilon(\boldsymbol{\mathsf{x}})=0$ and
$f(x)=(x_{\pi^{-1}(1)},\dots x_{\pi^{-1}(p)})$.

######  Example 5.

Let, $\boldsymbol{\mathsf{x}}$ has a multivariate Gaussian distribution. Then,
starting with any orthonormal basis
$B=\\{\boldsymbol{v}_{1},\dots,\boldsymbol{v}_{d},\dots\boldsymbol{v}_{p}\\}$
(for $d<p$), one may decompose,
$\boldsymbol{\mathsf{x}}=\sum_{i=1}^{d}\phi_{i}(\boldsymbol{\mathsf{x}})\boldsymbol{v}_{i}+\sum_{i=(d+1)}^{p}\psi_{i}(\boldsymbol{\mathsf{x}})\boldsymbol{v}_{i}$.
Here,
$\phi_{i}(\boldsymbol{\mathsf{x}})=\boldsymbol{v}_{i}^{\intercal}\boldsymbol{\mathsf{x}}$
are the embedding functions and $f(x_{1},\dots
x_{d})=\sum_{i=1}^{d}x_{i}\boldsymbol{v}_{i}$ is the reconstruction function.
Due to multivariate normality of $\boldsymbol{\mathsf{x}}$ and orthogonality
of $B$, the error
$\sum_{i=(d+1)}^{p}\psi_{i}(\boldsymbol{\mathsf{x}})\boldsymbol{v}_{i}$ is
independent of the $\sigma$-field generated by the embedding functions.

Note that, the same example also holds if $\boldsymbol{\mathsf{x}}$ follows a
mixture of Gaussian distributions. Any linear dimensionality reduction
technique that computes the embedding functions with respect to some
orthogonal basis (including PCA) falls under this category.

######  Theorem 1.

If there exists two measurable functions $f$ and $g$ such that

$\boldsymbol{\mathsf{x}}=f(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))+\epsilon_{1}(\boldsymbol{\mathsf{x}})=g(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))+\epsilon_{2}(\boldsymbol{\mathsf{x}})$

i.e. both can be regarded as reconstruction function corresponding to the
reduction $(\phi_{1},\dots\phi_{d})$, then both reconstructions are same
almost surely.

The implication of this result is that to determine such a dimensional
reduction, it is enough to specify the embedding functions
$\phi_{1},\phi_{2},\dots\phi_{d}$ only, since in view of Theorem 1, any
reconstruction function will output the same decoded version almost surely.

###  2.2 Dimensionality reduction algorithms

Now, turning to the real world datasets, we do not have the random variable
$\boldsymbol{\mathsf{x}}$, but instead we have i.i.d. data
$\boldsymbol{\mathsf{x}}_{1},\dots\boldsymbol{\mathsf{x}}_{n}$, each being
$p$-dimensional. Therefore, any dimensionality reduction algorithm basically
defines a loss function (convex or nonconvex) [51] which is a function of the
original datapoints, the embedding outputs
$\mathsf{y}_{ij}=\phi_{j}(\boldsymbol{\mathsf{x}}_{i})$ for $j=1,2,\dots d$
and $i=1,2,\dots n$, the reconstruction function $f$, and the errors
$\epsilon$. However, in view of Theorem 1, the specification of $f$ is not
necessary when $\mathsf{y}_{ij}$ is known, and also when the reconstruction
and the errors are known, specification of the original dataset is not needed.
Therefore, one may consider a loss function as
$L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i}):i=1,\dots n;j=1,\dots
d\\},\epsilon(\boldsymbol{\mathsf{x}}_{1}),\dots\epsilon(\boldsymbol{\mathsf{x}}_{n})\right)$
for $n$ observations. However, such specification requires different loss
functions to be specified for each choice of the intrinsic dimension $d$
(which is not known a-priori), as the loss function above is a function from
$\mathbb{R}^{nd}\times\mathbb{R}^{np}$ to $\mathbb{R}$. This requires special
analysis for each of these loss functions, hence we extend this notion and
consider only a single loss function $L:\mathbb{R}^{2np}\rightarrow\mathbb{R}$
given by

$L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i}):i=1,\dots n;j=1,\dots
p\\},\epsilon(\boldsymbol{\mathsf{x}}_{1}),\dots\epsilon(\boldsymbol{\mathsf{x}}_{n})\right)$

while we perform the minimization of this loss subject to the constraints that
$\phi_{j}(\boldsymbol{\mathsf{x}})=0$ for all $j=(d+1),\dots p$ and for all
$\boldsymbol{\mathsf{x}}$, in order to estimate the embedding functions.

######  Definition 2.

A dimensionality reduction algorithm $\mathcal{A}_{n}$ is a function from
$\\{1,2,\dots p\\}\times\mathbb{R}^{np}$ to a restricted class $\mathcal{C}$
of functions (say smooth, linear, piecewise linear etc.) along with a
specification of the loss function $L:\mathbb{R}^{2np}\rightarrow\mathbb{R}$.
This means, the algorithm takes the first argument as the number of embedding
functions (i.e the lower dimension $d$) and the second argument as the $n$
datapoints each of dimension $p$. The algorithm then chooses the best possible
reduction by minimizing the associated loss function for any $d=1,2,\dots p$,

$(\widehat{\phi}_{1},\dots\widehat{\phi}_{d})=\arg\min_{\phi_{1},\dots\phi_{d}\in\mathcal{C}}L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i}):i=1,\dots
n;j=1,\dots
p\\},\epsilon(\boldsymbol{\mathsf{x}}_{1}),\dots\epsilon(\boldsymbol{\mathsf{x}}_{n})\right),\\\
\text{subject to the constraints }\phi_{j}(\boldsymbol{\mathsf{x}})=0\ \forall
j=(d+1),\dots p,\ \forall\boldsymbol{\mathsf{x}}\in\mathbb{R}^{p}$ (1)

The embedded lower dimensional representation is denoted by
$\boldsymbol{\mathsf{Y}}=\mathcal{A}_{n}(d,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})$,
where $\boldsymbol{\mathsf{Y}}$ is a $n\times d$ matrix, whose rows are
$\boldsymbol{\mathsf{y}}_{i}=(\widehat{\phi}_{1}(\boldsymbol{\mathsf{x}}_{i}),\dots\widehat{\phi}_{d}(\boldsymbol{\mathsf{x}}_{i}))$.

The loss function should not be any arbitrary one, there are some obvious
restrictions.

  1. 1.

$L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)=L\left(\\{a_{j}+\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)$,
i.e. a translation to the embedding functions preserves the shape and the
structure of the data, hence the loss function should not change. Note that,
the translation amount may be different in different coordinates.

  2. 2.

$L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)=L\left(\\{(T\circ\phi)_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)$,
where $T:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ is an orthonormal
transformation. Since, orthonormal transformations simply rotate the original
data, the structure and shape of the data is preserved.

  3. 3.

$L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)=L\left(\\{a\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)$.
If all the axis are scaled by the same amount, the loss function should not
change. If the scaling factors are different, the shape of the function may
distort.

  4. 4.

For any other transformation that distorts the shape of the data, the loss
function should output something different.

  5. 5.

$L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)\geq
0$ and is equal to $0$ if and only if
$\epsilon(\boldsymbol{\mathsf{x}}_{i})=0$ for all $i=1,\dots n$.

Note that, the rotational invariance imply permutation invariance, hence the
optimization problem given in Eq. (1) can also be re-framed by constraining
$\phi_{j}(\boldsymbol{\mathsf{x}})=0$ for any $(p-d)$ functions, while the
optimization is done on the rest.

###  2.3 Nonparametric estimator of out-of-the-sample reduction

Many dimensionality reduction algorithms (which specifically uses non-convex
loss functions), only outputs the embedded lower dimensional representations
$\boldsymbol{\mathsf{y}}_{1},\dots\boldsymbol{\mathsf{y}}_{n}$ (each being $d$
dimensional) of the training dataset
$\boldsymbol{\mathsf{x}}_{1},\dots\boldsymbol{\mathsf{x}}_{n}$ (each being $p$
dimensional) without explicitly specifying the embedding functions
$\phi_{1},\dots\phi_{d}$, in accordance to Defn. 2. However, if one wants to
find the embedding of a test datapoint $\boldsymbol{\mathsf{x}}$, one needs to
know either a nonparametric or parametric form of $\phi_{j}$’s. The following
theorem precisely provides that under the general setup considered so far
assuming the search space $\mathcal{C}$ for an algorithm $\mathcal{A}_{n}$ is
subset of some Reproducible Kernel Hilbert Space (RKHS) (see Appendix A.1).

######  Theorem 2.

For any dimensionality reduction algorithm $\mathcal{A}_{n}$, given that the
restricted class of functions in the search space $\mathcal{C}$ is a subset of
some RKHS for some kernel $k$, then the embedding functions $\phi_{j}$ for any
$j=1,2,\dots d$ can be evaluated as,

$\phi_{j}(\boldsymbol{\mathsf{x}})=\boldsymbol{\mathsf{y}}_{j}^{\intercal}\boldsymbol{K}^{-1}\boldsymbol{K}(\boldsymbol{\mathsf{x}})$
(2)

where,

$\boldsymbol{\mathsf{y}}_{j}=\begin{pmatrix}\mathsf{y}_{1j}\\\ \dots\\\
\mathsf{y}_{nj}\end{pmatrix},\
\boldsymbol{K}=\begin{pmatrix}k(\boldsymbol{\mathsf{x}}_{1},\boldsymbol{\mathsf{x}}_{1})&k(\boldsymbol{\mathsf{x}}_{1},\boldsymbol{\mathsf{x}}_{2})&\dots&k(\boldsymbol{\mathsf{x}}_{1},\boldsymbol{\mathsf{x}}_{n})\\\
k(\boldsymbol{\mathsf{x}}_{2},\boldsymbol{\mathsf{x}}_{1})&k(\boldsymbol{\mathsf{x}}_{2},\boldsymbol{\mathsf{x}}_{2})&\dots&k(\boldsymbol{\mathsf{x}}_{2},\boldsymbol{\mathsf{x}}_{n})\\\
\vdots&\vdots&\ddots&\vdots\\\
k(\boldsymbol{\mathsf{x}}_{n},\boldsymbol{\mathsf{x}}_{1})&k(\boldsymbol{\mathsf{x}}_{n},\boldsymbol{\mathsf{x}}_{2})&\dots&k(\boldsymbol{\mathsf{x}}_{n},\boldsymbol{\mathsf{x}}_{n})\end{pmatrix},\
\text{ and
}\boldsymbol{K}(\boldsymbol{\mathsf{x}})=\begin{pmatrix}k(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{x}}_{1})\\\
\dots\\\ k(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{x}}_{n})\end{pmatrix},$

and $\mathsf{y}_{ij}=\phi_{j}(\boldsymbol{\mathsf{x}}_{i})$ for any
$i=1,2,\dots n;j=1,2,\dots d$. This gives us a way to obtain nonparametric
estimator of the low dimensional embedding on the basis of observed embedding
for training samples, for any dimensionality reduction algorithm.

###  2.4 Building the reconstruction function

This section concerns with the derivation of the reconstruction function $f$.
Theorem 1 says that given the specific embedding functions, the reconstruction
function is unique almost surely, it does not provide a constructive way to
produce it. This is difficult to obtain in its generality, because the
independence of $\epsilon$ and
$\sigma\left(\\{\phi_{1},\dots\phi_{d}\\}\right)$ is hard to verify.
Therefore, we shall assume that each coordinate function $f_{l}$ of the
reconstruction function $f$ belongs to a RKHS. This simplification actually
allows one to reduce the independence condition to zero-covariance operator in
RKHS domain (see Theorem 10 for details).

######  Theorem 3.

For any dimensionality reduction algorithm $\mathcal{A}_{n}$ with embedding
functions
$\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}})$ and
the training dataset
$(\boldsymbol{\mathsf{x}}_{1},\dots\boldsymbol{\mathsf{x}}_{n})$, if

  1. 1.

Each $\phi_{j}(\cdot)$ is a member of some RKHS $\mathcal{H}_{x}$ generated by
a universal kernel function $k_{x}$ given by
$\phi_{j}(\boldsymbol{\mathsf{x}})=\sum_{i=1}^{n}\alpha_{ij}k_{x}(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{x}}_{i})$
with
$\mathbb{E}(k_{x}(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{x}}))<\infty$.

  2. 2.

Each coordinate function of $f$ i.e. $f_{l}$ is a member of some RKHS
$\mathcal{H}_{y}$ generated by a universal kernel function $k_{y}$ with
$\mathbb{E}(k_{y}(\boldsymbol{\mathsf{y}},\boldsymbol{\mathsf{y}}))<\infty$.

  3. 3.

The probability measure $P_{x,y}$ defining the joint distribution of
$(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{y}})$ has the support as a
compact subset of $\mathbb{R}^{p}\times\mathbb{R}^{d}$.

Then the coordinate functions of the reconstruction function $f$ can be
expressed as,

$f_{l}(\boldsymbol{\mathsf{y}})=n^{-1}\sum_{i=1}^{n}\boldsymbol{\mathsf{x}}_{il}+\sum_{i=1}^{n}\beta_{il}\left(k_{y}(\boldsymbol{\mathsf{y}},\boldsymbol{\mathsf{y}}_{i})-n^{-1}\sum_{j=1}^{n}k_{y}(\boldsymbol{\mathsf{y}}_{j},\boldsymbol{\mathsf{y}}_{i})\right),\
l=1,2,\dots p$

where
$\boldsymbol{\mathsf{y}}_{i}=(\mathsf{y}_{i1},\dots\mathsf{y}_{id})^{\intercal}$
and $\mathsf{y}_{ij}=\phi_{j}(\boldsymbol{\mathsf{x}}_{i})$ and
$\beta_{il}=\boldsymbol{A}\left(\boldsymbol{A}^{\intercal}\boldsymbol{A}\right)^{-1}\boldsymbol{c}_{l}$,
and the matrix $\boldsymbol{A}_{n\times d}$ has entries
$a_{rj}=(\boldsymbol{G}_{y}\boldsymbol{G}_{x}\boldsymbol{\alpha})_{rj}$,

$\boldsymbol{\alpha}=\begin{bmatrix}\alpha_{11}&\dots&\alpha_{1d}\\\
\vdots&\ddots&\vdots\\\ \alpha_{n1}&\dots&\alpha_{nd}\end{bmatrix}_{n\times
d}$

and $\boldsymbol{G}_{x},\boldsymbol{G}_{y}$ are gram matrices corresponding to
RKHS $\mathcal{H}_{x}$ and $\mathcal{H}_{y}$ as given in Eq. (3)-Eq. (4), and
$\boldsymbol{c}_{l}=(c_{1l},\dots c_{dl})$ where

$c_{jl}=n^{-1}\sum_{i=1}^{n}\mathsf{x}_{il}\mathsf{y}_{ij}-\left(n^{-1}\sum_{i=1}^{n}\mathsf{x}_{il}\right)\left(n^{-1}\sum_{i=1}^{n}\mathsf{y}_{ij}\right)$

###  2.5 Trustability and Consistency index

A dimensionality reduction algorithm $\mathcal{A}_{n}$ is characterized by the
loss function $L(\phi_{1},\dots\phi_{p},\epsilon)$ according to Definition 2,
hence one may measure the trustability of the algorithm based on the loss
function itself.

Consider the following situation: starting with any algorithm
$\mathcal{A}_{n}$ and a $p$-dimensional dataset
$\boldsymbol{\mathsf{X}}_{n\times p}$, we ask the algorithm to reduce the data
into a $p$-dimensional output $\boldsymbol{\mathsf{Y}}_{n\times p}$. Clearly,
the original data itself becomes the “best” possible reduction achieving $0$
loss as $\epsilon(\boldsymbol{\mathsf{x}}_{i})=0$ for all $i=1,\dots n$. Since
the algorithm minimizes the loss function given in Eq. (1), which is invariant
under translation, uniform scaling and rotation, we should have
$\boldsymbol{\mathsf{Y}}\approx\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}+\lambda\boldsymbol{P}\boldsymbol{\mathsf{X}}$
for some mean vector $\boldsymbol{\mu}$, real number $\lambda$ and orthogonal
matrix $\boldsymbol{P}$. Therefore, the quantity
$\min_{\boldsymbol{\mu},\lambda,\boldsymbol{P}}\|\boldsymbol{\mathsf{Y}}-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{P}\boldsymbol{\mathsf{X}}\|_{2}^{2}$
could be regarded a measure of how trustworthy the algorithm is in preserving
shape and structure of the original dataset. The above minimization problem
can be exactly solved by the Procrustes problem (see Appendix A.3), which
results in the construction of a Trustability Index for any dimensionality
reduction algorithm.

######  Definition 3.

The Trustability Index of any dimensionality reduction algorithm
$\mathcal{A}_{n}$ for any dataset $\boldsymbol{\mathsf{X}}_{n\times p}$ is
given by,

$\text{TI}(\mathcal{A}_{n},\boldsymbol{\mathsf{X}}):=\sum\text{sing}(\boldsymbol{\Sigma}_{\boldsymbol{\mathsf{YY}}})-\left(\sum\text{sing}(\boldsymbol{\Sigma}_{\boldsymbol{\mathsf{XX}}})\right)^{-1}\left(\sum\text{sing}\boldsymbol{\Sigma}_{\boldsymbol{\mathsf{XY}}}\right)^{2}$

where $\sum\text{sing}(\boldsymbol{A})$ denotes the sum of the singular values
of the matrix $\boldsymbol{A}$, and

$\boldsymbol{\Sigma}_{\boldsymbol{\mathsf{XX}}}=(\boldsymbol{\mathsf{X}}-\overline{\boldsymbol{\mathsf{X}}})^{\intercal}(\boldsymbol{\mathsf{X}}-\overline{\boldsymbol{\mathsf{X}}}),\
\boldsymbol{\Sigma}_{\boldsymbol{\mathsf{XY}}}=(\boldsymbol{\mathsf{X}}-\overline{\boldsymbol{\mathsf{X}}})^{\intercal}(\boldsymbol{\mathsf{Y}}-\overline{\boldsymbol{\mathsf{Y}}}),\
\boldsymbol{\Sigma}_{\boldsymbol{\mathsf{YY}}}=(\boldsymbol{\mathsf{Y}}-\overline{\boldsymbol{\mathsf{Y}}})^{\intercal}(\boldsymbol{\mathsf{Y}}-\overline{\boldsymbol{\mathsf{Y}}})$

where
$\boldsymbol{\mathsf{Y}}=\mathcal{A}_{n}(p,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})$
and $\overline{\boldsymbol{\mathsf{X}}},\overline{\boldsymbol{\mathsf{Y}}}$
are the $n\times p$ matrices whose each row is
$n^{-1}\boldsymbol{1}_{n}^{\intercal}\boldsymbol{\mathsf{X}}$ and
$n^{-1}\boldsymbol{1}_{n}^{\intercal}\boldsymbol{\mathsf{Y}}$ respectively.

######  Remark 1.

The more trustable an algorithm is, the less is the value of the index. This
means, the most trustworthy algorithm attains the index value equal to $0$.

######  Example 6.

The algorithm that has a $0$-$1$ loss function, outputs a random subset of
size $d$ of the available $p$ variables. Thus, with $d=p$, it outputs a
permutation of columns of $\boldsymbol{\mathsf{X}}$, hence attains TI equal to
0. Similarly, for Principal Component Analysis, since
$\boldsymbol{\mathsf{Y}}$ is simply a orthogonal transformation of
$\boldsymbol{\mathsf{X}}$, it is also most trustable algorithm.

Naturally, the above algorithms are not reasonably well if $d<p$. Therefore,
along with Trustability Index, we need another index which is competing in
nature. The problem in the aforementioned algorithms is that they have a large
variability in its output under different transformations of the original
data.

######  Theorem 4.

Let, $\boldsymbol{z}=T(\boldsymbol{x})$ be some invertible transformation
$T:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$. For a dimensionality reduction
algorithm $\mathcal{A}_{n}$, with loss function $L$ and restricted class of
functions $\mathcal{C}$, if

  1. 1.

For any $\phi\in\mathcal{C}$, both $\phi\circ T$ and $\phi\circ T^{-1}$ belong
to $\mathcal{C}$.

  2. 2.

$\epsilon(\boldsymbol{\mathsf{x}}_{i})=\epsilon(\boldsymbol{\mathsf{z}}_{i})=\boldsymbol{\mathsf{\epsilon}}_{i}$
for all samples $\boldsymbol{\mathsf{x}}_{1},\dots\boldsymbol{\mathsf{x}}_{n}$
and transformed samples
$\boldsymbol{\mathsf{z}}_{i}=T(\boldsymbol{\mathsf{x}}_{i})$.

Then,

$\min_{\phi_{1},\dots\phi_{d}\in\mathcal{C}}L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right)=\min_{\psi_{1},\dots\psi_{d}\in\mathcal{C}}L\left(\\{\psi_{j}(\boldsymbol{\mathsf{z}}_{i})\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right)$

where both the above minimization are performed subject to the conditions that
$\phi_{j}(\boldsymbol{\mathsf{x}})=\psi_{j}(\boldsymbol{\mathsf{z}})=0$ for
all $j=(d+1),\dots p$.

Theorem 4 opens up the possibility that even if the underlying intrinsic
manifold is changed, as long as the errors in reconstruction remain the same,
the loss also remains the same. For example, in Figure 7, we demonstrate two
different 2D datasets, both of which have an underlying one-dimensional
manifold. Clearly, the underlying one-dimensional manifolds can be transformed
into each other employing a simple one-to-one transformation, hence a
desirable algorithm should output the same one-dimensional coordinate
representation for both these datasets. In Theorem 4, the assumption that the
errors remain the same even after the transformation is crucial. Therefore,
even if $T(\boldsymbol{x})=\lambda\boldsymbol{x}$ for any $\lambda>0$ is a
one-to-one transformation, it cannot be readily applied to the data as it may
shrink (or inflate) both the errors and the extent of the low dimensional
manifold simultaneously. We, on the other hand, want to apply this
transformation only to the underlying manifold and keeping the reconstruction
error fixed. Therefore, we may employ the following scheme.

Figure 7: Two high dimensional data which have equivalent low dimensional
underlying manifold

  1. 1.

From
$\boldsymbol{\mathsf{Y}}=\mathcal{A}_{n}(d,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})$,
obtain the reconstruction function $f$ as shown in Section 2.4, and let the
reconstructed version of $\boldsymbol{\mathsf{X}}$ is denoted by
$\boldsymbol{\mathsf{X}}_{\mathcal{A},d}$.

  2. 2.

Fix any one-to-one transformation $T:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$,
and obtain
$\boldsymbol{\mathsf{Z}}=(\boldsymbol{\mathsf{z}}_{1},\dots\boldsymbol{\mathsf{z}}_{n})$,
(the $n$ transformed datapoints) such that
$T((\boldsymbol{\mathsf{X}}_{\mathcal{A},d})_{i})=\boldsymbol{\mathsf{z}}_{i}$
for any $i=1,2,\dots n$.

  3. 3.

Create new dataset
$\tilde{\boldsymbol{\mathsf{X}}}=\boldsymbol{\mathsf{Z}}+(\boldsymbol{\mathsf{X}}-\boldsymbol{\mathsf{X}}_{\mathcal{A},d})$.

  4. 4.

By means of Theorem 4, the output of the algorithm $\mathcal{A}_{n}$ applied
on this new dataset i.e.
$\mathcal{A}_{n}(d,\tilde{\boldsymbol{\mathsf{X}}})(\tilde{\boldsymbol{\mathsf{X}}})$
should be simply a translated, uniformly scaled and rotated version of
$\boldsymbol{\mathsf{Y}}$.

Again, one can use Procrustes’s problem (see Appendix A.3) to measure the
discrepancy and obtain the consistency index, which measures the variability
of the output.

######  Definition 4 (Consistency Index).

The Consistency Index for any algorithm $\mathcal{A}_{n}$ for any dataset
$\boldsymbol{\mathsf{X}}$ at intrinsic dimension $d$, is given by,

$\text{CI}(\mathcal{A}_{n},\boldsymbol{\mathsf{X}},d)=\max_{T\text{ is
invertible}}\min_{\boldsymbol{\mu},\lambda,\boldsymbol{P}}\|\mathcal{A}_{n}(d,\tilde{\boldsymbol{\mathsf{X}}})(\tilde{\boldsymbol{\mathsf{X}}})-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{P}\mathcal{A}_{n}(d,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})\|_{2}^{2}$

where
$\tilde{\boldsymbol{\mathsf{X}}}=T(\boldsymbol{\mathsf{X}}_{\mathcal{A},d})+(\boldsymbol{\mathsf{X}}-\boldsymbol{\mathsf{X}}_{\mathcal{A},d})$.
Similar to the trustability index, lower values of consistency index implies
that algorithm is more consistent in its outputs.

######  Example 7.

Let us consider another simple dimensionality reduction algorithm, which
always outputs the same constant, i.e.
$\phi_{j}(\boldsymbol{\mathsf{x}})=\text{constant}$ for any
$\boldsymbol{\mathsf{x}}$. Then, it is possible to achieve the highest
consistency in output (i.e. the consistency index is equal to $0$), however,
the algorithm is not at all trustable, since at $d=p$, the algorithm does not
output a rotated, uniformly scaled and translated version of the original
dataset.

In its generality, it is not possible to compute the consistency Index for an
algorithm, since the maximization happens over all invertible functions $T$.
Rather, we shall look at the maximization subject to the constraint that the
invertible function $T$ lie in some RKHS $\mathcal{H}$ containing
$\mathbb{R}^{p}$-valued functions. Note that, it is not true that all
functions of the RKHS is invertible (for example, if $f\in\mathcal{H}$, then
$(-f)\in\mathcal{H}$ and hence $f+(-f)\equiv 0\in\mathcal{H}$, which is non-
invertible). Also, we would like the invertible functions to be such that it
does not scale the dataset too small or too large. Therefore, we restrict our
attention to the subclass of functions (i.e. a subset of the RKHS
$\mathcal{H}$ generated by a matrix-valued kernel $k$),

$\mathcal{H}^{\ast}=\left\\{f:f\in\mathcal{H},f(\cdot)=\sum_{i=1}^{\infty}\sum_{j=1}^{p}\beta_{ij}k(\cdot,\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j},0\leq\beta_{ij}\leq
1,\sum_{i=1}^{\infty}\sum_{j=1}^{p}\beta_{ij}=1,\text{ and }f\text{ is
invertible}\right\\}.$

Even then, since the index itself depends on the particular algorithm
$\mathcal{A}_{n}$ we choose, unless we have an understanding of how the
algorithm behaves, it is hard to evaluate the maximum value analytically.
Since $\mathcal{H}^{\ast}$ is an infinite set, it is also not computationally
possible to maximize by searching over the whole space $\mathcal{H}^{\ast}$.
Under certain simple conditions, this maximization over all
$T\in\mathcal{H}^{\ast}$ can be reduced to maximization over only a finite set
of functions, which make this problem computationally tractable. The following
theorem proves this fact.

######  Theorem 5.

If
$S(T)=\min_{\boldsymbol{\mu},\lambda,\boldsymbol{P}}\|\mathcal{A}_{n}(d,\tilde{\boldsymbol{\mathsf{X}}})(\tilde{\boldsymbol{\mathsf{X}}})-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{P}\mathcal{A}_{n}(d,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})\|_{2}^{2}$,
where $\tilde{\boldsymbol{\mathsf{X}}}$ is as given in Definition 4, is a
convex function in $T$, then
$\max_{T\in\mathcal{H}^{\ast}}S(T)=\max_{T\in\partial\mathcal{H}^{\ast}}S(T)$,
where

$\partial\mathcal{H}^{\ast}=\left\\{k(\cdot,\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j}:i=1,\dots
n;\ j=1,\dots p\right\\}.$

Since, $\partial\mathcal{H}^{\ast}$ is now a finite set, therefore, we can use
it to compute the consistency Index for any algorithm $\mathcal{A}_{n}$
approximately. Therefore, we have the following approximate version of
consistency Index, called Tractable Consistency Index.

######  Definition 5.

Tractable Consistency Index for any algorithm $\mathcal{A}_{n}$ for any
dataset $\boldsymbol{\mathsf{X}}$ at intrinsic dimension $d$ is given by,

$\text{TCI}(\mathcal{A}_{n},\boldsymbol{\mathsf{X}},d)=\max_{T\in\partial\mathcal{H}^{\ast}}\min_{\boldsymbol{\mu},\lambda,\boldsymbol{P}}\|\mathcal{A}_{n}(d,\tilde{\boldsymbol{\mathsf{X}}})(\tilde{\boldsymbol{\mathsf{X}}})-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{P}\mathcal{A}_{n}(d,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})\|_{2}^{2}$

where
$\tilde{\boldsymbol{\mathsf{X}}}=T(\boldsymbol{\mathsf{X}}_{\mathcal{A},d})+(\boldsymbol{\mathsf{X}}-\boldsymbol{\mathsf{X}}_{\mathcal{A},d})$,
and
$\partial\mathcal{H}^{\ast}=\left\\{k(\cdot,\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j}:i=1,\dots
n;\ j=1,\dots p\right\\}$ where $k(\cdot,\cdot)$ is a positive definite
matrix-valued kernel function from $\mathbb{R}^{p}\times\mathbb{R}^{p}$ to
$\mathbb{R}^{p}\times\mathbb{R}^{p}$.

##  3 Localized Skeletonization and Dimensionality Reduction (LSDR)

###  3.1 Manifold Approximation Graph

As discussed by Maaten and Hinton [50], one cannot hope to preserve all
pairwise distances in the lower dimensional space, since the same
$\epsilon$-ball in $\mathbb{R}^{p}$ can hold more equidistant points than in
$\mathbb{R}^{d}$ (with $d<p$). This phenomenon is popularly referred to as the
“Crowding problem”. Therefore, the authors argue that only preserving the
local distances can be achieved. A visualisation of this argument is given in
Figure 8(b), where the two red points shown have very different distances in
Euclidean sense and in the distance through the manifold sense. As we notice,
this distance through the manifold can be measured by connecting several
smaller localized distances, where a point is connected to its local
neighbourhood and only these small distances are summed up to calculate the
larger distances. Hence, if one can obtain a graph $G$ connecting these
datapoints by edges whenever two points are in close proximity, then the
corresponding graph would be an approximation of the underlying manifold.
Consequently, the distance between any two points through the manifold can be
computed by the length of the shortest path between those two points through
the graph. Delaunay tessellation and Minimum Cost Spanning Tree (MCST) are two
integral components of creation of such approximation graphs.

(a) Delaunay triangulation for points on a circle

(b) Euclidean distance vs Manifold distance

(c) Euclidean MCST

(d) Manifold approximating graph $G_{\alpha}(P)$

Figure 8: Effect of Delaunay triangulation, Minimum cost spanning tree, and
the manifold approximating graph $G_{\alpha}(P)$

Delaunay triangulation (also known as a Delone triangulation, in short DT) for
a given set $P$ of discrete points in a general position on a plane, is a
triangulation $DT(P)$ such that no point in $P$ is inside the circumcircle of
any triangle in $DT(P)$. It can be shown that this triangulation avoids narrow
triangles and maximizes the sum of the minimum angles of all triangles [9],
hence provide a good triangular approximation of the underlying mesh. Delaunay
[10] extended the idea of triangulation into tessellation which can be
obtained for any set of points $P$ lying in general position in any
dimensional Euclidean vector space. Figure 8(a) shows one such example of DT
for points distributed on a circle on $\mathbb{R}^{2}$. One could simply use
DT as the required graph $G$, but there are some longer edges that connect
points lying on the opposite parts of the circle as DT outputs a graph with
only triangles. However, if Delaunay triangulation is modified to have
polygonal tessellation as well, then such tessellation would allow us to
create the required graph.

Let, $P=\\{\boldsymbol{x}_{1},\dots\boldsymbol{x}_{n}\\}$ be a set of points
lying in the Euclidean space $\mathbb{R}^{p}$, and its Delaunay triangulation
is denoted as $DT(P)=(G_{DT}(P),E_{DT}(P))$. Consider a point
$\boldsymbol{x}_{i}$ and let $e_{1},\dots e_{k}\in E_{DT}(P)$ be the
associated edges with one endpoint equal to $\boldsymbol{x}_{i}$. Then,
considering a multivariate normal distribution centered at
$\boldsymbol{x}_{i}$ with variance $\sigma I_{p}$, one can show that
$\dfrac{\|e_{j}\|^{2}}{\sigma^{2}}$ has a $\chi^{2}_{p}$ distribution.
Similarly,
$\dfrac{1}{\sigma^{2}}\sum_{j=1}^{k}\|e_{j}\|^{2}\sim\Gamma(kp/2,1/2)$
distribution. Therefore, one can show that

$\dfrac{\|e_{j}\|^{2}}{\sum_{j^{\prime}=1}^{k}\|e_{j^{\prime}}\|^{2}}\sim\text{Beta}\left(\dfrac{p}{2},\dfrac{(k-1)p}{2}\right)$

Hence, one can choose some $\alpha\approx 0.95$ and use the above statistic
and compare it with the threshold
$\text{Beta}\left(\frac{p}{2},\frac{(k-1)p}{2},\alpha\right)$ which denotes
the $\alpha$-th quantile of
$\text{Beta}\left(\frac{p}{2},\frac{(k-1)p}{2}\right)$ distribution, and
remove the particular edge $e_{j}$ if the statistic exceeds this threshold.
This is similar to a statistical hypothesis testing setup where the null
hypothesis assumes that the points in the locality of $\boldsymbol{x}_{i}$
which are joined by DT follows a multivariate normal distribution. When the
null hypothesis is false, the high value of the test statistic would indicate
the presence of a longer edge $e_{j}$ in DT, hence the removal of this edge
make the graph $DT(P)$ a better approximation to the manifold.

The edge removal process will only keep the small edges in $DT(P)$. However,
if there are two or more clusters present in the data, then some global
structure or characteristic need to be retained (Otherwise, similar problems
as described in Section 1.2 and Figure 5 may occur). This global structure can
only be preserved through some large edges, which may be deleted in the
aforementioned edge removal process. In this case, the Euclidean minimum cost
spanning tree (MCST) provides us with useful insights. Euclidean MCST is a
subgraph of the Delaunay tessellation $DT(P)$ such that it is a spanning tree
with the sum of its edge lengths equal to the minimum possible value. The
Euclidean MCST for the same circular data is shown in Figure 8(c). The
following result provides a theoretical justification of why MCST is a
reasonable choice to preserve global structure between several clusters.

######  Theorem 6.

Let us assume we have two groups of points given by,

$\boldsymbol{x}=\boldsymbol{\mu}_{1}\boldsymbol{1}_{C_{1}}+\boldsymbol{\mu}_{2}\boldsymbol{1}_{C_{2}}+\epsilon$

such that $|C_{1}|=m,|C_{2}|=n$ and the errors $\epsilon$ are independently
distributed. Assume that for any $\delta>0$,

$P(\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\|_{2}>\delta)<1$

for any $\boldsymbol{x}_{i},\boldsymbol{x}_{j}$ lying in the same class. Then,
the probability that the euclidean MCST of these set of points has two edges
connecting different classes tends to zero as $m,n\rightarrow\infty$, namely
the MCST would have two subtrees, each of them being a spanning tree for each
of the class, and there will be only one edge that connects two classes.

######  Remark 2.

The condition required for Theorem 6 is satisfied if $\epsilon$ has a
continuous distribution with support containing $0$. It simply requires that
the points within the same cluster can be probabilistically arbitrarily close,
which if not met, indicates that the cluster itself can be broken into more
than one sub-clusters. Also, Theorem 6 holds for more than two clusters as
well, where between any two clusters there is at most a single edge of MCST.

Theorem 6 shows that the edges which are part of the MCST are extremely
crucial to approximate the underlying manifold of a point set $P$. Therefore,
we start with the Delaunay triangulation $DT(P)$ of the point set $P$, then
perform the edge removal process by comparing the edge length based test
statistic with the cut-off from beta distribution, and only remove the edge it
is not a part of the MCST. This whole procedure yields us a graph
$G_{\alpha}(P)$ for the point set. An example of such $G_{\alpha}(P)$ is shown
in Figure 8(d). We call such a graph as $\alpha$-approximation of a manifold
generated by the point set $P$.

###  3.2 Skeletonization

The next step in our algorithm is to obtain a skeletal representation of the
manifold from its graph approximation. This skeletal representation would
allow one to obtain points lying on the underlying low-dimensional manifold
(i.e. without any errors). In mathematical terms, a skeletal representation of
a point set
$P=\\{\boldsymbol{x}_{1},\boldsymbol{x}_{2},\dots\boldsymbol{x}_{n}\\}$ is
given by a subset of points
$Sk(P)=\\{\boldsymbol{x}_{i_{1}},\dots\boldsymbol{x}_{i_{k}}\\}$ such that
$\boldsymbol{x}_{i_{j}}=f(\phi_{1}(\boldsymbol{x}_{i_{j}}),\dots,\phi_{d}(\boldsymbol{x}_{i_{j}}))$
i.e. $\epsilon(\boldsymbol{x}_{i_{j}})=0$ for all $j=1,2,\dots k$ (see
Definition 1). Let us now consider a single point
$\boldsymbol{x}=f(\phi_{1}(\boldsymbol{x}),\dots\phi_{d}(\boldsymbol{x}))+\epsilon(\boldsymbol{x})$
where $\boldsymbol{x}\in\mathbb{R}^{p}$. Let $\boldsymbol{v}\in\mathbb{R}^{p}$
be a vector such that
$\boldsymbol{v}^{\intercal}\nabla\phi_{j}(\boldsymbol{x})=0$ for all
$j=1,2,\dots d$, i.e. $v$ is locally orthogonal to each of the embedding
functions. Applying Taylor’s theorem, we then get
$\phi_{j}(\boldsymbol{x}+\boldsymbol{v})\approx\phi_{j}(\boldsymbol{x})+\boldsymbol{v}^{\intercal}\nabla\phi_{j}(\boldsymbol{x})=\phi_{j}(\boldsymbol{x})$,
i.e. the embedding remains approximately constant on the line
$\boldsymbol{x}+\lambda\boldsymbol{v}$ where $\lambda\in\mathbb{R}$ is a
parameter. Observe that, in this case,
$\epsilon(\boldsymbol{x}+\lambda\boldsymbol{v})\approx\epsilon(\boldsymbol{x})+\lambda\boldsymbol{v}$,
which shows that the error is locally linear on that orthogonal line. This
realization then allows us to define boundaries of the manifold in
mathematical terms, a point $\boldsymbol{x}_{b}$ is a boundary point if either
$\boldsymbol{x}_{b}=\displaystyle\max_{\boldsymbol{x}+\lambda\boldsymbol{v}\in
P}(\boldsymbol{x}+\lambda\boldsymbol{v})$ or
$\boldsymbol{x}_{b}=\displaystyle\min_{\boldsymbol{x}+\lambda\boldsymbol{v}\in
P}(\boldsymbol{x}+\lambda\boldsymbol{v})$ for some $\boldsymbol{x}\in Sk(P)$.
The set of all boundary points are denoted as $B(P)$.

Now consider a point $\boldsymbol{x}\in Sk(P)$ and let $\boldsymbol{x}_{b}\in
B(P)$ be a corresponding boundary point. Assume that $\boldsymbol{x}_{b}$ is a
part of a $(p-1)$-dimensional hyperplane $\mathcal{P}$ constituted of some
edges in the graph $G_{\alpha}(P)$ obtained from Delaunay tessellation. If
this hyperplane is part of two $(p+1)$-simplices $S_{1}$ and $S_{2}$, they
must lie in the opposite side of $\mathcal{P}$, therefore, the ray parallel to
the vector $(\boldsymbol{x}_{b}-\boldsymbol{x})$ starting from
$\boldsymbol{x}$ must intersect $S_{1}$ and $S_{2}$ both. Clearly, for one of
this intersection point, the distance from the intersection to
$\boldsymbol{x}$ must be larger than $\|\boldsymbol{x}_{b}-\boldsymbol{x}\|$
which contradicts the extremality of the boundary point. Therefore, we have
the following result.

######  Theorem 7.

Any boundary point of the manifold in $\mathbb{R}^{p}$ must lie on a $(p-1)$
dimensional hyperplane which is a part of atmost one $(p+1)$-simplex.

Due to Theorem 7, it is possible to obtain the boundary points by simply
verifying the condition for each vertex in $G_{\alpha}(P)$. Figure 9(a) shows
the boundary of $G_{\alpha}(P)$ computed in this way for a point set arranged
in a spiral. However, the significance of this procedure is that it enables us
to detect the boundary without any reference to the skeletal representation,
hence one can obtain skeletal representation from the boundary itself. To see
this, note that the skeletal points are given by those where the error is
equal to $0$, and the boundary points are essentially where the error attains
its maximum value in magnitude. Since the error is locally linear as we have
shown, one can move away from the boundary inside the point cloud which should
necessarily decrease the error. Therefore, the point that is the most distance
away from the boundaries should have the smallest error in magnitude, which
may serve our purpose as a skeletal point.

(a) $G_{\alpha}(P)$ and its boundary

(b) $G_{\alpha}(P)$ and its skeletal representation

Figure 9: A point set $P$ in the form of a one dimensional manifold in a 2D
plane, its approximation graph $G_{\alpha}(P)$ with its boundary and skeletal
representation.

In order to put these ideas into a working algorithm, for every point
$\boldsymbol{x}\in P$, we compute its distance from the boundary as

$d_{G}(\boldsymbol{x},B(P))=\min_{\boldsymbol{x}_{b}\in
B(P)}d_{G}(\boldsymbol{x},\boldsymbol{x}_{b})$

where $d_{G}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$ is the sum of the length
of edges in the shortest path between $\boldsymbol{x}_{i}$ and
$\boldsymbol{x}_{j}$ through the graph $G_{\alpha}(P)$. Then, we iterate over
all points $\boldsymbol{x}\in P$ again and mark them as a skeletal point if

$d_{G}(\boldsymbol{x},B(P))=\max_{\boldsymbol{w}\in
NN_{k}(\boldsymbol{x})}d_{G}(\boldsymbol{w},B(P))$

where $NN_{k}(\boldsymbol{x})$ denotes the set of $k$-nearest neighbours of
$\boldsymbol{x}$. When $k$ is sufficiently small, the above criterion assures
that the identified skeletal points lie in the more inner region than their
neighbours. This idea of the skeleton as the most distant points from its
boundary dates back to Blum, 1967 [3] and is later proved to have many
optimality criterion [44]. Such a skeletal representation is shown in Figure
9(b) where the skeletal points are identified with a larger size and different
colour. As seen from Figure 9(b), the skeleton points lie approximately in the
middle of the one-dimensional spiral manifold.

###  3.3 Multidimensional Scaling and Kernel embedding

Multidimensional scaling (MDS) attempts to find a low dimensional
representation $\\{\boldsymbol{y}_{1},\dots\boldsymbol{y}_{n}\\}$ with each
$\boldsymbol{y}_{i}\in\mathbb{R}^{d}$ such that the usual Euclidean distance
between $\boldsymbol{y}_{i}$ and $\boldsymbol{y}_{j}$ approximately equals
$q_{ij}$ for all $1\leq i<j\leq n$, for pre-specified choice of $q_{ij}$ [4].
Therefore, for pre-specified values of the distances, it attempts to find
points in a $d$-dimensional space that preserves the distances as much as
possible, by minimizing the stress function

$\text{Stress}_{D}(\boldsymbol{y}_{1},\dots\boldsymbol{y}_{n})=\left(\sum_{1\leq
i<j\leq
n}(q_{ij}-\|\boldsymbol{y}_{i}-\boldsymbol{y}_{j}\|)^{2}\right)^{1/2}.$

This variation is known as the metric MDS. Since the distances are invariant
under translation and rotational transformation of $\boldsymbol{y}_{i}$’s, the
solution to metric MDS is not unique. However, one can show that a condition
similar to the desired values of trustability index mentioned in Section 2.5
can be achieved by metric MDS.

######  Theorem 8.

Let us assume for a dataset
$\\{\boldsymbol{x}_{1},\dots\boldsymbol{x}_{n}\\}$, the metric MDS is applied
on $q_{ij}=\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\|$ to obtain an embedding
$\boldsymbol{y}_{i}\in\mathbb{R}^{p}$. Then,
$\boldsymbol{y}_{i}=\boldsymbol{v}+\boldsymbol{P}\boldsymbol{x}_{i}$ for all
$i=1,2,\dots n$ for some $\boldsymbol{a}\in\mathbb{R}^{p}$ and some orthogonal
matrix $\boldsymbol{P}$.

Theorem 8 shows the connection between metric MDS and the trustability index.
In other words, it says that as long as the original distances can be
recovered in some proposed DR algorithm, a metric MDS step then ensures that
the DR algorithm remains optimal with respect to the trustability index.
However, the distances which are needed to pass into the metric MDS step is
the Euclidean distance. In contrast, if we pass an approximation to the
geodesic distance in the metric MDS setup, then one can show under certain
reasonable assumptions that the corresponding DR algorithm will be optimal
with respect to the consistency index defined earlier.

######  Theorem 9.

Let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ be a smooth co-ordinate
function of a Riemannian manifold $(\mathcal{M},g)$ such that a unit-velocity
curve in $\mathbb{R}^{d}$ induces a constant velocity curve in $\mathcal{M}$.
Let, $\phi:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ be a local diffeomorphism
such that $\phi\circ f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{p}$ be the
coordinate function of another Riemannian manifold
$(\mathcal{M}^{\prime},g^{\prime})$ with

$g(v,w)=cg^{\prime}(\phi^{\ast}v,\phi^{\ast}w),\qquad\forall
v,w\in\mathcal{M}$

where $\phi^{\ast}$ is the push-forward and $c$ is some positive real
constant. Then, for any
$\\{\boldsymbol{\theta}_{1},\dots\boldsymbol{\theta}_{n}\\}$ with each
$\boldsymbol{\theta}_{i}\in\mathbb{R}^{d}$, the $d$-dimensional output
$\\{\boldsymbol{y}_{1},\dots\boldsymbol{y}_{n}\\}$ of metric MDS of
$\\{f(\boldsymbol{\theta}_{1}),\dots f(\boldsymbol{\theta}_{n})\\}$ and the
$d$-dimensional output $\\{\boldsymbol{z}_{1},\dots\boldsymbol{z}_{n}\\}$ of
metric MDS of $\\{\phi\circ f(\boldsymbol{\theta}_{1}),\dots\phi\circ
f(\boldsymbol{\theta}_{n})\\}$ are related by

$\boldsymbol{z}_{i}=\boldsymbol{v}+\lambda\boldsymbol{Py}_{i},\qquad
i=1,2,\dots n$

for some fixed $\boldsymbol{v}\in\mathbb{R}^{p},\lambda\in\mathbb{R}$ and some
orthogonal matrix $\boldsymbol{P}$.

From the skeletal points obtained in Section 3.2, one may assume that these
points lie on the low dimensional manifold in $\mathbb{R}^{p}$. The
assumptions made in Theorem 9 enables one to approximate the geodesic distance
on this manifold by summing up local Euclidean distances through the
approximation graph $G_{\alpha}$. Therefore, if one can use the length of the
shortest paths between two skeletal points on $G_{\alpha}$ as an approximation
to the geodesic distance, then Theorem 9 ensures that the output of the metric
MDS of these distances would ensure a dimension reduction for the skeletal
points and that DR would be optimal with respect to the consistency index.
This is precisely what we will do to find the low dimensional embedding of the
skeletal points.

However, the metric MDS would be unable to properly find low dimensional
embedding for the non-skeletal points, since they do not lie exactly on the
manifold. Therefore, given that the skeletal points sufficiently cover the
entire structure of the manifold, one may use the technique mentioned in
Section 2.3 to obtain a low-dimensional embedding of non-skeletal points.
While metric MDS preserves the global structure as much as possible by
considering only the skeletal points, the local structure is recovered by the
kernel embedding technique, hence the whole procedure is expected to balance
the local and global structures found in the high dimensional data nicely.

Data: $\boldsymbol{x}_{1},\dots\boldsymbol{x}_{n}$ where each
$\boldsymbol{x}_{i}\in\mathbb{R}^{p}$

Input: $\alpha,k$, embedding dimension $d$ and some kernel $K(\cdot,\cdot)$

Output: $\boldsymbol{y}_{1},\dots\boldsymbol{y}_{n}$ where each
$\boldsymbol{y}_{i}\in\mathbb{R}^{d}$ ($d<<p$)

begin Step 1

Let $P=\\{\boldsymbol{x}_{1},\dots\boldsymbol{x}_{n}\\}$;

Compute $p$-dimensional Delaunay tessellation $DT(P)$;

Compute Euclidean MCST of the point set $P$;

foreach  _point $\boldsymbol{x}_{i}$_ do

Let $e_{1},\dots e_{d(\boldsymbol{x}_{i})}$ be the edges with one endpoint
equal to $\boldsymbol{x}_{i}$ where $d(\boldsymbol{x}_{i})$ is the degree of
$\boldsymbol{x}_{i}$;

foreach  _edge $e_{j}$_ do

if  _$e_{j}\notin E(MCST(P))$_ then

Compute
$T_{j}=\dfrac{\|e_{j}\|^{2}}{\sum_{j^{\prime}=1}^{d(\boldsymbol{x}_{i})}\|e_{j^{\prime}}\|^{2}}$;

if  _$T_{j} >\text{Beta}(p/2,(d(\boldsymbol{x}_{i})-1)/2,\alpha)$_ then remove
$e_{j}$ from the graph ;

end if

end foreach

end foreach

Call the resulting graph $G_{\alpha}(P)$;

Use Dijkstra’s algorithm to calculate all graph-based pairwise distances;

Let $d_{G}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})$ be the graph-based distance
between $\boldsymbol{x}_{i}$ and $\boldsymbol{x}_{j}$;

end

begin Step 2

Intialize $B(P)=\\{\\}$;

foreach  _edges $e\in E(G_{\alpha}(P))$_ do

$Count=$ Number of $(p+1)$-simplices with $e$ as an edge;

if  _$\text{Count}\leq 1$_ then insert edge $e$ into $B(P)$;

end foreach

foreach  _point $\boldsymbol{x}_{i}$_ do

Compute $d_{B}(\boldsymbol{x}_{i})=\min_{\boldsymbol{x}\in
B(P)}d_{G}(\boldsymbol{x}_{i},\boldsymbol{x})$;

end foreach

Initialize $Sk(P)=\\{\\}$;

foreach  _point $\boldsymbol{x}_{i}$_ do

Calculate the $k$-nearest neighbours $NN_{k}(\boldsymbol{x}_{i})$ according to
$d_{G}$ metric;

if  _$d_{B}(\boldsymbol{x}_{i})=\max_{\boldsymbol{v}\in
NN_{k}(\boldsymbol{x}_{i})}d_{B}(\boldsymbol{v})$_ then Mark
$\boldsymbol{x}_{i}$ as skeletal point and add it to $Sk(P)$;

end foreach

end

begin Step 3

Let $Sk(P)=\\{\boldsymbol{x}_{i_{1}},\dots\boldsymbol{x}_{i_{r}}\\}$;

Apply metric MDS on the distance matrix $Q$ with entries
$q_{jj^{\prime}}=d_{G}(\boldsymbol{x}_{i_{j}},\boldsymbol{x}_{i_{j^{\prime}}})$
to obtain $d$-dimensional representation
$\boldsymbol{y}_{i_{1}},\dots\boldsymbol{y}_{i_{r}}$;

foreach  _$i=1,2,\dots n$_ do

$\boldsymbol{y}_{i}=\sum_{j=1}^{r}\boldsymbol{y}_{i_{j}}K(\boldsymbol{x}_{i},\boldsymbol{x}_{i_{j}})/\sum_{j=1}^{r}K(\boldsymbol{x}_{i},\boldsymbol{x}_{i_{j}})$;

end foreach

end

Algorithm 1 Localized Skeletonization and Dimensionality Reduction (LSDR)

The final algorithm, which we call “Localized Skeletonization and Dimensional
Reduction” (LSDR) is described in Algorithm 1. While the algorithm in its
generality is described by three parameters $\alpha,k$ and the kernel of
choice, we found that the standard choice of $\alpha\approx 0.9$ or $0.95$
(the usual choice of statistical significance), $k=3$ works pretty well. The
required kernel could be chosen as the Gaussian kernel with a properly chosen
bandwidth in some data-dependent way. The choice of the bandwidth $\sigma$
should be such that if the one considers a $\sigma$-ball around the skeletal
points, no point should lie outside of them. In other words, one
recommendation could be to take

$\sigma=\max_{\boldsymbol{x}_{i}\in
Sk(P)}\sqrt{d_{B}^{2}(\boldsymbol{x}_{i})+\min_{\boldsymbol{x}_{j}\in
Sk(P),i\neq j}d_{G}^{2}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})}$

##  4 Applications

###  4.1 Simulation Studies

####  4.1.1 Performance indices

In this section, we shall compute the trustability index and the Tractable
consistency index for various dimensionality reduction algorithms. In
particular, we consider the algorithms Principal Component Analysis (PCA),
Kernel Principal Component Analysis (kPCA), Locally Linear Embedding (LLE),
Hessian Locally Linear Embedding (HLLE), Laplacian Eigenmaps (LE),
t-distributed Stochastic Neighbour Embedding (tSNE), Uniform Manifold
Approximation and Projection (UMAP) etc. To compute the indices, we shall
consider various artificial and real-life datasets as discussed below.

  1. 1.

In simulation setup S1, three clusters of datapoints are generated, with each
cluster of datapoints following $10$-dimensional Gaussian distribution with a
different location and common covariance matrix.

  2. 2.

In simulation setup S2 the data has been generated uniformly from the
$10$-dimensional hypercube $[0,1]^{10}$.

  3. 3.

In setup S3 the datapoints lie approximately on the $2$-D surface of the unit
sphere in $\mathbb{R}^{3}$.

  4. 4.

In simulation setup S4, the datapoints lie approximately on a swiss roll
manifold in $\mathbb{R}^{3}$.

  5. 5.

Sonar data is a popular benchmark dataset for dimensionality reduction [16].
The dataset contains $60$ variables (each between $0$ and $1$) representing
the angular deviations of different signals when bouncing off a metal or rock
object.

  6. 6.

WBCD or Wisconsin Breast Cancer Dataset from UCI Machine Learning Repository
[12] contains $32$ attributes of $569$ patients with tumour, including a class
variable indicating whether the tumour is malignant or benign.

  7. 7.

COIL2000 data set [52] used in the COIL 2000 Challenge contains information on
customers of an insurance company. The data consists of 86 variables and
includes product usage data and socio-demographic data.

Since generally the trustability index increases as the number of datapoints
$n$, it is meaningful to look at
$\text{TI}(\mathcal{A}_{n},\boldsymbol{\mathsf{X}})/n$. In each of the
simulation setups S1, S2, S3, S4, we consider $n=200$ datapoints. The
normalized trustability and tractable consistency indices have been computed
for the aforementioned datasets and for the aforementioned algorithms, which
have been summarized in Table 1 and Table 2 respectively. As seen from these
tables, PCA, LLE, Hessian LLE, and Laplacian Eigenmaps provides very trustable
outputs, however, they are less consistent. Therefore, their output can be
used for subsequent analysis of regression and classification purposes, since
these algorithms do not create new patterns in the data. On the other hand,
tSNE and UMAP are less trustable but more consistent. Their output should
primarily be used for visualisation purposes, but not for subsequent analysis
if trustworthiness is a requirement. These algorithms might create new
superficial patterns in the data which were not originally present. Also,
Laplacian Eigenmaps is generally better than Hessian LLE for both trustability
and consistency, hence should be preferred for any kind of data. Finally, our
proposed algorithm LSDR is found to be extremely trustable, but without
compromising its consistency index to much extent.

Dataset | PCA | kPCA | LLE | HLLE | LE | tSNE | UMAP | LSDR
---|---|---|---|---|---|---|---|---
S1 | 0 | 138.727 | 10.112 | 0.084 | 0.083 | 128.366 | 91.807 | 0.995
S2 | 0 | 168.667 | 9.997 | 0.822 | 0.018 | 32.007 | 4.885 | 1.037
S3 | 0 | 32.781 | 2.991 | 0.011 | 0.014 | 4.917 | 11.701 | 0.012
S4 | 0 | 112.678 | 3 | 0.021 | 0.034 | 142.131 | 60.023 | 0.159
Sonar | 0 | 56.897 | 59.984 | 59.983 | 0.288 | 2.771 | 9.826 | 0.944
WBCD | 0 | 30.867 | 30.021 | 33.309 | 0.052 | 39.205 | 64.206 | 1.746
COIL2000 | 0 | 118.322 | 86.001 | 92.116 | 0.172 | 11.819 | 38.311 | 2.193
Table 1: Normalized Trustability Index ($n^{-1}\text{TI}(\mathcal{A}_{n},\boldsymbol{\mathsf{X}})$) of various algorithms on various datasets Dataset | PCA | kPCA | LLE | HLLE | LE | tSNE | UMAP | LSDR
---|---|---|---|---|---|---|---|---
S1 | 58434.22 | 9.31 | 232.511 | 203.115 | 153.313 | 23.963 | 11.818 | 6.93
S2 | 143.233 | 2.946 | 6.305 | 8.818 | 8.003 | 3.116 | 2.291 | 2.445
S3 | 67.414 | 10.065 | 15.337 | 29.212 | 29.252 | 13.535 | 2.904 | 4.196
S4 | 72.212 | 5.193 | 34.611 | 31.995 | 27.212 | 3.998 | 3.184 | 1.008
Sonar | 58.414 | 2.854 | 18.884 | 21.229 | 18.818 | 3.971 | 4.498 | NA
WBCD | 1.234$\times 10^{6}$ | 3.214 | 48.913 | 74.227 | 70.212 | 2.105 | 2.816 | NA
COIL2000 | 109783.7 | 4.009 | 69.913 | 75.200 | 71.200 | 3.211 | 2.877 | NA
Table 2: Normalized Tractable Consistency Index
($n^{-1}\text{TCI}(\mathcal{A}_{n},\boldsymbol{\mathsf{X}})$) of various
algorithms (Due to computational complexity limitations, we could not compute
consistency index of LSDR for Sonar, WBCD, COIL2000 data; see Section 5 for
details)

####  4.1.2 Performance of LSDR

We have performed various simulation to compare the performance of LSDR with
that of tSNE and UMAP, the two state of the art dimensionality reduction
algorithms. To compare them on common ground, simulated setups are chosen from
the demonstrations by Google engineers [54, 7] who have provided an extensive
analysis of tSNE and UMAP over different sorts of simulated data. They have
also shown the effect of several hyperparameters of these algorithms and how
they behave. Due to limited time and space constraints, we shall discuss only
a few here. More discussion of other simulation setups is provided in the
Appendix.

Our first simulation setup is the popular spiral data in 2-dimension, which is
governed by a single parameter $\theta$ denoting the slope of a point with
respect to the origin. The original data, its one-dimensional embeddings by
tSNE, UMAP and LSDR as a function of the true $\theta$ is shown in Figure 11.
As seen from the plots, tSNE and UMAP, even under different parameter values
recover $\theta$ only partially, hence results in discontinuous graphs. On the
other hand, LSDR outputs a continuous curve even under different choices of
its parameters and is extremely well correlated with the true value of
$\theta$. However, there is a kernel effect at the tail of the manifold, where
the correlation between $\theta$ and the embedding weakens severely.

(a) Original 3D Wooly Mammoth Data

(b) 2D Embedding by tSNE of Wooly Mammoth Data

(c) 2D Embedding by UMAP of Wooly Mammoth Data

(d) 2D Embedding by LSDR of Wooly Mammoth Data

Figure 10: Analysis of Wooly Mammoth Data

(a) Original Spiral Dataset

(b) 1-dimensional embedding by tSNE vs the true parameter for different
perplexity

(c) 1-dimensional embedding by UMAP vs the true parameter for different
neighbours

(d) 1-dimensional embedding by LSDR vs the true parameter for different NN
parameter $k$

Figure 11: Analysis of Spiral Data

Another interesting property of LSDR is that it preserves more global
structure than tSNE and UMAP. For example, in the usual Gaussian cluster data
[54] where there are three clusters with different distances between them,
tSNE usually fails to distinguish them in small perplexity parameter. As the
perplexity grows, three clusters separate from each other, however, the one-
dimensional embedding never shows that there are regions where no data is
available. Whereas UMAP overcomes this problem and shows a blank region
between clusters, however, the clusters become equidistant despite being
different distances away from each other. As shown in Figure 12, LSDR
preserves this global structure nicely.

(a) Original Spiral Dataset

(b) 1-dimensional embedding by tSNE for different perplexity

(c) 1-dimensional embedding by UMAP for different neighbours

(d) 1-dimensional embedding by LSDR for different bandwidths $\sigma$ in
radial kernel

Figure 12: Analysis of Three Cluster Data

Another more close to real-life example is the Wooly Mammoth Data provided by
Google PAIR team [7]. The original dataset contains points from the surface of
a 3D Wooly mammoth, and a possible best realization of its 2D embedding would
be its skin. As shown in Figure 10, both tSNE and UMAP retains the local
structure of the data extremely well, however, some parts of the embedding
becomes discontinuous. It is quite clear that if the mammoth is skinned and
its skin is put flat onto a plane, it would be continuous with the parts
corresponding to the legs and trunks sticking out. The 2D embedding obtained
by LSDR seems to outperform tSNE and UMAP in this direction by preserving the
continuity as well as depicting a very relevant output. LSDR only fails to
position the head properly between two frontal legs, which should have been
exactly in the middle of those two legs in an ideal dimensionality reduction.

More details about various aspects of comparison between tSNE, UMAP and LSDR
are detailed in the Appendix C.

##  5 Conclusions and Future Aspects

Analysis of the Trustability index and the Tractable consistency index of
various algorithms provide insights into how these algorithms behave in
general. In essence, the output of tSNE or UMAP may be used in general for
visualisation purposes. However, if the post-reduction analysis requires
classification or clustering on the reduced data, then Laplecian Eigenmaps
provide a good choice of balance. Also, since the graphical algorithms are
generally more trustable, one should use graph-based approximation to the
manifold to preserve the global structure, and this primarily gave the
intuition of LSDR.

Turning to our proposed algorithm, LSDR outperforms both tSNE and UMAP in
terms of its preservation of global structures. Use of metric MDS with the
manifold approximation graph $G_{\alpha}$ provides us with near-optimality
properties with respect to both the trustability and consistency indices (see
Theorem 8 and 9). LSDR preserves the global structure by using skeletal
representation, while the local structure is provided by the kernel embedding
based usage. In a sense, LSDR is similar to SVM, where the skeletal points act
as support points, and any point on the space can be represented in low
dimension using the representations of these support points. One of the major
problem in the trustability of tSNE and UMAP is that they are iterative in
nature, therefore, the final embedding depends a lot on the initialization of
the points. In contrast, our proposed LSDR is not an iterative algorithm,
hence overcomes this difficulty, but retaining the same level of performance
for these setups. The choice of hyperparameters is one of the deciding factors
of the performance of LSDR. For example, if $k$ is chosen very large, then
there will naturally be fewer skeletal points, and consequently, these points
would not be sufficient to describe the manifold’s skeleton properly.
Therefore, a choice of smaller $k$ (about 3-5) and a higher bandwidth for the
kernel is required to ensure that all the points in the dataset are
sufficiently represented.

On the other hand, there are some obvious limitations. The most expensive part
in LSDR is the Delaunay tessellation which has a worst-case complexity of
$O(n^{p/2})$. Therefore, if $p$ is large, it is not computationally feasible
in its current form. One could then proceed by making chunks of the whole
dataset with respect to close proximities and then perform Delaunay
tessellation (this would reduce effective $n$). Another choice could be to use
Gabriel graph or Beta skeleton [44, 8] which are mathematically proven to be a
good approximation to the underlying manifold. However, their implementations
were not readily available in R, hence with the limited scope of time and
space, we could not explore it further. Another problem occurs if there is no
error in the data, i.e. the data lies in an exact $d$-dimensional manifold in
$\mathbb{R}^{p}$. In this case, LSDR would detect every point as a point on
the boundary $B(P)$, and theoretically, every point would become a skeletal
point. While in a realistic scenario, there are always errors, such a case
might occur if $n<p$. In this case, an initial PCA (or metric MDS) step is
needed to be performed to reduce the dimension while preserving the pairwise
distances. Since in a $d$-dimensional plane, only $(d+1)$ distances are needed
to be known to define the position of a point, and since there are $(n-1)$
pairwise distances available, it is possible to embed all these $n$ points in
$\mathbb{R}^{n-2}$ without changing their pairwise distances. After this
initial adjustment, our proposed method LSDR can be used. Finally, LSDR
requires its input to become numerical quantities, hence it cannot reduce the
dimensionality of the points if only a similarity metric is given. For
example, LSDR cannot be used for textual data. However, this problem can again
be solved by using an initial metric MDS step to obtain numerical
representations in $\mathbb{R}^{n-2}$ by preserving the similarities (or
distances).

## Acknowledgements

I would like to thank my supervisor, Dr. Smarajit Bose, ISRU, Indian
Statistical Institute, Kolkata, for his helpful comments and enthusiastic
guidance throughout the project.

##  Appendix A Prerequisites

###  A.1 Reproducible Kernel Hilbert Spaces (RKHS)

Reproducible Kernel Hilbert Spaces (RKHS) is a Hilbert space of functions
which is generated by a basis of kernel functions. If
$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ be a kernel function,
then the RKHS generated by its is the set of functions

$\mathcal{H}=\left\\{\sum_{i=1}^{\infty}\alpha_{i}k(\cdot,x_{i}):x_{i}\in\mathcal{X}\text{
and }\sum_{i=1}^{\infty}\alpha_{i}^{2}k(x_{i},x_{i})<\infty\right\\},$

is a Hilbert space equipped with the inner product

$\left\langle\sum_{i=1}^{\infty}\alpha_{i}k(\cdot,x_{i}),\sum_{j=1}^{\infty}\beta_{j}k(\cdot,y_{j})\right\rangle_{\mathcal{H}}=\sum_{i=1}^{\infty}\sum_{j=1}^{\infty}\alpha_{i}\beta_{j}k(x_{i},y_{j})$

Clearly, $\mathcal{H}$ is a set of functions $f$ that maps elements of
$\mathcal{X}$ to $\mathbb{R}$. For any choice of positive definite kernel
function, the existence of such RKHS is guaranteed by Moore–Aronszajn theorem
[2]. To evaluate any function $f\in\mathcal{H}$ at any $x\in\mathcal{X}$, one
can simply compute $\langle f,k(\cdot,x)\rangle_{\mathcal{H}}$ since,

$\langle
f,k(\cdot,x)\rangle_{\mathcal{H}}=\left\langle\sum_{i=1}^{\infty}\alpha_{i}k(\cdot,x_{i}),k(\cdot,x)\right\rangle_{\mathcal{H}}=\sum_{i=1}^{\infty}\alpha_{i}k(x,x_{i})=f(x).$

This means, the evaluation functional $T_{x}:\mathcal{H}\rightarrow\mathbb{R}$
that maps a function $f\rightarrow f(x)$ is a linear functional. This turns
out to be a defining property of RKHS in the sense that, any Hilbert space
$\mathcal{H}$ on which evaluation functional is a bounded linear functional is
generated by some positive definite kernel $k$ [2, 27].

######  Example 8.

Let $\mathcal{H}$ be all functions $f$ on $\mathbb{R}$ such that the support
of the Fourier transform of $f$ is contained in $[-a,a]$. Then,

$k(x,y)=\dfrac{\sin(a(y-x))}{a(y-x)}$

and $\langle f,g\rangle=\displaystyle\int fg$, defines a RKHS.

######  Example 9.

Another popular example of RKHS is the class of functions
$f:\mathbb{R}^{p}\rightarrow\mathbb{R}$ for any $p$, generated by the radial
basis kernel

$k(\boldsymbol{x},\boldsymbol{y})=\exp\left(-\lambda\|\boldsymbol{x}-\boldsymbol{y}\|^{2}\right),\qquad\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{p}$

######  Example 10.

Finally, an example of non-continuous kernel is the Bregman kernel defined by,

$k(x,y)=\begin{cases}1&\text{ if }x=y\\\ 0&\text{ otherwise }\end{cases}$

This kernel could be a good choice if there are categorical variables in the
dataset. This kernel reproduces a RKHS $\mathcal{H}$ containing all square
integrable real-valued functions.

While the RKHS is introduced in a fixed variable setup, it can generalized to
stochastic setup easily by replacing $x_{i}$’s by $\mathcal{X}$-valued random
variables $\mathsf{x}_{i}:\Omega\rightarrow\mathcal{X}$. Now, if one considers
the linear functional $\mathcal{E}:\mathcal{H}\rightarrow\mathbb{R}$ that maps
$f\in\mathcal{H}$ to $\mathbb{E}(f(\boldsymbol{\mathsf{x}}(\omega))$, and
assumes that
$\mathbb{E}(k(\boldsymbol{\mathsf{x}}(\omega),\boldsymbol{\mathsf{x}}(\omega))<\infty$,
then it follows that $\mathcal{E}$ is a bounded linear functional. Riesz
representation theorem [46] guarantees that there is a unique member
$f_{0}\in\mathcal{H}$ such that $\langle
f,f_{0}\rangle_{\mathcal{H}}=\mathcal{E}(f)$, which one may define as
$\mathbb{E}(k(\cdot,\boldsymbol{\mathsf{x}}(\omega))$.

Analogously, if
$\mathbb{E}(k(\boldsymbol{\mathsf{x}}(\omega),\boldsymbol{\mathsf{x}}(\omega))<\infty$,
one may define the variance operator of $\boldsymbol{\mathsf{x}}$ in
$\mathcal{H}$ as,

$\Gamma_{xx}=\mathbb{E}(k(\cdot,\boldsymbol{\mathsf{x}}(\omega))\otimes
k(\cdot,\boldsymbol{\mathsf{x}}(\omega)))-\mathbb{E}(k(\cdot,\boldsymbol{\mathsf{x}}(\omega)))\otimes\mathbb{E}(k(\cdot,\boldsymbol{\mathsf{x}}(\omega)))$

where $\otimes$ is the symbol for tensor product (see Li [27] for details).
One can further generalize this concept to define covariance operator between
functions lying in two different RKHS. Let, $\mathcal{H}_{x}$ and
$\mathcal{H}_{y}$ be two RKHS generated by two different kernels $k_{x}$ and
$k_{y}$. Then, the covariance operator is similarly given by

$\Gamma_{xy}=\mathbb{E}(k_{x}(\cdot,\boldsymbol{\mathsf{x}}(\omega))\otimes
k_{y}(\cdot,\boldsymbol{\mathsf{y}}(\omega)))-\mathbb{E}(k_{x}(\cdot,\boldsymbol{\mathsf{x}}(\omega)))\otimes\mathbb{E}(k_{y}(\cdot,\boldsymbol{\mathsf{y}}(\omega))).$

This means, for any $f(\boldsymbol{\mathsf{x}})\in\mathcal{H}_{x}$ and
$g(\boldsymbol{\mathsf{y}})\in\mathcal{H}_{y}$,
$\text{Var}(f(\boldsymbol{\mathsf{x}}))=\langle
f,\Gamma_{xx}f\rangle_{\mathcal{H}_{x}}$ and
$\text{Cov}(f(\boldsymbol{\mathsf{x}},g(\boldsymbol{\mathsf{y}})=\langle
f,\Gamma_{xy}g\rangle_{\mathcal{H}_{x}}$. In particular, according to Theorem
12.1 of Li [27], when one has i.i.d. sample
$\boldsymbol{\mathsf{x}}_{1},\dots\boldsymbol{\mathsf{x}}_{n}$ for
$\boldsymbol{\mathsf{x}}$ and
$\boldsymbol{\mathsf{y}}_{1},\dots\boldsymbol{\mathsf{y}}_{n}$ for
$\boldsymbol{\mathsf{y}}$, these operators are estimated as

$\widehat{\Gamma}_{xx}=\dfrac{1}{n}\boldsymbol{G}_{x},\
\widehat{\Gamma}_{xy}=\dfrac{1}{n}\boldsymbol{G}_{x}\boldsymbol{G}_{y},\
\widehat{\Gamma}_{yy}=\dfrac{1}{n}\boldsymbol{G}_{y}$

where the matrix $G_{x}$ and $G_{y}$ have entries given by,

$\displaystyle(G_{x})_{i,j}$
$\displaystyle=k_{x}(\boldsymbol{\mathsf{x}}_{i},\boldsymbol{\mathsf{x}}_{j})-\dfrac{1}{n}\sum_{i}k_{x}(\boldsymbol{\mathsf{x}}_{i},\boldsymbol{\mathsf{x}}_{j})-\dfrac{1}{n}\sum_{j}k_{x}(\boldsymbol{\mathsf{x}}_{i},\boldsymbol{\mathsf{x}}_{j})+\dfrac{1}{n^{2}}\sum_{i}\sum_{j}k_{x}(\boldsymbol{\mathsf{x}}_{i},\boldsymbol{\mathsf{x}}_{j})$
(3) $\displaystyle(G_{y})_{i,j}$
$\displaystyle=k_{y}(\boldsymbol{\mathsf{y}}_{i},\boldsymbol{\mathsf{y}}_{j})-\dfrac{1}{n}\sum_{i}k_{y}(\boldsymbol{\mathsf{y}}_{i},\boldsymbol{\mathsf{y}}_{j})-\dfrac{1}{n}\sum_{j}k_{y}(\boldsymbol{\mathsf{y}}_{i},\boldsymbol{\mathsf{y}}_{j})+\dfrac{1}{n^{2}}\sum_{i}\sum_{j}k_{y}(\boldsymbol{\mathsf{y}}_{i},\boldsymbol{\mathsf{y}}_{j}).$
(4)

######  Definition 6 (Universal Kernels).

A continuous kernel $k(\cdot,\cdot)$ on a compact metric space
$(\mathcal{X},d)$ is called universal if and only if the RKHS $\mathcal{H}$
induced by the kernel is dense in $C(\mathcal{X})$, the space of continuous
functions on $\mathcal{X}$, with respect to the infinity norm
$\|f-g\|_{\infty}$.

######  Remark 3.

Existence of universal kernels allows one to approximate any continuous
function by some elements of RKHS. Steinwart [45] showed that the radial
kernels
$k(\boldsymbol{x},\boldsymbol{x}^{\prime})=\exp(-\lambda\|\boldsymbol{x}-\boldsymbol{x}^{\prime}\|^{2})$
is an example of universal kernel.

The following theorem provides a easy condition to assure independence of two
random variables, by using zero-covariance structure in RKHS reproduced by
universal kernels.

######  Theorem 10 (Gretton et al. [19]).

If $\mathcal{H}_{x}$ and $\mathcal{H}_{y}$ are RKHS generated from universal
kernels $k_{x},k_{y}$ on the compact metric spaces $\mathcal{X}$ and
$\mathcal{Y}$, then $\boldsymbol{\mathsf{x}}$ and $\boldsymbol{\mathsf{y}}$
are independent with respect to the joint measure $P_{x,y}$ on
$(\mathcal{X}\times\mathcal{Y},\Lambda_{x}\times\Lambda_{y})$ (where
$\Lambda_{x}$ and $\Lambda_{y}$ are Borel $\sigma$-algebras on $\mathcal{X}$
and $\mathcal{Y}$) if and only if the covariance operator $\Gamma_{xy}=0$.

######  Remark 4.

The theorem can be generalized to non-compact spaces also, without the choice
of universal kernel. However, the RKHS should then contain a sufficiently rich
class of functions to represent all higher order moments. The theorem for a
more lenient assumption about kernels has been proven by Fukumizu [14] via a
notion of characteristic kernels.

###  A.2 Vector Valued RKHS

As defined in the previous section, RKHS is created to structure the space of
functions from $\mathcal{X}$ to $\mathbb{R}$. A generalization to vector
valued function (in general operator valued function) is possible [5]. Let,
$\mathcal{X}$ be the domain space, and $\mathcal{W}$ be the range which is
also a Hilbert space. Let, $L(\mathcal{W})$ denotes the space of all bounded
linear operators on $\mathcal{W}$.

######  Definition 7.

A function $k:\mathcal{X}\times\mathcal{X}\rightarrow L(\mathcal{W})$ is said
to be an operator-valued positive definite kernel if

  1. 1.

For each pair $(x,y)\in\mathcal{X}\times\mathcal{X}$, $k(x,y)$ is a self-
adjoint operator in $L(\mathcal{W})$, i.e. for any
$w_{1},w_{2}\in\mathcal{W}$, $\langle
k(x,y)w_{1},w_{2}\rangle_{\mathcal{W}}=\langle
w_{1},k(x,y)w_{2}\rangle_{\mathcal{W}}$

  2. 2.

For any $x_{1},\dots x_{n}\in\mathcal{X}$ and any $w_{1},\dots
w_{n}\in\mathcal{W}$,

$\sum_{i=1}^{n}\sum_{j=1}^{n}\langle
w_{i},k(x_{i},x_{j})w_{j}\rangle_{\mathcal{W}}\geq 0$

Similar to defining the RKHS, we consider the basis functions given by
$k_{x}w(\cdot)=k(\cdot,x)w$ defined by $k_{x}w(y)=k(x,y)w$ for any
$x\in\mathcal{X},w\in\mathcal{W}$. The RKHS $\mathcal{H}$, given by the
closure of $\mathcal{H}_{0}$, where

$\mathcal{H}_{0}=\text{span}\left\\{k_{x}w(\cdot):x\in\mathcal{X},w\in\mathcal{W}\right\\}$

equipped with the inner product,

$\left\langle\sum_{i=1}^{\infty}\alpha_{i}k_{x_{i}}w_{i},\sum_{i=1}^{\infty}\beta_{i}k_{y_{i}}z_{i}\right\rangle_{\mathcal{H}}=\sum_{i=1}^{\infty}\sum_{j=1}^{\infty}\alpha_{i}\beta_{j}\langle
w_{i},k(x_{i},y_{j})z_{j}\rangle_{\mathcal{W}}$

The reproducing property here is that $\langle
f(x),w\rangle_{\mathcal{W}}=\langle f,k_{x}w\rangle_{\mathcal{H}}$ for any
$f\in\mathcal{H}$. In particular, if $\mathcal{W}=\mathbb{R}^{p}$ the
$p$-dimensional Euclidean vector space, then we can choose $w_{i}$’s as the
basis of the vector space, hence, the above infinite summands in expressing
the functions in the RKHS can be reduced to a finite sum over choices of
$w_{i}$, as we can simply take $w_{i}=\boldsymbol{e}_{i}$, for $i=1,2,\dots
p$, where $e_{i}$ is the $i$-th Euclidean basis vector. Therefore, the
functions in the RKHS can be expressed as,

$f(\cdot)=\sum_{i=1}^{\infty}\sum_{j=1}^{p}\alpha_{ij}k(x_{i},\cdot)\boldsymbol{e}_{j}$

###  A.3 Procrustes Analysis

Procrustes Analysis refers to a form of statistical analysis which is used to
analyze the shape of the distribution of datasets. The name “Procrustes”
refers to a bandit in greek mythology, who used to cut or stretch the limbs of
his victims to fit them into a bed.

There are many variations of Procrustes problem [17]. The simplest variation
asks one to find the optimal rotation matrix $\boldsymbol{P}$ such that the
matrix $\boldsymbol{A}$ can be transformed into $\boldsymbol{B}$ using the
orthogonal transformation $\boldsymbol{T}$. That means, the problem asks to
find
$\widehat{\boldsymbol{P}}=\arg\min_{\boldsymbol{P}}\|\boldsymbol{A}-\boldsymbol{PB}\|_{2}$,
where $\|\cdot\|_{2}$ denotes the usual Frobenius norm subject to
$\boldsymbol{P}^{\intercal}\boldsymbol{P}=\boldsymbol{I}$. Here, we consider a
general problem instead [43], that allows translation, uniform scaling and
optimal rotation to approximately transform a matrix $\boldsymbol{A}_{n\times
p}$ to $\boldsymbol{B}_{n\times p}$. Therefore, we consider the generalized
problem,

$\min_{\boldsymbol{\mu},\lambda,\boldsymbol{P}}\|\boldsymbol{A}-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{PB}\|_{2}^{2}$

where $\boldsymbol{\mu}$ is a $(p\times 1)$ column vector,
$\lambda\in\mathbb{R}$ is a scaling factor and $\boldsymbol{P}$ is an
orthogonal matrix. Hence, we consider a Lagrangian
$L=\|\boldsymbol{A}-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{PB}\|_{2}^{2}+\eta\|\boldsymbol{P}^{\intercal}\boldsymbol{P}-\boldsymbol{I}\|_{2}^{2}$,
where $\eta$ is the Lagrangian tuning parameter. Now, we compute the
derivatives of $L$ with respect to $\boldsymbol{\mu},\lambda$ and
$\boldsymbol{P}$ as follows.

$\displaystyle\dfrac{\partial L}{\partial\boldsymbol{\mu}}$
$\displaystyle=-2((\boldsymbol{A}-\lambda\boldsymbol{PB})^{\intercal}-\boldsymbol{\mu}\boldsymbol{1}_{n}^{\intercal})\boldsymbol{1}_{n}$
(5) $\displaystyle\dfrac{\partial L}{\partial\lambda}$
$\displaystyle=2\text{trace}\left(\boldsymbol{PB}((\boldsymbol{A}-\lambda\boldsymbol{PB})^{\intercal}-\boldsymbol{\mu}\boldsymbol{1}_{n}^{\intercal})\right)$
(6) $\displaystyle\dfrac{\partial L}{\partial\boldsymbol{P}}$
$\displaystyle=-2\lambda\left(A-\lambda\boldsymbol{PB}-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}\right)\boldsymbol{B}^{\intercal}+2\eta\boldsymbol{P}(\boldsymbol{P}^{\intercal}\boldsymbol{P}-\boldsymbol{I})$
(7)

Setting these derivatives given in Eq. (5)-(7) equal to zero to obtain the
normal equations, we obtain that,

$\displaystyle\widehat{\mu}$
$\displaystyle=\dfrac{1}{n}\left(\boldsymbol{A}-\widehat{\lambda}\boldsymbol{\widehat{P}B}\right)^{\intercal}\boldsymbol{1}_{n}$
(8) $\displaystyle\widehat{\lambda}$
$\displaystyle=\dfrac{\text{trace}\left(\widehat{\boldsymbol{P}}\boldsymbol{BA}^{\intercal}(\boldsymbol{I}-n^{-1}\boldsymbol{1}_{n}\boldsymbol{1}_{n}^{\intercal})\right)}{\text{trace}\left(\boldsymbol{BB}^{\intercal}(\boldsymbol{I}-n^{-1}\boldsymbol{1}_{n}\boldsymbol{1}_{n}^{\intercal})\right)}$
(9)

and $\boldsymbol{\widehat{P}}$ is the closest orthogonal matrix to
$\dfrac{1}{\widehat{\lambda}}\boldsymbol{\tilde{A}\tilde{B}}^{\intercal}$,
where $\tilde{\boldsymbol{A}}$ and $\tilde{\boldsymbol{B}}$ are demean
versions of $\boldsymbol{A}$ and $\boldsymbol{B}$ respectively. Thus, the
optimal rotation matrix is given by
$\widehat{\boldsymbol{P}}=\boldsymbol{UV}^{\intercal}$ such that
$\dfrac{1}{\widehat{\lambda}}\boldsymbol{\tilde{A}\tilde{B}}^{\intercal}$ has
a singular value decomposition (SVD) $\boldsymbol{U\Sigma V}^{\intercal}$.
Clearly, since this solution is independent of $\lambda$, one may obtain
$\widehat{\boldsymbol{P}}$ by simply performing SVD of
$\boldsymbol{\tilde{A}\tilde{B}}^{\intercal}$. The optimal value of
$\widehat{\lambda}$ and $\widehat{\boldsymbol{\mu}}$ can then be retrieved by
means of Eq. (8) and Eq. (9).

Therefore, the following steps calculates the optimal transformations from
$\boldsymbol{A}$ to $\boldsymbol{B}$.

  1. 1.

For each column in $\boldsymbol{A}$, compute its average and subtract it from
the corresponding entries. Similarly, demand the columns of $\boldsymbol{B}$.
Denote this demeaned version as $\tilde{\boldsymbol{A}}$ and
$\tilde{\boldsymbol{B}}$ respectively.

  2. 2.

Let $\boldsymbol{M}=\tilde{\boldsymbol{A}}\tilde{\boldsymbol{B}}^{\intercal}$,
and let its SVD is denoted as $\boldsymbol{M}=\boldsymbol{U\Sigma
V}^{\intercal}$. Optimal rotation is
$\widehat{\boldsymbol{P}}=\boldsymbol{UV}^{\intercal}$.

  3. 3.

Finally, the optimal scaling factor is simply,

$\widehat{\lambda}=\dfrac{\text{trace}\left(\widehat{\boldsymbol{P}}\tilde{\boldsymbol{B}}\tilde{\boldsymbol{A}}^{\intercal}\right)}{\text{trace}\left(\tilde{\boldsymbol{B}}\tilde{\boldsymbol{B}}^{\intercal}\right)}$

Note that, in particular,
$R=\|\boldsymbol{A}-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{PB}\|_{2}^{2}=\|\tilde{\boldsymbol{A}}-\lambda\boldsymbol{P\tilde{B}}\|_{2}^{2}$
is simply a quadratic function in $\lambda$, whose minimum is given by

$\displaystyle\min\|\tilde{\boldsymbol{A}}-\lambda\widehat{\boldsymbol{P}}\tilde{\boldsymbol{B}}\|_{2}^{2}$
$\displaystyle=\text{trace}(\tilde{\boldsymbol{A}}^{\intercal}\tilde{\boldsymbol{A}})-\dfrac{\text{trace}(\tilde{\boldsymbol{B}}^{\intercal}\widehat{\boldsymbol{P}}^{\intercal}\tilde{\boldsymbol{A}})^{2}}{\text{trace}(\tilde{\boldsymbol{B}}^{\intercal}\tilde{\boldsymbol{B}})}$
$\displaystyle=\text{trace}(\tilde{\boldsymbol{A}}^{\intercal}\tilde{\boldsymbol{A}})-\dfrac{\text{trace}(\tilde{\boldsymbol{A}}\tilde{\boldsymbol{B}}^{\intercal}\boldsymbol{U}^{\intercal}\boldsymbol{V})^{2}}{\text{trace}(B^{\intercal}B)}$
$\displaystyle=\text{trace}(\tilde{\boldsymbol{A}}^{\intercal}\tilde{\boldsymbol{A}})-\dfrac{\text{trace}(\boldsymbol{U\Sigma
V}^{\intercal}\boldsymbol{U}^{\intercal}\boldsymbol{V})^{2}}{\text{trace}(\tilde{\boldsymbol{B}}^{\intercal}\tilde{\boldsymbol{B}})}$
$\displaystyle=\text{trace}(\tilde{\boldsymbol{A}}^{\intercal}\tilde{\boldsymbol{A}})-\dfrac{\text{trace}(\boldsymbol{\Sigma})^{2}}{\text{trace}(\tilde{\boldsymbol{B}}^{\intercal}\tilde{\boldsymbol{B}})}$

where $\boldsymbol{\Sigma}$ is the diagonal matrix consisting of the singular
values of the matrix
$\tilde{\boldsymbol{A}}\tilde{\boldsymbol{B}}^{\intercal}$.

##  Appendix B Proofs

###  B.1 Proof of Theorem 1

Let us denote the errors as $\epsilon_{i}$ instead of
$\epsilon_{i}(\boldsymbol{\mathsf{x}})$ for $i=1,2$. Note that,
$\boldsymbol{\mathsf{x}}=f(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))+\epsilon_{1}=g(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))+\epsilon_{2}$
implies that,

$f(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))-g(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))=\epsilon_{2}-\epsilon_{1}$

Therefore, $(\epsilon_{2}-\epsilon_{1})$ is
$\mathcal{F}=\sigma(\\{\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}})\\})$-measurable.
Now,

$\displaystyle\epsilon_{2}-\epsilon_{1}$
$\displaystyle=\mathbb{E}(\epsilon_{2}-\epsilon_{1}\mid\mathcal{F}),\qquad\text{a.s.}$
$\displaystyle=\mathbb{E}(\epsilon_{2}\mid\mathcal{F})-\mathbb{E}(\epsilon_{1}\mid\mathcal{F})$
$\displaystyle=\mathbb{E}(\epsilon_{2})-\mathbb{E}(\epsilon_{1}),\qquad\text{a.s.}$
$\displaystyle\qquad\text{since, both }\epsilon_{1}\text{ and
}\epsilon_{2}\text{ are independent of }\mathcal{F}$ $\displaystyle=0$

Thus,
$f(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))=g(\phi_{1}(\boldsymbol{\mathsf{x}}),\dots\phi_{d}(\boldsymbol{\mathsf{x}}))$
almost surely.

###  B.2 Proof of Theorem 2

For any DR algorithm $\mathcal{A}_{n}$, the outputted embedding functions
$\widehat{\phi}_{1},\dots\widehat{\phi}_{d}$ are implicitly minimizers of the
loss function as described in Eq. (1). A representer theorem by Schölkopf,
Herbrich and Smola [40] can then be applied to represent these functions as a
finite linear combination of kernel functions, since the search space
$\mathcal{C}$ is a subset of some RKHS.

######  Theorem 11 (Schölkopf, Herbrich and Smola).

Suppose we are given a nonempty set $\mathcal{X}$, a positive definite real-
valued kernel $k$ on $\mathcal{X}\times\mathcal{X}$, a training sample
$(x_{1},y_{1}),\dots,(x_{m},y_{m})\in\mathcal{X}\times\mathbb{R}$, a strictly
monotonically increasing real-valued function $g$ on $[0,\infty]$, an
arbitrary cost function
$c:(\mathcal{X}\times\mathbb{R}^{2})^{m}\rightarrow\mathbb{R}\cup\\{\infty\\}$,
and a class of functions

$\mathcal{F}=\left\\{f\in\mathbb{R}^{\mathcal{X}}:f(\cdot)=\sum_{i=1}^{\infty}\beta_{i}k(\cdot,z_{i}),\beta_{i}\in\mathbb{R},z_{i}\in\mathcal{X},\|f\|<\infty\right\\}$

Then, for any $f\in\mathcal{F}$ minimizing the regularized risk functional

$c\left((x_{1},y_{1},f(x_{1})),\dots(x_{m},y_{m},f(x_{m}))\right)+g(\|f\|)$

admits a representation of the form

$f(\cdot)=\sum_{i=1}^{m}\alpha_{i}k(\cdot,x_{i})$

Now, since each of the functions $\phi_{j}$ belongs to Reproducible Kernel
Hilbert Space (RKHS) for some positive definite kernel $k$, one may transform
the constrained optimization problem for any algorithm into an unconstrained
problem by introducing Lagrangian multipliers as

$\min
L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\epsilon(\boldsymbol{\mathsf{x}}_{i})\\}_{i=1}^{n}\right)+\sum_{j=(d+1)}^{p}\eta_{j}\|\phi_{j}\|$

In view of a simple extension of Theorem 11 for more than one functions $f$,
we see that each $\phi_{j}$ can be represented as a finite linear combination
of kernel basis functions. Mathematically, this means, there exists constants
$\alpha_{ij}$ such that,

$\phi_{j}(\boldsymbol{\mathsf{x}})=\sum_{i=1}^{n}\alpha_{ij}k(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{x}}_{i}),\qquad
i=1,\dots n;j=1,\dots d$

Putting
$\mathsf{y}_{i^{\prime}j}=\phi_{j}(\boldsymbol{\mathsf{x}}_{i^{\prime}})$, we
obtain the following system of equations,

$\mathsf{y}_{i^{\prime}j}=\sum_{i=1}^{n}\alpha_{ij}k(\boldsymbol{\mathsf{x}}_{i^{\prime}},\boldsymbol{\mathsf{x}}_{i}),\qquad
i,i^{\prime}=1,\dots n;j=1,\dots d$ (10)

For fixed $j$, the system of equations given by Eq. (10) is a system of $n$
linear equations in $n$ variables. Also, the system is uniquely solvable since
the kernel matrix i.e. the matrix $\boldsymbol{K}$ with entries
$\\{k(\boldsymbol{\mathsf{x}}_{i},\boldsymbol{\mathsf{x}}_{i^{\prime}})\\}_{i,i^{\prime}=1}^{n,n}$
is a positive definite matrix. Therefore,

$\begin{pmatrix}\alpha_{1j}\\\ \alpha_{2j}\\\ \dots\\\
\alpha_{nj}\end{pmatrix}=\boldsymbol{K}^{-1}\begin{pmatrix}\mathsf{y}_{1j}\\\
\mathsf{y}_{2j}\\\ \dots\\\ \mathsf{y}_{nj}\end{pmatrix}$

Substituting these values to the form of $\phi_{j}(\boldsymbol{\mathsf{x}})$
yields the expression given in the theorem.

###  B.3 Proof of Theorem 3

We have that the $l$-th coordinate of the errors is given by
$\epsilon_{l}(\boldsymbol{\mathsf{x}})=\mathsf{x}_{\cdot
l}-f_{l}(\mathsf{y}_{\cdot 1},\dots\mathsf{y}_{\cdot d})$, where
$\mathsf{x}_{\cdot l}$ denotes the $l$-th variable in the dataset
$\boldsymbol{\mathsf{X}}$. Since we are concerned about the independence of
$\epsilon_{l}(\boldsymbol{\mathsf{x}})$ with the embeddings
$\phi_{j}(\boldsymbol{\mathsf{x}})$, which is usually hard to verify, the
assumption on RKHS actually allows one to simplify independence to only zero
covariance in RKHS domain. Since $\epsilon_{l}(\boldsymbol{\mathsf{x}})$ is
independent of $\phi_{j}(\boldsymbol{\mathsf{x}})$ if and only if
$\epsilon_{l}(\boldsymbol{\mathsf{x}})-c_{l}$ is independent of
$\phi_{j}(\boldsymbol{\mathsf{x}})$ for any choice of $c_{l}\in\mathbb{R}$,
one may assume without loss of generality that
$\epsilon_{l}(\boldsymbol{\mathsf{x}})$ has mean $0$.

In view of Theorem 10, to make sure that $\epsilon(\boldsymbol{\mathsf{x}})$
is independent of $\sigma\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}):j=1,2\dots
d\\}\right)$, it is enough to ensure that
$\text{Cov}\left(\epsilon_{l}(\boldsymbol{\mathsf{x}}),\phi_{j}(\boldsymbol{\mathsf{x}})\right)=0$
for all $l=1,2,\dots p;j=1,\dots d$. Since,
$\epsilon_{l}(\boldsymbol{\mathsf{x}})=\mathsf{x}_{\cdot
l}-f_{l}(\mathsf{y}_{\cdot 1},\dots\mathsf{y}_{\cdot d})$, we require,

$\text{Cov}(\boldsymbol{\mathsf{x}}_{\cdot
l},\phi_{j}(\boldsymbol{\mathsf{x}}))=\text{Cov}\left(f_{l}(\mathsf{y}_{\cdot
1},\dots\mathsf{y}_{\cdot d}),\phi_{j}(\boldsymbol{\mathsf{x}})\right),\qquad
l=1,2,\dots p;\ j=1,\dots d.$

Now, letting
$f_{l}(\cdot)=\sum_{i=1}^{n}\beta_{il}k_{y}(\cdot,\boldsymbol{\mathsf{y}}_{i})$
where $k_{y}$ is a kernel on $\mathbb{R}^{d}\times\mathbb{R}^{d}$, and
$\phi_{j}(\cdot)=\sum_{i=1}^{n}\alpha_{ij}k_{x}(\cdot,\boldsymbol{\mathsf{x}}_{i})$,
we obtain that,

$\text{Cov}\left(f_{l}(\boldsymbol{\mathsf{y}}),\phi_{j}(\boldsymbol{\mathsf{x}})\right)=\sum_{r=1}^{n}\beta_{rl}\sum_{s=1}^{n}({\Gamma}_{xy})_{rs}\alpha_{sj}$

where $\Gamma_{xy}$ is the covariance operator between RKHS $\mathcal{H}_{x}$
and $\mathcal{H}_{y}$. Therefore, we can now put the estimates of covariances
instead, to obtain the system of equations,

$\sum_{r=1}^{n}\beta_{rl}\left[\sum_{s=1}^{n}(\widehat{\Gamma}_{xy})_{rs}\alpha_{sj}\right]=\dfrac{1}{n}\sum_{i=1}^{n}\mathsf{x}_{il}\phi_{j}(\boldsymbol{\mathsf{x}}_{i})-\left(n^{-1}\sum_{i=1}^{n}\mathsf{x}_{il}\right)\left(n^{-1}\sum_{i=1}^{n}\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\right)$
(11)

Note that, in Eq. (11), the quantities $\alpha_{sj}$’s are known (or can be
estimated) using Theorem 2, and $\phi_{j}(\boldsymbol{\mathsf{x}}_{i})$’s are
also known for the training samples. For a fixed $l$, the system of linear
equations given in Eq. (11) contains $n$ unknown variables
$\beta_{1l},\dots\beta_{rl}$ in $d$ equations. Since generally, $d<<n$, the
system is underdetermined. Clearly, the system thus has more than one
solution, and it is clear that one should find the solution that minimizes the
number of nonzero $\beta_{rl}$’s, so that minimum number of parameters are
used to describe the reconstruction function. To ease the computation, we seek
for a solution to Eq. (11) such that the norm of
$\beta=(\beta_{1l},\dots\beta_{nl})^{\intercal}$ is minimized for each
$l=1,2,\dots n$.

Therefore, the problem boils down to finding $\boldsymbol{\beta}$ with minimum
$\boldsymbol{\beta}^{\intercal}\boldsymbol{\beta}$ such that
$\boldsymbol{A\beta}=\boldsymbol{c}$. Introducing the Lagrangian multiplier
$\lambda$, we consider the unconstrained minimization problem
$\boldsymbol{\beta}^{\intercal}\boldsymbol{\beta}+\lambda\left(\boldsymbol{A\beta}-\boldsymbol{c}\right)$.
Solving it yields the solution,
$\boldsymbol{\beta}=\boldsymbol{A}(\boldsymbol{A}^{\intercal}\boldsymbol{A})^{-1}\boldsymbol{c}$.
This leads us to the result stated in Theorem 3.

###  B.4 Proof of Theorem 4

Note that,

$\displaystyle\min_{\phi_{1},\dots\phi_{d}\in\mathcal{C}}L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right)$
$\displaystyle=\quad$
$\displaystyle\min_{\phi_{1},\dots\phi_{d}\in\mathcal{C}}L\left(\\{\phi_{j}(T^{-1}(\boldsymbol{\mathsf{z}}_{i}))\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right)$
$\displaystyle\geq\quad$
$\displaystyle\min_{\psi_{1},\dots\psi_{d}\in\mathcal{C}}L\left(\\{\psi_{j}(\boldsymbol{\mathsf{z}}_{i})\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right),\qquad\text{since,
}\phi_{j}\circ T^{-1}\in\mathcal{C}\text{ for any }\phi_{j}\in\mathcal{C}$

Conversely,

$\displaystyle\min_{\psi_{1},\dots\psi_{d}\in\mathcal{C}}L\left(\\{\psi_{j}(\boldsymbol{\mathsf{z}}_{i})\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right)$
$\displaystyle=\quad$
$\displaystyle\min_{\psi_{1},\dots\psi_{d}\in\mathcal{C}}L\left(\\{\phi_{j}(T(\boldsymbol{\mathsf{x}}_{i}))\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right)$
$\displaystyle\geq\quad$
$\displaystyle\min_{\phi_{1},\dots\phi_{d}\in\mathcal{C}}L\left(\\{\phi_{j}(\boldsymbol{\mathsf{x}}_{i})\\}_{i,j=1}^{n,p},\\{\boldsymbol{\mathsf{\epsilon}}_{i}\\}_{i=1}^{n}\right),\qquad\text{since,
}\psi_{j}\circ T\in\mathcal{C}\text{ for any }\psi_{j}\in\mathcal{C}$

This completes the proof.

###  B.5 Proof of Theorem 5

Since the set of functions in $\mathcal{H}^{\ast}$ is constrained by choice of
$\beta_{ij}$’s such that $0\leq\beta_{ij}\leq 1$ and their sum
$\sum_{i=1}^{\infty}\sum_{j=1}^{p}\beta_{ij}=1$, we may again introduce
Lagrangian multipliers to write the problem of obtaining Consistency Index as,

$\max_{T\in\mathcal{H}^{\ast}}\min_{\boldsymbol{\mu},\lambda,\boldsymbol{P}}\|\mathcal{A}_{n}(d,\tilde{\boldsymbol{\mathsf{X}}})(\tilde{\boldsymbol{\mathsf{X}}})-\boldsymbol{1}_{n}\boldsymbol{\mu}^{\intercal}-\lambda\boldsymbol{P}\mathcal{A}_{n}(d,\boldsymbol{\mathsf{X}})(\boldsymbol{\mathsf{X}})\|_{2}^{2}+\eta(\sum_{i,j}\beta_{ij}-1)+\sum_{i,j}\eta_{ij}\beta_{ij}$

Again, we can use a generalized version of representer theorem [33, 32] can be
used to claim that the maximizer $T$, if belongs to such a space
$\mathcal{H}^{\ast}$, can be represented as linear span of only kernel
functions evaluated at $n$ datapoints
$\boldsymbol{\mathsf{x}}_{1},\dots\boldsymbol{\mathsf{x}}_{n}$. This means,

$\widehat{T}(\boldsymbol{\mathsf{x}})=\sum_{i=1}^{n}\sum_{j=1}^{p}\beta_{ij}k(\boldsymbol{\mathsf{x}},\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j}$
(12)

where $k(\cdot,\cdot)$ is the reproducing kernel on
$\mathbb{R}^{p}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}\times\mathbb{R}^{p}$
and $\boldsymbol{e}_{j}$’s are the Euclidean basis vector of $\mathbb{R}^{p}$,
and $\widehat{T}$ is the argument of the maximizer in the definition of
Consistency Index. However, the number of such functions is uncountable, since
there are uncountably many choices for the real numbers $\beta_{ij}$’s.
Therefore, although the form of the maximizer function is known by Eq. (12),
the maximum value cannot be found by a computer so easily.

So we note that, $\partial\mathcal{H}^{\ast}\subseteq\mathcal{H}^{\ast}$,
hence clearly,
$\max_{T\in\mathcal{H}^{\ast}}S(T)\geq\max_{T\in\partial\mathcal{H}^{\ast}}S(T)$.

On the other hand,

$\displaystyle\max_{T\in\mathcal{H}^{\ast}}S(T)$
$\displaystyle\leq\max_{0\leq\beta_{ij}\leq
1,\sum_{i,j}\beta_{ij}=1}S\left(\sum_{i=1}^{n}\sum_{j=1}^{p}\beta_{ij}k(\cdot,\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j}\right),\qquad\text{this
is a bigger class}$ $\displaystyle\leq\max_{0\leq\beta_{ij}\leq
1,\sum_{i,j}\beta_{ij}=1}\left[\sum_{i=1}^{n}\sum_{j=1}^{p}\beta_{ij}S(k(\cdot,\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j})\right],\qquad\text{since
S(T) is convex}$ $\displaystyle\leq\max_{0\leq\beta_{ij}\leq
1,\sum_{i,j}\beta_{ij}=1}\left[\sum_{i=1}^{n}\sum_{j=1}^{p}\beta_{ij}\max_{T\in\partial\mathcal{H}^{\ast}}S(T)\right],\qquad\text{since,
}k(\cdot,\boldsymbol{\mathsf{x}}_{i})\boldsymbol{e}_{j}\in\partial\mathcal{H}^{\ast}$
$\displaystyle=\max_{T\in\partial\mathcal{H}^{\ast}}S(T)$

###  B.6 Proof of Theorem 6

Let, $A,B,C,D$ be 4 points such that $A,B\in C_{1}$ and $C,D\in C_{2}$ and
lying on the MCST. Figure 13 shows one such general configuration.

Figure 13: Four arbitrary points with two connections between two distinct
clusters in Euclidean MCST

Let us compute the conditional probability $P(E|A,B,C,D)$ where $E$ is the
required event that the MCST has two edges connecting two classes. Since the
MCST is spanning graph and is a tree, if we remove $AC$, then the graph
becomes disconnected. In other words, either $A$, $B$ are connected or $C$,
$D$ are connected. Otherwise, inclusion of $AC$ would lead to a cycle.

Without loss of generality, assume that $A$, $B$ are disconnected. Since the
tree is spanning, atleast one of the connected component of $A$ or $B$
(intersected with $C_{1}$) must be atleast as large as $(m-2)/2$. Say, the
connected component of $A$ is larger than $(m-2)/2$. Then, starting with the
initial MCST, removing $AC$ and connecting $A$ with any points from the
connected component of $B$ will again lead to a spanning tree. The former is a
minimum cost spanning tree, if for any such newly formed tree, the cost
increases, i.e. the new edge length is more than the length of $AC$.
Therefore, we have,

$\displaystyle P(E|A,B,C,D)$ $\displaystyle\leq
4P(\|A-X_{i}\|_{2}>|AC|,\forall X_{i}\in C_{1},X_{i}\rightsquigarrow B\mid
A,B,C,D)$ $\displaystyle\leq
4\prod_{i=1}^{(\min(m,n)-2)/2}P(\|A-X_{i}\|_{2}>|AC|\mid A,B,C,D)$
$\displaystyle\qquad\qquad\qquad\text{conditional on A},\|X_{i}-A\|_{2}\text{
are independent}$
$\displaystyle=4\alpha^{(\min(m,n)-2)/2}\qquad(\alpha<1,\text{ by
assumption})$ $\displaystyle\rightarrow 0$

Now, we can apply DCT to show that the unconditional probability also goes to
zero.

###  B.7 Proof of Theorem 8

Let us assume without loss of generality that
$\dfrac{1}{n}\sum_{i}\boldsymbol{x}_{i}=\dfrac{1}{n}\sum_{i}\boldsymbol{y}_{i}=0$.
Let, $\boldsymbol{X}$ and $\boldsymbol{Y}$ be the $n\times p$ matrices with
each row being $\boldsymbol{x}_{i}$ and $\boldsymbol{y}_{i}$ respectively.
Then, the usual double centering mechanism (see Borg and Groenen [4] for
details) will indicate that

$\boldsymbol{XX}^{\intercal}=-\dfrac{1}{2}\boldsymbol{CQC}=\boldsymbol{YY}^{\intercal}$

where $\boldsymbol{C}=\boldsymbol{I}_{n}-\dfrac{1}{n}\boldsymbol{J}_{n}$ with
$\boldsymbol{J}_{n}$ being the $n\times n$ matrix with all entries being equal
to $1$, and $\boldsymbol{Q}$ is the matrix with the entries
$q_{ij}^{2}=\|\boldsymbol{x}_{i}-\boldsymbol{x}_{j}\|^{2}$.

Now, a simultenous singular value decomposition of $\boldsymbol{X}$ and
$\boldsymbol{Y}$ would yield that

$\boldsymbol{X}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}_{1}^{\intercal},\qquad\boldsymbol{Y}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}_{2}^{\intercal}$

where the same $\boldsymbol{U},\boldsymbol{V}_{1},\boldsymbol{V}_{2}$ are
unitary matrices and $\boldsymbol{D}$ is a diagonal matrix. The same
$\boldsymbol{U}$ and $\boldsymbol{D}$ works since
$\boldsymbol{XX}^{\intercal}=\boldsymbol{YY}^{\intercal}$. This shows that
$\boldsymbol{X}=\boldsymbol{Y}\boldsymbol{V}_{2}\boldsymbol{V}_{1}^{\intercal}$,
i.e. the matrix
$\boldsymbol{P}=\boldsymbol{V}_{2}\boldsymbol{V}_{1}^{\intercal}$ works as the
required orthogonal transformation.

###  B.8 Proof of Theorem 9

Since $g$ is a geodesic of $\mathcal{M}$ and a unit-velocity curve in
$\mathbb{R}^{d}$ induces a constant velocity curve in $\mathcal{M}$, it
follows that (see Lee [26] for details)

$g(f(\boldsymbol{\theta}_{i}),f(\boldsymbol{\theta}_{j}))=c_{1}\|\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{j}\|,\qquad\forall
1\leq i<j\leq n.$

Therefore, $\\{\boldsymbol{y}_{1},\dots\boldsymbol{y}_{n}\\}$ is the output of
metric MDS for the distance metric
$q_{ij}=c_{1}\|\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{j}\|$. On the
other hand, due to the definition of $g^{\prime}$, the geodesic on
$\mathcal{M}^{\prime}$, it follows that

$g^{\prime}(\phi\circ f(\boldsymbol{\theta}_{i}),\phi\circ
f(\boldsymbol{\theta}_{j}))=cc_{1}\|\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{j}\|,\qquad\forall
1\leq i<j\leq n$

i.e. $\\{\boldsymbol{z}_{1},\dots\boldsymbol{z}_{n}\\}$ is the output of
metric MDS for the distance metric
$q^{\prime}_{ij}=cc_{1}\|\boldsymbol{\theta}_{i}-\boldsymbol{\theta}_{j}\|=cq_{ij}$.
Now, proceeding along the lines of the proof of Theorem 8, one can show the
required result.

##  Appendix C More Simulated Examples

Several simulated examples, their embeddings using tSNE, UMAP and LSDR, under
different hyperparameters such as perplexities, n-neighbours, kernel
bandwidths etc. are described in Figure 14-24. These examples are taken from
the analysis of Google PAIR engineers [7, 54], where an extensive comparison
of tSNE and UMAP is already available. These are some of the datasets where
tSNE and UMAP are found to have better performances, and yet in these
datasets, LSDR achieves a similar level of performances compared to them. In
particular, in linked circle data (see Figure 18) and in trefoil knot data
(see Figure 24), both tSNE and UMAP only sees the local structure, and hence
outputs two unlinked circles and a large circle respectively. However, there
does not exist a continuous transformation that can map those two linked
circles into two unlinked circles, and similarly, the trefoil knot, being a
non-trivial knot, is not isomorphic to a single large circle. These notions
can only be achieved if we look at a global scale, and LSDR takes care of
that.

On the other hand, in the two long linear cluster data (see Figure 18), it is
possible to separate the two clusters by a single line, hence even a one-
dimensional embedding should be able to perform a reduction of higher quality.
However, tSNE and UMAP does not necessarily show this importance in their
first reduced component. In other words, it is not necessarily true that the
first reduced variable obtained by tSNE or UMAP is necessarily better than its
second reduced component. But, LSDR turns out to be capable of using such
structure in its first reduced output only, while the secondary component (or
reduced variable) provides an orthogonal direction to look at. This phenomenon
can be explained by the fact that the metric MDS step in LSDR (which provides
the low dimensional embeddings for the skeletal points) actually perform an
eigendecomposition, which creates several orthogonal eigen-directions in the
order of the most explainable variation of distances through the manifold
approximation graph $G_{\alpha}$.

(a) Original Grid Dataset

(b) 2D embedding of Grid data by tSNE for different perplexities

(c) 2D embedding of Grid data by UMAP for different n-neighbours

(d) 2D embedding of Grid data by tSNE for different kernel bandwidth

Figure 14: Analysis of Grid Data

(a) Original data with two long linear clusters

(b) 2D Embedding of Two linear clusters data via tSNE for different
perplexities

(c) 2D Embedding of Two linear clusters data via UMAP for different
n-neighbours

(a) 2D Embedding of Two linear clusters data via LSDR for different
n-neighbours $k$ and different kernel bandwidth combination Figure 16:
Analysis of two long linear clusters data

(a) Original data with two linked circle in 3D

(b) 2D Embedding of Two linked circles data via tSNE for different
perplexities

(c) 2D Embedding of Two linked circles data via UMAP for different
n-neighbours

(a) 2D Embedding of Two linked circles data via LSDR for different
n-neighbours $k$ and different kernel bandwidth combination Figure 18:
Analysis of two linked circles data

(a) Original data with two unlinked circle in 3D

(b) 2D Embedding of Two unlinked circles data via tSNE for different
perplexities

(c) 2D Embedding of Two unlinked circles data via UMAP for different
n-neighbours

(a) 2D Embedding of Two unlinked circles data via LSDR for different choice of
$\alpha$-approximation and different kernel bandwidth combination, fixing
$k=5$ Figure 20: Analysis of two unlinked circles data

(a) Original data with circular clusters in 3D

(b) 2D Embedding of circular clusters data via tSNE for different perplexities

(c) 2D Embedding of circular clusters data via UMAP for different n-neighbours

(a) 2D Embedding of circular clusters data via LSDR for different choice of
n-neighbours $k$ and different kernel bandwidth combination Figure 22:
Analysis of circular clusters data

(a) Original trefoil knot data in 3D

(b) 2D Embedding of trefoil knot data via tSNE for different perplexities

(c) 2D Embedding of trefoil knot data via UMAP for different n-neighbours

(a) 2D Embedding of trefoil knot data via LSDR for different choice of
n-neighbours $k$ and different kernel bandwidth combination Figure 24:
Analysis of trefoil knot data

## References

  * [1] Mikhail Belkin and Partha Niyogi.  Semi-supervised learning on riemannian manifolds.  Machine learning, 56(1):209–239, 2004.
  * [2] Alain Berlinet and Christine Thomas-Agnan.  Reproducing kernel Hilbert spaces in probability and statistics.  Springer Science & Business Media, 2011.
  * [3] Harry Blum et al.  A transformation for extracting new descriptors of shape, volume 4.  MIT press Cambridge, 1967.
  * [4] Ingwer Borg and Patrick JF Groenen.  Modern multidimensional scaling: Theory and applications.  Springer Science & Business Media, 2005.
  * [5] Claudio Carmeli, Ernesto De Vito, Alessandro Toigo, and Veronica Umanitá.  Vector valued reproducing kernel hilbert spaces and universality.  Analysis and Applications, 8(01):19–61, 2010.
  * [6] Stephan K Chalup and Andreas Mitschele.  Kernel methods in finance.  In Handbook on information technology in finance, pages 655–687. Springer, 2008.
  * [7] Andy Coenen and Adam Pearce.  Understanding umap, 2019.
  * [8] Carlos Correa and Peter Lindstrom.  Towards robust topology of sparsely sampled data.  IEEE Transactions on Visualization and Computer Graphics, 17(12):1852–1861, 2011.
  * [9] M. de Berg, O. Cheong, M. van Kreveld, and M. Overmars.  Computational Geometry: Algorithms and Applications.  Springer Berlin Heidelberg, 2008.
  * [10] Boris Delaunay et al.  Sur la sphere vide.  Izv. Akad. Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk, 7(793-800):1–2, 1934.
  * [11] R Dennis Cook.  Save: a method for dimension reduction and graphics in regression.  Communications in statistics-Theory and methods, 29(9-10):2109–2121, 2000.
  * [12] Dheeru Dua and Casey Graff.  UCI machine learning repository, 2017.
  * [13] Greg Friedman.  An elementary illustrated introduction to simplicial sets.  arXiv preprint arXiv:0809.4221, 2008.
  * [14] Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Schölkopf.  Kernel measures of conditional dependence.  In NIPS, volume 20, pages 489–496, 2007.
  * [15] Keinosuke Fukunaga.  Introduction to statistical pattern recognition.  Elsevier, 2013.
  * [16] R Paul Gorman and Terrence J Sejnowski.  Analysis of hidden units in a layered network trained to classify sonar targets.  Neural networks, 1(1):75–89, 1988.
  * [17] John C Gower, Garmt B Dijksterhuis, et al.  Procrustes problems, volume 30.  Oxford University Press on Demand, 2004.
  * [18] John Greene.  Feature subset selection using thornton’s separability index and its applicability to a number of sparse proximity-based classifiers.  In Proceedings of Annual Symposium of the Pattern Recognition Association of South Africa, 2001.
  * [19] Arthur Gretton, Ralf Herbrich, Alexander Smola, Olivier Bousquet, Bernhard Schölkopf, et al.  Kernel methods for measuring independence.  2005\.
  * [20] Amy F Iezzoni and Marvin P Pritts.  Applications of principal component analysis to horticultural research.  HortScience, 26(4):334–338, 1991.
  * [21] Jagadeesh Jagarlamudi and Raghavendra Udupa.  Revisiting dimensionality reduction techniques for nlp.
  * [22] N. Jayant.  Signal Compression: Coding of Speech, Audio, Text, Image and Video.  Selected topics in electronics and systems. World Scientific, 1997.
  * [23] John NR Jeffers.  Two case studies in the application of principal component analysis.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 16(3):225–236, 1967.
  * [24] Nandakishore Kambhatla and Todd K Leen.  Dimension reduction by local principal component analysis.  Neural computation, 9(7):1493–1516, 1997.