\section{INTRODUCTION}
\cite{price1}
reported that citations networks had a heavy-tailed distribution following a power law.
He also described a ``cumulative advantage'' that explains the power law~\citep{price2}.
This same networks with a power-law distribution of connectivity is nowadays called ``scale-free''
and the same explanation for the distribution is called ``preferential attachment'',
as the result of a rediscovery of the property by~\cite{barabasi1}.
This power-law rediscovery, and the rediscovery of small-world networks by~\cite{small},
is usually considered the birth of the complex networks field,
in the turn of the millennium.
In the last decade and a half,
most of the articles in physics report advances in complex networks.

Another interesting fact about the meaning of the field
is related to the essence of these structures.
Components (be them human agents or not), in diverse contexts, often
exhibit a power-law distribution of activity.
If their activity yields links,
the power-law distribution of connectivity is one of the byproducts of activity.
Indeed, connectivity and activity present high correlation in many of such scenarios~\citep{fabbri1,fabbri2,fabbri3}.

The complex network literature stresses concentration and often revolves around the hub.
More accurately, it glamorizes the most connected nodes as
``both the strength and the weakness of scale-free networks''~\citep{wikipedia,networks},
and places them as the most important vertexes.
In this work, we argue that this is naive for a number of reasons.
Therefore, this document is countercurrent with respect to available literature:
not only the focus here is to observe scale-free networks egalitarian aspects,
but also to dilute the hubs boast.
No academic writing was found by the authors to expose this simple and pertinent content.

Section~\ref{sec:can} presents a glimpse at standard theory related to power laws.
Section~\ref{sec:three} states the aspects of scale-free networks that are egalitarian.
Section~\ref{sec:delusion} considers also the hubs, intermediary and peripheral vertices
and a probable unbalance current literature has in considering the importance of these sectors.
Section~\ref{sec:con} is dedicated to conclusions and future work.

\section{CANONICAL BACKGROUND}\label{sec:can}
Networks with a power-law distribution of connectivity (degree or strength\footnote{Degree of a vertex
is the number of links connected to it.
Strength of a vertex is the sum of the weight of the links connected to it.})
are called \emph{scale-free}.
In other words, be $p(k)$ the probability that an arbitrary vertex has degree $k$,
than, for a scale-free network, one can assume:
\begin{equation}
p(k) \sim k^{-\gamma}
\end{equation}
\noindent where $\gamma > 1$ is constant and typically $\gamma \in [2,3]$.
This same distribution is called, under certain conditions, e.g. the Pareto distribution or the Zipf's law.
In fact, they all report that nature is performing a distribution of resources in
a way that can be modeled by power laws.
The next sections exposes this content very briefly.

\subsection{Zipf's law}
This law is unique among the ones mentioned in this article because it
relates the frequency of an element to its rank in an ordering where the
most frequent comes first.
More specifically, a word with rank $r$ has the frequency (or probability of occurring):
\begin{equation}\label{eq:pl}
p(r) \sim r^{-\gamma}
\end{equation}
\noindent where $\gamma \approx 1$.
If the definition given above is followed strictly, the law entails
that:

\begin{equation}
	p(r) = \frac{r^{-\gamma}}{\sum_{n=1}^N n^{-\gamma}}
\end{equation}
\noindent where N is the vocabulary size.
Mandelbrot proposed a shifted version of this law that
follows more closely the distributions found in real languages:

\begin{equation}
	p(r) \sim (r+s)^{-\gamma}
\end{equation}
\noindent where $s\approx 2.7$.

There are many non-trivial considerations to be made
about the Zipf's law but, overall, it seems to hold universally
in all languages, including those that are extinct and untranslated~\citep{zipf}.

\subsection{Pareto distribution}
The Pareto distribution is defined in many ways throughout the literature.
\cite{paretoWikipedia}, for example, defines the distribution through
its survival function and presents 4 variants and a generalized form
with an extra parameter.
In its most basic form, the Pareto distribution can be regarded as
having the probability density function:
\begin{equation}
	p(k) = \frac{\gamma k_m^\gamma}{k^{\gamma+1}}
\end{equation}
\noindent where $k_m$ is the minimum value for $k$ to yield $p(k)>0$.
This distribution is used to describe many observable phenomena,
including in social, financial and geophysical theories.

\subsection{Human perception: Steven's law}
The relation between the magnitude of a physical stimulus and
its perceived intensity is often given by the Steven's law:
\begin{equation}
	S(I) = kI^\gamma
\end{equation}
\noindent where $S(I)$ is the perceived intensity (or strength)
of intensity $I$ of stimulus type that entails $\gamma$ and
$k$ is a proportionality constant that accounts for the arbitrary
units used.

\subsection{Many other laws}
Notice that this last law (Steven's) is not a probability law as in the other cases above.
In fact, many other phenomena are modeled as power laws, especially natural phenomena,
such as Newton's laws of gravitation and inertia.
This calls for a consideration of the expressive capabilities of power laws and
for an honest inquiry of what is so special about these functions that
can account for its ubiquity~\citep{fabbri3}.

\subsection{Scale-free networks}
Now back to the main subject of this article.
A scale-free network is one where the distribution of connectivity
follows the Equation~\ref{eq:pl}.
These networks also present other non-trivial properties,
such as the 'small-world-ness', high clusterization and
abundant motifs~\citep{networks}.
Most interesting for this discussion, these networks are not isolated
entities, are not static and the power-law distribution is directly related to
an uniform distribution as will be exposed in the next section.

\section{THREE EGALITARIAN ASPECTS}\label{sec:three}
This is the most important section of this article.
We strived to keep everything as simple as possible to allow
both this document to be reasonable to a large audience
and for the reader to promptly grasp the contents.
The kernel of our subject is: given the ubiquity of power-laws,
and that they arise from an uniform distribution,
and that, in the context of social networks, these structures
are not isolated nor static,
what should be the emphasis given to the association of scale-free
networks to inequality?
In the same context: does it explain the ubiquity of inequality
but states, at the same time, paradoxically, that equality is also prevalent?
We should examine each of the egalitarian aspects of scale-free networks:
\begin{itemize}
	\item according to the exposition in~\cite{fabbri3},
		a power-law is entailed by a uniform distribution.
		The same quantity $C$ of resources
		is distributed uniformly across all concentrations of resources:
		$\forall k : p_U(k)=C \Rightarrow p(k) = C.k^{-\gamma}$
		where $\gamma$ is the dimensionality of the resource.
                In other words, the uniform distribution of a $\gamma$-dimensional
		resource result is a power law with coefficient $\gamma$.
		The more resources is allocated per component (e.g. per agent), the less numbered
		are the components.
	\item The networks are not static!
		Some of them might be (e.g. those yield by bone cavities)
		but social networks, which are at the center of our discussion, are not.
		Social networks are in constant change and, in fact, stable hubs were only observed
		in very small networks~\citep{barabasi}.
		Although the overall distribution is somewhat stable (or invariant) in the power-law outline,
		the participants move constantly between the periphery to intermediary and hubs sectors~\citep{fabbri1,fabbri2,versinus}.
		Thus, the scale-free (social) networks are egalitarian in the sense that
		the participants present only transient concentrations of resources.
	\item The networks are not isolated!
		This has deep consequences for social networks.
		For example, you, the reader, are a hub in some nucleus of your own family.
		Some friend of yours is a peripheral in your family because he/she only knows you.
		In the same way, everyone is a hub, a peripheral and an intermediary\footnote{A
		sound method for deriving the hubs, intermediary and peripheral sectors,
		by comparing the real network against the Erdös-Rényi model, is described
		in~\cite{fabbri1} and~\cite{fabbri2}.} at the same time,
		but in a (vast) number of different networks.
\end{itemize}

\section{EXALTATION OF HUBS AND DELUSIONS OF GRANDEUR}\label{sec:delusion}
The hubs, defined as the nodes with the greatest number of links in a scale-free network,
are reported in many studies as being the most relevant components.
In fact, selective attacks to the hubs (generally modeled as removing them)
cause the networks to loose many of the properties that make scale-free networks
appealing: the low average-distance between nodes, the fast transmission of information
and other sorts of transportation~\citep{networks}.
Nevertheless, aging of hubs is a recognized phenomenon,
i.e. the descent of the most connected hubs is observed
and even expected~\citep{aging}.
Also, diseases spread rapidly by means of the hubs
and there are many other important negative contributions
of the hubs and important positive contributions of the
other vertices.
In social networks terms, often used for such discussions:
\begin{itemize}
    \item Hubs are but a few vertexes, which are replaced constantly.
    \item Hubs present the most trivial behavior:
	    they interact as much as possible, in every situation, with every other agent.
		That is a prerequisite for they to be hubs.
    \item Hubs tend to present corrupt behavior,
	    as they apply huge amounts of time and energy (resources)
		on the network where it is a hub and frequently depend on such activity for basic provision.
		This might also be seen as a tendency of the hub to have a committed attitude,
		but given that hubs are also ubiquitous in non-professional networks,
		this might be understood as corrupt behavior by the non-warned.
    \item World input to the network is done mostly by peripheral and intermediary agents,
	    as they do not apply all their energy on the network.
    \item Authorities are often intermediaries or less active hubs,
	    especially in cases where they deliver quality, not quantity.
		This can be noticed by the fact that authorities usually exploit their authority
		in order to achieve their goals by means of the dedication of other agents.
	\item Network structure is yield by intermediary vertices~\citep{fabbri2},
	    as they are the only with non-trivial behavior:
		peripheral vertexes interact only a few times; hubs interact everywhere possible.
\end{itemize}

In summary, the constant exaltation of the hubs in current
complex networks literature seems to be related to the
fact that the field is still incipient.
A psychologist might also add that it is related to
grandiose delusions which affect a great fraction of
our society and thus leaks into the scientific tradition~\citep{grandiose}.
A social scientist might add that there is a constant
pressure performed by the more powerful and rich groups
to legitimate themselves and the inequality as a natural tendency~\citep{socialB,marcuse}.
An anthropologist might mention societies where the inequality
is symbolized in rituals as something to be avoided
and extreme cases in which the existence of leaders is denied~\citep{clastres}.
It is also appropriate to remember the notes of~\cite{heisenberg}
about nature, physics and human social power.
Anyway, there is an evident consonance between the
constant ennoblement of the hubs, the stress on the association
of scale-free networks and inequality,
and our observation that these associations should be
considered with parsimony.

\section{CONCLUSIONS AND FURTHER WORK}\label{sec:con}
The observations made in this brief article have direct consequences for
the purposes of the mathematical and agent-based models of complex networks.
How can our models reflect the diversity input by peripherals
and the greater authority held by the intermediary sector or the less connected hubs?
Considering our current scientific literature:
should one build a survey of the
excesses about inequality and exaltation of hubs in recent complex networks books, article and software?
More fundamentally, are the issues reported in this article sound
and relevant?
Given that the complex networks field is recent and that it is as driven
by social constructs and conveniences as any knowledge field,
we need to choose a very lucid path
in 1) asking more questions than stating fresh conclusions as facts;
2) establishing a deeper dialog with the traditions that have studied
these phenomena for centuries in advance, such as social scientists and philosophers.

Further work should involve attempts to answer the questions and follow the recommendations above,
further confirm or refute the arguments in this documents, and:
\begin{itemize}
	\item develop models to reflect the fundamental roles performed by nodes beyond hubs (i.e.
		peripherals and intermediary).
	\item Deepen the model exposed in~\cite{fabbri3} and submit it to peer review processes in order
		to understand the extent of its validity.
	\item Study the dynamics of individual agents in social networks with respect to the concentration of links,
		as suggested in~\cite{fabbri1,fabbri2,versinus}.
	\item Study the integration of the networks, maybe by using the database in~\cite{losd} which
		contains a number of networks with shared participants.
	\item Deepen the absorption of the literature and the partnership with academics which produce
		knowledge within the complex networks field, both to better understand the convictions
		and to spread our concerns.
	\item Consider an antithesis: the formalism of institutionalized posts
		might impose both a less dynamic variation of participant roles
		and a less extreme concentration of resources.
\end{itemize}

\subsection*{\textit{Acknowledgements}}
The authors thank CNPq for the funding received while researching the topic of this article,
the researchers of IFSC/USP and ICMC/USP for the recurrent collaboration in every situation
where we needed directions for investigation.


\section{Introduction}

With recent advances in data collection and storage capabilities, researchers across a wide variety of fields are facing larger datasets with increasing dimensionality, arising from images, videos, text documents, etc. Dimensionality reduction methods have been developed and applied as pre-processing tools to deal with such high-dimensional datasets.

In dimensionality reduction, we assume that the observed data in high-dimensional space lie on (or near) an embedded submanifold with lower dimensionality. With this fundamental assumption, it is often possible to represent the high-dimensional data in a low-dimensional space without losing important information from the original data.

Unless stated otherwise, in this paper all vectors are assumed to be column vectors. For a vector $\mathbf{x}$, its transpose is denoted by $\mathbf{x}^{\prime}$. Formally, we state the dimension reduction problem as follows. Suppose that there are $n$ points in a $q$-dimensional space $\mathbb{R}^{q}$, denoted individually by $\left\{ \mathbf{y}_{1},\ldots,\mathbf{y}_{n} \right\}$, or collectively by an $n\times q$ matrix $ \mathbf{Y} $ with the $j$-th row being $\mathbf{y}_{j}^{\prime}$. Furthermore, we assume $\mathbf{y}_{1},\ldots,\mathbf{y}_{n}$ are mapped into a higher-dimensional space $\mathbb{R}^{p}$ by an unknown isometry $\varphi:\mathbb{R}^{q} \rightarrow \mathcal{M}^{q} \subset \mathbb{R}^{p}$ ($p>q$) possibly with added random error:

\begin{equation}
\label{eq1}
\mathbf{x}_{j}=\varphi(\mathbf{y}_{j})+\pmb{\epsilon}_{j}\,,
\end{equation}
where $\pmb{\epsilon}_{j}\in \mathbb{R}^{p},j=1,\ldots,n$ are i.i.d. with mean 0, and $\mathcal{M}^{q}$ is a differentiable embedded submanifold.

We represent $\left\{ \mathbf{x}_{1},\ldots,\mathbf{x}_{n} \right\}$ collectively as an $n\times p$ matrix $ \mathbf{X} $. The vectors $\left\{ \mathbf{y}_{1},\ldots,\mathbf{y}_{n} \right\}$ are not observed, and we only observe $\left\{ \mathbf{x}_{1},\ldots,\mathbf{x}_{n} \right\}$. We say $\left\{ \mathbf{x}_{1},\ldots,\mathbf{x}_{n} \right\}$ lie on (or near) the manifold $\mathcal{M}^{q}$ with \textsl{intrinsic dimensionality} $q$, or we say the intrinsic dimensionality of $\mathbf{X}$ is $q$. By ``case $j$'' we shall mean the index of the corresponding points $\mathbf{x}_{j}$ and $\mathbf{y}_{j}$.

The purpose of dimensionality reduction algorithms is to reconstruct the inverse mapping $\psi=\varphi^{-1}$ up to an arbitrary isometry, and to reconstruct $\mathbf{y}$ by $ \widehat{\mathbf{y}}=\widehat{\psi}(\mathbf{x}) $, again up to an arbitrary isometry. Some dimension reduction methods do not reconstruct $\psi$ directly, but just $\mathbf{y}$. We denote a given dimensionality reduction method as a mapping $\widehat{\psi}:\mathbb{R}^{p}\rightarrow \mathbb{R}^{q}$, and we also write the low-dimensional configuration $\widehat{\mathbf{Y}}$ as $\widehat{\mathbf{Y}}=\widehat{\psi}(\mathbf{X})$ for convenience.

Dimension reduction algorithms can be broadly classified into two types: linear and nonlinear, depending on whether $\widehat{\psi}$ is a linear mapping or not. When $\mathcal{M}$ is assumed to be a $q$-flat in $\mathbb{R}^{p}$, it is natural to use a linear method. Well-known linear techniques include principal component analysis \citep{pca1} and metric multidimensional scaling \citep{mds1}. Nonlinear dimension reduction methods cover all other cases, and they will typically be used when $\mathcal{M}$ is assumed to have some extrinsic curvature, which means $\mathcal{M}$ is not a flat. Over the last two decades, many nonlinear methods have been introduced in the literature. These include Kernel PCA \citep{kpca}, Local Linear Embedding (LLE) \citep{lle}, Isometric Feature Mapping (Isomap) \citep{isomap}, Laplacian Eigenmap \citep{LEM}, Maximum Variance Unfolding (MVU) \citep{mvu} and Local Tangent Space Alignment (LTSA) \citep{ltsa}, among others. Summaries
and surveys can be found in many books and papers \citep{sum1,sum2,sum3}.

There is no universal agreement on how to assess and compare the performance of these different methods. In the supervised learning problems, such as regression or classification, a performance measure can be defined on the response variable (such as prediction error or classification error). But dimensionality reduction, as we formulated here, is an unsupervised learning problem without such a measure. Therefore, a performance measure is needed to quantify the quality of different dimension reduction algorithms. Ideally, we expect this
\begin{itemize}
\item to be easily interpretable,

\item to be applicable to most algorithms and datasets,

\item to have the property that any tuning parameters of the measure can be easily selected by a simple method.
\end{itemize}
In addition to these broad criteria, other considerations are necessary.  For example, data which arise in high dimensional contexts will often have interpoint distance distributions which are highly skewed, even locally for nearest neighbor distances.  This is part of the so-called ``curse of dimensionality.''  When this occurs, standard multivariate methods, such as arise from Gaussian assumptions used in low dimensions, are not suitable.  Distances with heavily skewed distributions can be transformed to more symmetric variables on an ad hoc basis.  An alternative, which we choose here, is to use ranking methods for nearest neighbor distances.  Therefore we can add the following requirement.
\begin{itemize}
\item We expect this method to avoid the problems associated with the curse of dimensionality.
\end{itemize}
In this paper, we develop a local rank-correlation measure to accomplish these goals.

In equation \eqref{eq1}, we have assumed that $\varphi$ is an isometry. It follows that the neighboring points in the input space should be mapped to neighbors in the output space, and vice versa for the inverse mapping $\psi$. This can be called ``neigborhood preservation.''
Early attempts to quantify the neigborhood preservation of a dimension reduction method were made in the study of Self-Organizing Maps \citep{SOM}. where measures such as the \textit{topographic product} \citep{qa1}, \textit{topographic function} \citep{qa2} and \textit{quantization error} \citep{qa6} were developed. Advantages and disadvantages of each method are discussed in detail by P\"{o}lzlbauer \citep{qa5}.

More recently, a few rank-based measures have been proposed. These include \textit{mean relative rank errors} (MRREs) \citep{rank2}, \textit{trustworthiness and continuity} (T\&C) \citep{rank1}, \textit{local continuity meta criterion} (LCMC) \citep{rank5}, and the \textit{agreement rate metric} (AR) \citep{rank4}.
Measures MRREs and T\&C are restricted to the interval $[0,\,1]$ and order outcomes of algorithms in such a way that larger values of these measures are desirable. Both of these measures try to quantify two types of errors that occur during the dimension reduction procedures,
\begin{itemize}
\item[(i)] non-neighboring points in $\mathbb{R}^{p}$ are mapped by $\widehat{\psi}$ to be neighboring points in $\mathbb{R}^{q}$,

\item[(ii)] neighboring points in $\mathbb{R}^{p}$ are mapped by $\widehat{\psi}$ to be non-neighboring points in $\mathbb{R}^{q}$.
\end{itemize}

These two types of errors create a discrepancy between nearest neighbor ranks in the input and output spaces. Therefore they can be measured by calculating the change of nearest neighbor ranks.
Measure ${\rm AR}_{_J}$ is the average size of the overlap of $J$-nearest neighborhoods between the low-dimensional reconstruction and the original data, while ${\rm LCMC}_{_J}$ accounts for the expected random overlap. Besides the neighborhood preservation measures, \cite{PA} proposed a Procrustes measure that evaluates how well each local neighborhood matches its corresponding embedding under an optimal linear transformation.

The rest of the paper is organized as follows. In Section 2, we introduce a new class of local rank correlation measures. In Section 3, we provide some applications on benchmark datasets. Finally, in Sections 4 and 5, we employ local rank correlations to choose suitable values of parameters used in modelling and tuning.

\section{Local rank correlation}

For an observed high-dimensional dataset $ \left\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n} \right\} \subset \mathcal{M}$ and a low-dimensional representation $\left\{ \widehat{\mathbf{y}}_{1},\ldots,\widehat{\mathbf{y}}_{n} \right\}$, we have the following notation:

\begin{itemize}
\item $\left\| \cdot \right\|$: the Euclidean norm.

$d^{\mathcal{M}}(\mathbf{x}_1,\,\mathbf{x}_2)$: the geodesic distance from $\mathbf{x}_1$ to $\mathbf{x}_2$ on the Riemannian manifold $\mathcal{M}$.

\bigskip

\item $s_{_{ij}}$: the rank of $\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|$ in ascending order, counting outward from $\mathbf{x}_{i}$, that is
$$ s_{_{ij}}=\#\left\{1\le k\le n\,\,:\,\left\|\mathbf{x}_{i}-\mathbf{x}_{k}\right\|\leq\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|\right\}$$

$r_{ij}$: the rank of $\left\|\mathbf{y}_{i}-\mathbf{y}_{j}\right\|$ in ascending order, counting outward from $\mathbf{y}_{i}$, that is
$$ r_{_{ij}}=\#\left\{1\le k\le n\,\,:\,\left\|\mathbf{y}_{i}-\mathbf{y}_{k}\right\|\leq\left\|\mathbf{y}_{i}-\mathbf{y}_{j}\right\|\right\}$$

$\widehat{r}_{_{ij}}$: the rank of $\left\|\widehat{\mathbf{y}}_{i}-\widehat{\mathbf{y}}_{j}\right\|$ in ascending order, counting outward from $\widehat{\mathbf{y}}_{i}$, that is
$$ r_{_{ij}}=\#\left\{1\le k\le n\,\,:\,\left\|\widehat{\mathbf{y}}_{i}-\widehat{\mathbf{y}}_{k}\right\|\leq\left\|\widehat{\mathbf{y}}_{i}-\widehat{\mathbf{y}}_{j}\right\|\right\}$$

\bigskip

\item $N^{^I}_{_J}(i)$: the index set of $J$-nearest neighbors of $\mathbf{x}_{i}$,
that is $N^{^I}_{_J}(i)=\left\{ j\,|\,1 \leq s_{_{ij}}\leq J \right\}$.

$N^{^O}_{_J}(i)$: the index set of $J$-nearest neighbors of $\mathbf{y}_{i}$, that is
$N^{^O}_{_J}(i)=\left\{ j|1 \leq \widehat{r}_{_{ij}}\leq J \right\}$.

Here the superscripts I/O stand for input/output spaces of the algorithms.
\bigskip

\item $\# A$: the cardinality of the set $A$.
\end{itemize}

In the problem setup \eqref{eq1}, where $\varphi$ is assumed to be an isometry, we can suppose that geodesic distances are well approximated by Euclidean distances on the local level.  Therefore, we might ideally assume that for any point $\mathbf{x}_{i}$, there exists a set $N(i)$ of neighboring cases such that the nearest neighbor ranks of the latent low-dimensional data $\mathbf{Y}$ are preserved in $\mathbf{X}$, i.e.
$$r_{_{ij}}=s_{_{ij}}.$$

Therefore, a low-dimensional representation $\widehat{\mathbf{Y}}$ can be said to have rank fidelity if $\widehat{\mathbf{Y}}$ also preserves such ranks, i.e.
$$\widehat{r}_{_{ij}}=s_{_{ij}}.$$
This identity assumes that there is a bijection between the data points in the input and output neighborhoods. In practice when perfect rank fidelity is not achieved then no bijection can be assumed. Therefore two types of errors could occur due to the mapping $\widehat{\psi}$.
\begin{itemize}
\item Output error: The changes of nearest neighbor ranks $\widehat{r}_{_{ij}}$ from the output space to the input space.

\item Input error: The changes of nearest neighbor ranks $s_{_{ij}}$ from the input space to the output space.
\end{itemize}

These two types of errors can be measured by a local rank correlation between the nearest neighbor distances in the input and output spaces. For a given $i$, and all $j$ in $N_{J}^{I}(i) \bigcup N_{J}^{O}(i)$, define the trimmed rank
\begin{align}
\label{ad1}
& S_{_{ij}}=
\begin{cases}
\delta_{_{ij}}, & j \in N_{J}^{I}(i) \bigcap N_{J}^{O}(i)\\
\frac{\zeta+J+1}{2}, & j \notin N_{J}^{I}(i) \bigcap N_{J}^{O}(i)\,,
\end{cases}\\
\label{ad2}
& \widehat{R}_{_{ij}}=
\begin{cases}
\widehat{\delta}_{_{ij}}, & j \in N_{J}^{I}(i) \bigcap N_{J}^{O}(i)\\
\frac{\zeta+J+1}{2}, & j \notin N_{J}^{I}(i) \bigcap N_{J}^{O}(i)\,,
\end{cases}
\end{align}
where
\begin{align*}
& \delta_{_{ij}}=\#\left\{k\in N^{I}_{_J}(i) \bigcap N^{O}_{_J}(i):\,\left\|\mathbf{x}_{i}-\mathbf{x}_{k}\right\|\leq\left\|\mathbf{x}_{i}-\mathbf{x}_{j}\right\|\right\}\\
& \widehat{\delta}_{_{ij}}=\# \left\{k\in N^{I}_{_J}(i) \bigcap N^{O}_{_J}(i):\,\left\|\widehat{\mathbf{y}}_{i}-\widehat{\mathbf{y}}_{k}\right\|\leq\left\|\widehat{\mathbf{y}}_{i}-\widehat{\mathbf{y}}_{j}\right\|\right\}\\
& \zeta=\# \left(N_{J}^{I}(i) \bigcap N_{J}^{O}(i) \right).
\end{align*}
The trimming in \eqref{ad1} and \eqref{ad2} is to make the ranks comparisons local. To measure the output error, we can define the following.

\theoremstyle{definition}
\newtheorem{defi}{Definition}
\begin{defi} {\it Local rank correlation for output error }:

For a given point $\mathbf{x}_i$, define the local output Spearman correlation as
\begin{align}
\label{rhoO}
\rho^{O}_{_J}(i,\mathbf{X},\widehat{\mathbf{Y}})=1-\frac{6\left( \underset{j \in N_{_J}^{O}(i)}{ \sum}\left\{\left(S_{_{ij}}-\widehat{r}_{_{ij}}\right)^{2}\right\}+U\right)}{J(J^{2}-1)}\,,
\end{align}
where $U=\left[(J-\zeta)^{3}-(J-\zeta)\right]/12$ is the adjustment made for the appearance of ties \citep{kendall1}. Define the local output Kendall correlation as
\begin{align}
\label{tauO}
\tau^{O}_{_J}(i,\mathbf{X},\mathbf{Y})=\frac{\underset{j<k \in N_{_J}^{O}(i)}{\sum}2\,\mathrm{sign}\left\{(S_{_{ij}}-S_{_{ik}})\cdot(\widehat{r}_{_{ij}}-\widehat{r}_{_{ik}})\right\}}{J(J-1)}.
\end{align}

For a given input dataset $\mathbf{X}$ and a given dimensionality reduction method $\widehat{\psi}: \mathbf{X} \mapsto \widehat{\psi}(\mathbf{X})$, an overall goodness measure can be defined by averaging the local correlation over all cases in the sample.
\begin{align}
\label{GO}
G^{O}_{_J}(\widehat{\psi},\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}\Gamma^{O}_{_J}\left( i,\mathbf{X},\widehat{\psi}(\mathbf{X}) \right),
\end{align}
where $\Gamma^{O}_{_J}$ can be either $\rho^{O}_{_J}$, or $\tau^{O}_{_J}$.
\end{defi}

The local rank correlations $\rho^{O}_{_J}(i)$ or $\tau^{O}_{_J}(i)$ measure the similarity, in terms of output errors, between the corresponding neighborhoods, $N_{_J}^{I}(i)$ and $N_{_J}^{O}(i)$. Similarly, we can define local rank correlations to measure the input error.

\begin{defi} {\it Local rank correlation for input error} :

Given an input dataset $\mathbf{X}$ and a low-dimensional representation $\widehat{\mathbf{Y}}$, the local Spearman correlation and local Kendall correaltion for the input error at the $i$-th case are defined as
\begin{align}
& \rho^{I}_{_J}(i,\mathbf{X},\widehat{\mathbf{Y}})=1-\frac{\underset{j \in N_{_J}^{I}(i)}{\sum}\left\{\left(s_{_{ij}}-\widehat{R}_{_{ij}}\right)^{2}\right\}+U}{\frac{1}{6}(J^{3}-J)}\label{rhoI}\,,\\
& \tau^{I}_{_J}(i,\mathbf{X},\widehat{\mathbf{Y}})=\frac{\underset{j<k,\, j,\,k \in N_{_J}^{I}(i)}{\sum}\mathrm{sign}\left\{(\widehat{R}_{_{ij}}-\widehat{R}_{_{ik}})\cdot(s_{_{ij}}-s_{_{ik}})\right\}}{\frac{1}{2}J(J-1)}\label{tauI}\,.
\end{align}

The overall goodness measure of a given method $\widehat{\psi}$ and input data $\mathbf{X}$ is defined as
\begin{align}
\label{GI}
G^{I}_{_J}(\widehat{\psi},\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}\Gamma^{I}_{_J}\left( i,\mathbf{X},\widehat{\psi}(\mathbf{X}) \right)\,,
\end{align}
where $\Gamma^{I}_{_J}$ can be either $\rho^{I}_{_J}$, or $\tau^{I}_{_J}$.
\end{defi}
\subsection{Remark}

The proposed local rank correlations have some nice properties. The higher values of local measures $\Gamma^{I}_{_J}(i)$ and $\Gamma^{O}_{_J}(i)$ indicate a higher degree of similarity between the original data and the low-dimensional configuration in the neighborhood of case $i$, while values close to 0, or negative values indicate that low-dimensional configuration fails to preserve the local structure of the input data in certain neighborhoods. Two special situations are:
\begin{itemize}
\item $\Gamma^{I}_{_J}(i)=\Gamma^{O}_{_J}(i)=1$ if all the ranking relationships of the observed data $\mathbf{X}$ in the neighborhood of case $i$ are preserved exactly in the corresponding neighborhood in the output data $\widehat{\mathbf{Y}}$.

\item The expected values $\mathrm{E}\left[\,\Gamma^{I}_{_J}(i)\,\right]$ and $\mathrm{E}\left[\,\Gamma^{O}_{_J}(i)\,\right]$ are both zero, for any case $i$, where the output $\widehat{\mathbf{Y}}$ is generated by an algorithm which is stochastically independent of the input data $\mathbf{X}$.
\end{itemize}
These two facts hold for both local Spearman and Kendall correlations.
Notice that the second situation is worse than we can have in practice. Moreover,  the local measures $\Gamma^{I}_{_J}(i)$ and $\Gamma^{O}_{_J}(i)$, can achieve negative values for some $i$. Nevertheless, the overall goodness measures $G_{_J}^{O}$ and $G_{_J}^{I}$, for sensible algorithms, will take values between 0 and 1. We remind the reader that the use of ranks is to protect against non-normality and extreme skewness of distance distributions in high dimensions.

The computational complexity is also of interest. To calculate the goodness measure, we first construct the $J$-nearest neighbor graph for both $\mathbf{X}$ and $\widehat{\mathbf{Y}}$. This step scales as $O(n^{2}p)$. In the next step, we calculate the local rank correlation in each neighborhood. This scales (in each neighborhood) as $O(J)$ for Spearman $\rho_{_J}$ and $O(J\log{J})$ for Kendall $\tau_{_J}$. Therefore, since $J\leq n$, the total complexity of calculating $G^{I}_{_J}$ (or $G^{O}_{_J}$) for $\rho_{_J}$ scales as $O(n^{2}p)$. The total complexity of calculating $G^{I}_{_J}$ (or $G^{O}_{_J}$) for $\tau_{_J}$ scales as  $O(n^{2}p+nJ\log J)$.

To use the proposed goodness measure $G_{_J}$ for assessing the performance of a dimension reduction method, four local measures can be chosen. We may choose either $\Gamma^{I}_{_J}$ or $\Gamma^{O}_{_J}$, and we may also  choose to use either Spearman $\rho_{_J}$ or Kendall $\tau_{_J}$. The measures $\Gamma^{I}_{_J}$ and $\Gamma^{O}_{_J}$ quantify different types of errors in dimension reduction. Although these two types of errors usually occur together, having both $G^{I}_{_J}$ and $G^{O}_{_J}$ provides more complete information about the performance of a given method.

\subsection{Choice of \texorpdfstring{$J$}{J}}

In the proposed measures, $J$ is a user-specified tuning parameter, which specifies the neighborhood size for local rank comparisons. Notice that some nonlinear dimensionality reduction methods start with a $K$-nearest neighbors graph, and $K$ is also a user-specified parameter. The choice of $J$ in the local rank correlation does not have to depend on the value of $K$. Ideally, $J$ is chosen sufficiently small that there exists a $q$-flat $F^q$ such that the data points in $N_{_J}^I(i)$ lie approximately on $F^q$. In addition, $J$ is ideally chosen sufficiently large that each $N_{_J}^I(i)$ is informative about the local geometrical characteristics of the data. In practice, we expect $G_{_J}(\widehat{\psi},\mathbf{X})$ to be large and stable for values of $J$ sufficiently small to satisfy the former condition. So we can choose $J$ by plotting $G_{_J}(\widehat{\psi},\mathbf{X})$ against $J$ and choosing the largest value of $J$ for which the stability is observed.

\subsection{Adjustments for output-normalized methods}

Normalization of the output, as used in algorithms such as Local Linear Embedding \citep{lle}, Laplacian Eigenmap \citep{LEM}, and Local Tangent Space Alignment \citep{ltsa}, affinely transforms the geometrical structure of neighborhoods, so that the distance relationships will not be preserved between the input and the output configurations. It is not adequate to check the correspondence of nearest neighbor ranks between $\mathbf{X}$ and $\widehat{\psi}(\mathbf{X})$ for those output-normalized methods. Instead, we will look for a transformation matrix $\widehat{\mathbf{A}}_{q\times q}$, and assess the performances of output-normalized methods by an adjusted measure
\begin{align}
\label{adj}
G^{\mathbf{A}}_{_J}(\widehat{\psi},\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}\Gamma_{_J}(i,\mathbf{X},\widehat{\psi}_{_\mathbf{A}}(\mathbf{X}))\,,
\end{align}
where $\widehat{\psi}_{_\mathbf{A}}(\mathbf{X})=\widehat{\psi}(\mathbf{X})\cdot \widehat{\mathbf{A}}$, and $\Gamma_{_J}$ can be $\rho^{I}_{_J}$, $\rho^{O}_{_J}$ or $\tau^{I}_{_J}$, $\tau^{O}_{_J}$.

It is hoped that after the affine transformation $\widehat{\mathbf{A}}$, $\widehat{\psi}_{_\mathbf{A}}(\mathbf{X})$ can preserve the proximities between neighboring points as much as possible, i.e. $\widehat{\mathbf{A}}$ will minimize the least squared error
$$\sum_{i=1}^{n}\sum_{j \in N^{^I}_{_J}(i)}\left[ (\mathbf{x}_{i}-\mathbf{x}_{j})^{\prime}(\mathbf{x}_{i}-\mathbf{x}_{j})-(\mathbf{y}_{i}-\mathbf{y}_{j})^{\prime}\mathbf{A}^{\prime}\mathbf{A}(\mathbf{y}_{i}-\mathbf{y}_{j}) \right]^{2}\,,$$
where $\mathbf{x}_{i}$ and $\mathbf{y}_{i}$ are the corresponding points in the original data and in the output of the algorithm $\widehat{\psi}$, respectively.

\section{Experiments}
\label{sec4}

In this section, we conduct numerical examples on some benchmark datasets to illustrate the effectiveness of the local rank correlation.

\theoremstyle{definition}
\newtheorem{example}{Example}
\begin{example} {\it The Swiss roll and the S-curve:}
In this experiment, $n=1000$ data points are generated randomly from two manifolds, the Swiss roll and the S-curve. They are both 2-dimensional manifolds embedded into $\mathbb{R}^{3}$ (Figure \ref{gen}). The data
points are colored to help readers recognize the structure of the manifolds. Among many dimension reduction methods, we choose four diverse methods, namely   ISOMAP, LTSA, MVU, and PCA. Figure \ref{swiss} and Figure \ref{Scurve} show four output configurations in $\mathbb{R}^2$ from these methods for the Swiss roll and the S-curve, respectively. We evaluate the performance of four methods by the local rank correlations $G^{I}_{_J}(\widehat{\psi},\mathbf{X})$ and $G^{O}_{_J}(\widehat{\psi},\mathbf{X})$ with both Spearman $\rho_{_J}$ and Kendall $\tau_{_J}$. Notice that LTSA is an output-normalized method, and therefore, its performance is assessed by the adjusted measures \eqref{adj}.

The goodness measures are calculated under different values of $J$. Figure \ref{swiss-lcorr} and Figure \ref{scurve-lcorr} show the values of $G_{_J}(\widehat{\psi},\mathbf{X})$ for each method as functions of $J$.

\begin{figure}[!ht]
	\centering
		\includegraphics[width=0.8\textwidth]{generate2.pdf}
	\caption{Data on the Swiss roll and the S-curve.}
	\label{gen}
\end{figure}

\begin{figure}[!ht]
	\centering
		\includegraphics[width=0.8\textwidth]{generate-swiss.pdf}
	\caption{Two-dimensional output configurations of different methods for the Swiss roll.}
	\label{swiss}
\end{figure}

\begin{figure}[!ht]
	\centering
		\includegraphics[width=0.8\textwidth]{generate-Scurve.pdf}
	\caption{Low-dimensional configurations of different methods for S-curve.}
	\label{Scurve}
\end{figure}

\begin{figure}[!ht]
	\centering
		\includegraphics[width=1\textwidth]{swiss-J.pdf}
	\caption{Local Spearman correlation as functions of $J$.}
	\label{swiss-lcorr}
\end{figure}

\begin{figure}[!ht]
	\centering
		\includegraphics[width=1\textwidth]{scurve-J.pdf}
	\caption{Local Spearman correlation as functions of $J$.}
	\label{scurve-lcorr}
\end{figure}

As can be seen from Figure \ref{swiss} and Figure \ref{Scurve}, in the 2-dimensional configurations from PCA in both the Swiss roll and the S-curve, the points with different colors are mixed together, because PCA fails to recover the nonlinear structure of the embedded data. Among three nonlinear methods, the configurations from MVU preserve the structure to some extent. Points with different colors are reasonably separated in the middle, but they mix a little at boundaries. Both LTSA and ISOMAP preserve the color level well, indicating a better embedding than MVU and PCA. These facts are all correctly reflected by the four goodness measures $G_{_J}(\widehat{\psi},\mathbf{X})$ in Figure \ref{swiss-lcorr} and Figure \ref{scurve-lcorr}. Also, both figures show that all the four measures are stable within a reasonable range of $J$.

We also compare the local rank correlation (LRC) with the goodness measures, MRREs, T\&C, and LCMC (all with $J=6$). The results are reported in Table \ref{tab-swiss} and Table \ref{tab-Scurve}.
\begin{table}[!ht]
\large
    \begin{center}
        \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
        \hline
\multirow{2}{*}{Methods} & \multicolumn{4}{|c|}{LRC} & \multicolumn{2}{|c|}{MRREs} & \multicolumn{2}{|c|}{T \& C} & \multirow{2}{*}{LCMC}\\ \cline{2-9}
& $\rho^{I}_{_J}$ & $\rho^{O}_{_J}$ & $\tau^{I}_{_J}$ & $\tau^{O}_{_J}$ & $M_{_O}$ & $M_{_I}$ & $T$ & $C$ & \\
            \hline \hline
ISOMAP & 0.787 & 0.782 & 0.701 & 0.698 & 0.999  &  0.999  &  0.999  &  0.999  & 0.894\\
LTSA & 0.988 & 0.978 & 0.981 & 0.975 & 0.999  &  0.998 &   0.993  &  0.998  &  0.609\\
MVU &  0.703 & 0.623 & 0.653 & 0.578 & 0.999  &  0.999  &  0.996  &  0.999  &  0.828\\
PCA &  0.594 & 0.198 & 0.483 & 0.171 & 0.998  &  0.997  &  0.883  &  0.995  &  0.415\\              \hline \hline
        \end{tabular}
    \end{center}
	\caption{Assessing ISOMAP, LTSA, MVU, PCA in Swiss Roll data ($J=6$).}
	\label{tab-swiss}
\end{table}
\begin{table}[!ht]
\large
    \begin{center}
        \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|}
        \hline
\multirow{2}{*}{Methods} & \multicolumn{4}{|c|}{LRC} & \multicolumn{2}{|c|}{MRREs} & \multicolumn{2}{|c|}{T \& C} & \multirow{2}{*}{LCMC}\\ \cline{2-9}
& $\rho^{I}_{_J}$ & $\rho^{O}_{_J}$ & $\tau^{I}_{_J}$ & $\tau^{O}_{_J}$ & $M_{_O}$ & $M_{_I}$ & $T$ & $C$ & \\
            \hline \hline
ISOMAP & 0.816 & 0.804 & 0.763 & 0.803 & 0.999 & 0.999  &  1.000  &  1.000  &  0.891\\
LTSA & 0.994 & 0.993 & 0.983 & 0.979 & 0.999  & 0.999 &  0.999  &  0.999  &  0.867\\
MVU & 0.721 & 0.646 & 0.695 & 0.617 & 0.999  & 0.998 &  0.993  &  0.997  &  0.754\\
PCA & 0.673 & 0.375 & 0.388 & 0.369 & 0.998  &  0.998  &  0.963   & 0.998  &  0.584\\              \hline \hline
        \end{tabular}
    \end{center}
	\caption{Assessing ISOMAP, LTSA, MVU, PCA in S-curve data ($J=6$).}
	\label{tab-Scurve}
\end{table}

We now turn to the analysis of the Swiss roll data.  We examine the output visually first.  As expected, PCA works poorly. This is a consequence of the many-to-one nature of a linear projection of the Swiss roll.  All three of the other algorithms separate colors well.  However, an ideal output should be perfectly rectangular or square.  Visually, we prefer LTSA to ISOMAP and prefer ISOMAP to MVU.

We next turn to the performance measures for the Swiss roll output. The LCMC performance measure has identified the failure of PCA to account for the nonlinearity of the Swiss roll.  However, it also ranks the outcome of MVU as superior to LTSA, which is not visually supported in Figure \ref{swiss}.  The ISOMAP algorithm is ranked highest by LCMC, but only slightly higher than MVU.  Despite the fact that the output from LTSA and ISOMAP are visually close, the LCMC measures are quite distinct.  It would seem reasonable to assume that if LCMC is picking up problems with these algorithms, the problems are not visually obvious.

The trustworthiness and continuity measures, present a more complex picture.  As expected, all algorithms perform reasonably well on the continuity criterion.  For example, the linear projection defining the outcome of PCA is a continuous mapping of the Swiss roll, and therefore satisfies the continuity criterion well with a high value of C.  That is also the case for the other algorithms.  The trustworthiness measure clearly separates out PCA as the least desirable algorithm, as expected.  Other algorithms perform very well and similarly to each other.  Once again, MVU is ranked higher than LTSA, in contrast to our visual interpretation.

Turning to MRREs, we see that very little separation can be seen among the algorithms.  Since these measures are not standardized, we must be wary of drawing too many conclusions from the proximity of these values to one.  Nevertheless, MVU is again ranked higher than LTSA and is not separated in performance from ISOMAP.

We are happy with the correspondence between LRC and our conclusions based upon visual inspection.  The LTSA algorithm is ranked highest by all four measures of LRC, followed by ISOMAP, MVU and PCA.  We find it especially helpful that there is no inconsistency in ranking between the four measures of LRC.

Our conclusions, both visual and quantitative, for the S-curve trend in a similar direction to the Swiss roll.  Visually, we rank ISOMAP and LTSA best with little to choose between them.  The MVU algorithm works very well, but performs poorly close to the boundary, where a mixing of colors can be seen.  The PCA algorithm is, once again, the worst of the four, but within an acceptable standard as the output is not from a many-to-one projection.

Next, we examine the performance measures for the S-curve output.  Note that LCMC performs in close correspondence with our visual analysis.  Similar remarks hold for T \& C.  Although continuity ranks PCA better than MVU, we should note that PCA is a continuous algorithm and that this is reasonable.  Trustworthiness rankings correspond to the visual rankings.  The discrimination provided by MRREs is, again, unclear.  Nevertheless it is consistent with visual conclusions.  It is reasonable to conclude from the Swiss roll data and S-curve data that MRREs are less useful for these two benchmark examples.  Under the LRC measures, LTSA is consistently the best, followed by ISOMAP, MVU and PCA in that order.  This agrees with the visual analysis.
\end{example}

\section{Choosing tuning parameters for algorithms}

Many nonlinear dimension reduction methods start with constructing the $K$-nearest neighbor graph. The success of graph-based nonlinear dimensionality reduction methods depends heavily on the selection of $K$.
If $K$ is chosen to be too small, the local geometric structure cannot be accurately represented in the neighborhood graph. On the other hand, if $K$ is chosen to be too large, the $K$-nearest neighbor graph will
contain \textit{shortcuts}, i.e. two points will be mistakenly considered as neighbors when they are in fact far away on the manifold. In practice, $K$ is usually chosen by experience or trial and error.

The local rank correlation can help in choosing $K$, because $G_{_J}(\widehat{\psi},\mathbf{X})$ measures the performance of $\widehat{\psi}$. We may apply the algorithm over a range of values of $K$, and calculate $G_{_J}(\widehat{\psi},\mathbf{X})$ as a function of $K$ (as shown in Figure \ref{fig51}). We then pick the $K$ that corresponds to the largest $G_{_J}(\widehat{\psi},\mathbf{X})$.

\begin{example}{\it Selecting neighborhood size $K$ in ISOMAP:}
Here we consider the performance of ISOMAP on the Swiss roll manifold. We demonstrate that it is risky to make a desultory choice of $K$, and how local rank correlation can solve this problem. The data are generated randomly on the Swiss roll manifold with sample size $n=1500$. The ISOMAP algorithm is applied on the data with different values of $K$, and Figure \ref{fig52} shows the respective low-dimensional configurations. In each case, the performance is evaluated by the local rank correlation and displayed as a function of $K$ in Figure \ref{fig51}.

\begin{figure}[!ht]
	\centering
		\includegraphics[width=0.8\textwidth]{2.pdf}
	\caption{Low-dimensional configurations with different values of $K$.}
	\label{fig52}
\end{figure}

\begin{figure}[!ht]
	\centering
		\includegraphics[width=1\textwidth]{swiss-k.pdf}
	\caption{Local rank correlation as a function of $K$ ($J=6$).}
	\label{fig51}
\end{figure}

In Figure \ref{fig51}, the left panel shows $G^{I}_{_J}$ and $G^{O}_{_J}$ with Spearman $\rho_{_J}$, and the right panel shows $G^{I}_{_J}$ and $G^{O}_{_J}$ with Kendall $\tau_{_J}$.
As can be easily noticed in Figure \ref{fig52}, the performance of ISOMAP gets better as $K$ increases from $K=7$ to $K=13$. A crucial change has happened at points $K=13$ and $K=14$. In these two situations, the neighborhood sizes only differ by 1 but the corresponding configurations suddenly become unsatisfactory (at $K=14$).
The fact is correctly captured by the local rank correlation and reflected in Figure \ref{fig51}. In all four measures, we observe a peak at $K=13$, and a steep drop at $K=14$. For nonlinear methods which contain the neighborhood size $K$ as a tuning parameter, it is often desirable to choose a relatively large value of $K$ to get a better embedding. On the other hand, a value of $K$ which is too large will invalidate the procedure.

\end{example}

\section{Estimating the intrinsic dimensionality of a manifold}

Another key parameter in dimension reduction algorithms is the intrinsic dimensionality $q$. The local rank correlation can be applied to help in estimating the intrinsic dimensionality. The idea is that if the dimensionality of the low-dimensional representation is chosen to be too small, important features of the original data might be ``collapsed'' onto the same dimensions so that the topological structure cannot be preserved very well. As the dimensionality of the representation increases, while remaining below the correct dimension $q$ of the manifold, the local rank correlation should increase. On the other hand, when the dimensionality of the representation is greater than $q$, the only additional information in the data provided by the additional dimensions would be noise. Therefore, provided that the noise is small, the local rank correlation would become stable at values larger than $q$. This is similar to the scree plot (\cite{ScreePlot}), but in reverse, used for choosing dimensionality of linear manifolds in PCA.
In practice, for a given dataset $\mathbf{X}$ and a chosen method $\widehat{\psi}$, we may apply the method with different dimensions, and evaluate the performance of $\widehat{\psi}$ by $G_{_J}(\widehat{\psi},\mathbf{X})$. We estimate the intrinsic dimensionality by plotting $G_{_J}(\widehat{\psi},\mathbf{X})$ as a function of dimension, and choosing the value $q$, beyond which $G_{_J}(\widehat{\psi},\mathbf{X})$ becomes stable.

\begin{example}{\it Estimating the intrinsic dimensionality of the sculpture face data}:
The sculpture face dataset (\cite{isomap}) includes 698 images, each image having $64\times 64$ pixels of a sculpture face while varying three free parameters: left-right pose, up-down pose, and lighting direction. So the data are originally in $\mathbb{R}^{64\times 64}$.  The sculpture face images are recorded as $64 \times 64$ vectors. Since images are taken from the same sculpture face by varying three parameters, the intrinsic dimensionality of the manifold on which these data vectors lie is three. We apply the ISOMAP algorithm with different values of $q$ having chosen the neighborhood size $K=8$. The local rank correlations are calculated as functions of $q$. In Figure \ref{fig54}, the left panel shows $G^{I}_{_J}$ and $G^{O}_{_J}$ with Spearman $\rho_{_J}$, and the right panel shows $G^{I}_{_J}$ and $G^{O}_{_J}$ with Kendall $\tau_{_J}$. As can be seen, all four curves become stable beyond $q=3$, based on which we estimate the intrinsic dimensionality to be $\widehat{q}=3$.

\begin{figure}[!ht]
	\centering
		\includegraphics[width=1\textwidth]{intrinsic.pdf}
	\caption{Local rank correlation as a function of dimensionality $q$ ($J=6$).}
	\label{fig54}
\end{figure}

\end{example}

\section{Conclusions}
To quantify the performance of a dimension reduction method, we introduced a family of local rank correlation measures, which are
easily interpretable and motivated by their robustness properties for nearest neighbor distributions in high dimensions. We found that the local rank correlation closely corresponds to our visual interpretation of the quality of the output in benchmark examples.
We also demonstrated that the local rank correlation can be applied to the problems of estimating the intrinsic dimensionality of the original data, and selecting appropriate values for  tuning parameters.

\bibliographystyle{apa}
\cleardoublepage

\phantomsection

\addcontentsline{toc}{chapter}{\textbf{References}}