
\chapter{High Dimensional Asymptotics\label{chapM:HDasy}}

A key insight of classical mathematical statistics is that asymptotic
analysis\index{asymptotic analysis} provides a useful tool for understanding
the fundamentals that underlie complicated statistical contexts. Limiting
operations as the sample size $n\rightarrow\infty$ have been a workhorse
approach that has revealed very many foundational statistical discoveries,
e.g. those related to laws of large numbers and the central limit
theorem. These have been so useful that some have developed the mistaken
idea that asymptotics should only be about growing sample sizes and
their impact on statistical analyses. A central theme of this chapter
is that in the modern era of very large and complex data sets, many
less familiar but very useful insights are available from the use
of more general limiting operations involving \emph{both} the sample
size $n$ and the dimension $d$.

As noted in \citet{Shen2016HDLSS,Shen2016GenFrame}, there are several
\emph{domains\index{asymptotic domains}} of quite different and important
asymptotic insights available, that are usefully indexed by the relationship
between the sample size $n$ and the dimension $d$. Important early
work on the asymptotics of high dimensional statistical methods was
done by \citet{Portnoy1984HDMest1,Portnoy1985HDMest2,Portnoy1988HDexpFam},
who studied various robust\index{robust methods} methodologies using
the limit as $n\rightarrow\infty$ together with the dimension $d$
growing at the rate $d\sim n^{1/2}$. A different domain is the large
literature on random matrix theory\index{random matrix theory|(},
where $n$ and $d$ grow at the same rate, some basics of which are
briefly reviewed in Section \ref{secM-A:RandomMatrix}. Perhaps surprising
is the idea that anything useful at all can be shown in yet another
domain where $d\rightarrow\infty$ for fixed $n$ (i.e. growing and
hence increasing noise, but a fixed number of data objects), as discussed
in Section \ref{secM-B:HDLSS}. A few ideas about the case where both
grow, but $d$ grows more quickly than $n$ are studied in Section
\ref{secM-C:HDMSS}.

It is worth noting that at the time of this writing, a number of the
ideas and views expressed here are directly contrary to the currently
prevailing opinion of the theoretical statistics community. The latter
is perhaps summarized in Section 8.3 of \citet{Wainwright2019HiDimStat}:
`does PCA still perform well in the high-dimensional regime $n<d$?
The answer to this question turns out to be a dramatic \textquotedblleft no\textquotedblright .'
That view notably clashes with the experience of practioners who actually
analyze data: PCA frequently gives useful visual insights for very
high dimensional data. Clear examples in this book include the Lung
Cancer data in Section \ref{subsecB.2.A:LungCancerData}, the Pan
Cancer data in Section \ref{subsecB.2.B:PanCancerExample}, and the
GWAS data in Section \ref{secO.B.B:GWAS}. The resolution of this
apparent paradox, as well as some others also based on such misinterpretations
of the asymptotic theory, can be found in Section \ref{secM-B:HDLSS}.

\section{Random Matrix Theory\label{secM-A:RandomMatrix}}

Random matrix theory is centered on the eigenvalues of the sample
covariance of a $d\times n$ matrix of independent standard Gaussian
random variables. Motivated by a model for the nuclei of heavy atoms
in nuclear physics, \citet{Wigner1947} first studied square $d\times d$
matrices of more general independent random variables (with common
variance) as $d\rightarrow\infty$, including providing a semi-circle
law for the distribution of the eigenvalues of the random matrix.
\citet{MarcenkoPastur1967} explored the asymptotic distribution of
the eigenvalues (studied around (\ref{eqP-O:EigenAnalofCov})) of
sample covariance matrices (defined in (\ref{eqB.F:SampleCovariance})).
Those matrices are based on i.i.d. $N_{d}\left(\boldsymbol{0}_{d,1},\sigma^{2}\boldsymbol{I}_{d}\right)$
(using the notation of (\ref{eqB.G:0sMatrix}), $\boldsymbol{0}_{d,1}$
is the $d\times1$ vector of zeros) random vectors $\widetilde{\boldsymbol{z}}_{1},\cdots,\widetilde{\boldsymbol{z}}_{n}\in\mathbb{R}^{d}$.
The limiting operation lets both $n\to\infty$ and $d\rightarrow\infty$
with $\frac{d}{n}\rightarrow c$. The density of the asymptotic eigenvalue
distribution is given in (\ref{eqM.I:MarcPastLSD}). That distribution
tends to \index{Wigner's semi-circle law}Wigner's semi-circle law
as $c\rightarrow0$ and also (ignoring zero eigenvalues) as $c\rightarrow\infty$.

Figure \ref{figM.D:Scree2MP} provides an illustration of the basic
ideas. For simplicity, all panels are based on a $d\times n$ matrix
$\widetilde{\boldsymbol{X}}$ of independent standard Gaussian random
variables where the common variance $\sigma^{2}$ is taken to be $1$.
An approximation of the sample covariance matrix, which is particularly
simple to work with, is the outer product $\widetilde{\boldsymbol{\varSigma}}=\frac{1}{n}\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{X}}^{t}$.
This is not quite the $d\times d$ sample covariance matrix (when
thinking of the columns of $\widetilde{\boldsymbol{X}}$ as data objects),
because the sample mean has not been subtracted. However, because
the expected value of each entry is $0$, this is still a reasonable
estimate of the underlying theoretical covariance matrix, $\boldsymbol{\varSigma}=\boldsymbol{I}_{d}$.
In fact it is the maximum likelihood estimate of $\boldsymbol{\varSigma}$
under the assumption that the sample mean is $\boldsymbol{0}_{d,1}$.
The left panel of Figure \ref{figM.D:Scree2MP} shows scree plot\index{scree plot|(}s
(recall from Figure \ref{figB-Y:Screeplots} this is a plot of the
sorted eigenvalues of $\widetilde{\boldsymbol{\varSigma}}$ as a function
of the PC index $j=1,\cdots,d$) in the case $d=100$, for several
values of $n$. Since $\widetilde{\boldsymbol{\varSigma}}$ estimates
$\boldsymbol{I}_{d}$, each eigenvalue can be viewed as an estimate
of $1$, shown as the horizontal dashed line. For the very large $n=10,000$
(shown in green) the eigenvalues are rather close to $1$, although
there is some natural variation as seen. That variation grows as $n$
is reduced to $1000$ (shown in red) and to $300$ (blue). There is
a growing asymmetry which becomes especially strong at $n=100$ (black).
The asymmetry can be thought of as being caused by the constraint
that eigenvalues must be non-negative (actually positive in this case).

While the left panel of Figure \ref{figM.D:Scree2MP} explores differing
matrix shapes, the case of constant matrix shape is considered in
the center panel. The shape is constant in the sense that $\frac{d}{n}=c$,
with $c=\frac{1}{5}$ there. While sampling variation is visible for
the magenta eigenvalues for $n=500$ and $d=100$, it diminishes for
the cyan case where $n=2000$ and $d=400$, and looks essentially
like a smooth curve for the yellow $n=8000$ and $d=1600$. To make
this apparent convergence occur, the horizontal axis has been rescaled
to $\frac{j}{d}$ for each case. The right panel shows a useful conceptual
viewpoint of the limit of this process, by studying the distribution
of the cyan ($n=2000$ and $d=400$) eigenvalues. The vertical axis
shows the binning of these into a 50 bin histogram. That histogram
approximates the $c=\frac{1}{5}$ version of the \index{Marchenko-Pastur distribution|(}Mar\v{c}enko-Pastur
distribution, which essentially reflects the limit of the apparent
curves in the center panel as $n,d\rightarrow\infty$ (with $\frac{d}{n}=\frac{1}{5}$).
Good overview of many aspects of Mar\v{c}enko-Pastur theory can be
found in the monographs \citet{BaiSilverstein2010RandMat} and \citet{Yao2015RandMat}.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterM/OODAbookChpMFigDScree2MP}

\caption{Scree plots motivating Mar\v{c}enko-Pastur distributions, using i.i.d.
standard Gaussian samples. Left panel illustrates differing distributional
shape over $n=100$ (black), $n=300$ (red), $n=1000$ (blue) and
$n=10,000$ (green), all for $d=100$. Center panel shows stability
of distribution for common $c=\frac{d}{n}=\frac{1}{5}$, in the cases
$n=500$ (magenta), $n=2000$ (cyan) and $n=8000$ (yellow), using
decreasing symbol size for clear visibility. Right panel shows how
binning of eigenvalues generates a histogram estimate of the Mar\v{c}enko-Pastur
distribution (shown as the thin curve) for the cyan $n=2000$, $d=400$
case. \label{figM.D:Scree2MP}}
\end{figure}
Generally in that limit as $n,d\rightarrow\infty$ with $\frac{d}{n}=c$
the vertical histograms, as illustrated in the right panel of Figure
\ref{figM.D:Scree2MP} converge to the \emph{Limiting Spectral Density}
(LSD) shown as the thin curve, which is the density of the Mar\v{c}enko-Pastur
distribution indexed by $c$. First consider the case $c<1$ (i.e.
$d<n$), where the sample covariance matrix $\widehat{\boldsymbol{\varSigma}}$
has no zero eigenvalues (with probability one). Figure \ref{figM.G:MarcPast}
shows the LSD for several values of $c$. In the more general case
where the Gaussian variance is $\sigma^{2}$, the LSD has the form
\begin{equation}
LSD_{c,\sigma^{2}}\left(x\right)=\frac{1}{2\pi c\sigma^{2}x}\sqrt{\left(b-x\right)\left(x-a\right)}1_{\left[a,b\right]}\left(x\right),\label{eqM.I:MarcPastLSD}
\end{equation}
where $1_{\left[a,b\right]}$ denotes the indicator function\index{indicator function}
(defined at (\ref{eqJ.A:LinearClasser})) on the interval $\left[a,b\right]$,
$a=\sigma^{2}\left(1-\sqrt{c}\right)^{2}$ and $b=\sigma^{2}\left(1+\sqrt{c}\right)^{2}$.
This probability distribution has the perhaps surprisingly simple
expected value $\sigma^{2}$ and variance $c\sigma^{4}$. The $LSD_{c,1}\left(x\right)$
curves (note the case $\sigma^{2}=1$) for several values of $c$
are plotted in Figure \ref{figM.G:MarcPast}. 

The red curve shows the case $c=\frac{1}{3}=\frac{100}{300}$, which
gives the LSD corresponding to the red scree plot in the left panel
of Figure \ref{figM.D:Scree2MP}. The lower density height on the
right corresponds to the increased spacing between red circles on
the left side of the scree plot. The upper end of the support of the
LSD, $b=\left(1+\sqrt{\frac{1}{3}}\right)^{2}\approx2.49$, also fits
well with the top red circle. The higher density to the left of $1.0$
is essentially caused by the eigenvalues being constrained to be positive.
Note that there is a distinct lower support point $a=\left(1-\sqrt{\frac{1}{3}}\right)^{2}\approx0.18$
to the LSD as well, i.e. even in the limit as $n,d\rightarrow\infty$
(with $c=\frac{d}{n}=\frac{1}{3}$) these eigenvalues are bounded
above $0$. Another aspect of the LSD that is visually clear is that
it is vertical at the support points, in the sense that its slope
is infinite there (recall the unit semi-circle has the same property
where it intersects that horizontal axis). This can be mathematically
checked by calculating $\lim_{x\downarrow a}$ and $\lim_{x\uparrow b}$
of the derivative of the LSD. 

\begin{figure}[H]
\includegraphics{ChapterM/OODAbookChpMFigGMarcPast}

\caption{Mar\v{c}enko-Pastur distribution LSDs for examples in the left panel
of Figure \ref{figM.D:Scree2MP} using same colors; black for $c=1$,
red for $c=\frac{1}{3}$, blue for $c=0.1$ and green for $c=0.01$.
Shows wide array of shapes of this distribution. \label{figM.G:MarcPast}}
\end{figure}
The black curve in Figure \ref{figM.G:MarcPast} is the LSD for $c=1$,
which explains the structure of the black scree plot in the left panel
of Figure \ref{figM.D:Scree2MP}. Note that the highest black circle
is around the upper support point $b=\left(1+\sqrt{1}\right)^{2}=4$.
This LSD is very strongly skewed and in fact has a pole at its lower
support point $a=\left(1-\sqrt{1}\right)^{2}=0$. That again can be
checked by calculating the $\lim_{x\downarrow0}LSD_{1,1}(x)$.

The case of $c=0.1=\frac{100}{1000}$ is the blue curve, which models
the blue scree plot in the left of Figure \ref{figM.D:Scree2MP}.
As the larger sample size $n=1000$ results in the blue circles having
much less variation, the support points have moved towards $1$ (indicated
by the vertical dashed line in Figure \ref{figM.G:MarcPast}). In
particular, they are now $a=\left(1-\sqrt{0.1}\right)^{2}\approx0.47$
and $b=\left(1+\sqrt{0.1}\right)^{2}\approx1.73$, which fits with
the top and bottom blue circles in the left panel of Figure \ref{figM.D:Scree2MP}.
Also note that the shape of the LSD has now changed substantially,
in particular it is concave everywhere, versus having the convex region
that was present for the smaller values of $c$.

The support points come even closer to $1$ in the case $c=0.01$
shown as the green curve, again predicting the structure of the green
scree plot in the left of Figure \ref{figM.D:Scree2MP}. Note that
now the shape is becoming rather symmetric, as it converges to a rescaled
version of Wigner's semi-circle\index{Wigner's semi-circle law} distribution
as discussed above. Also the support points $a=\left(1-\sqrt{0.01}\right)^{2}\approx0.81$
and $b=\left(1+\sqrt{0.01}\right)^{2}\approx1.21$ are now becoming
closer to equidistant to 1. 

A parallel Mar\v{c}enko-Pastur theory for the alternate matrix shape
case of $c>1$, i.e. $d>n$, can also be derived. In this case the
sample covariance matrix $\widehat{\boldsymbol{\varSigma}}$ (defined
at (\ref{eqB.F:SampleCovariance})) has (with probability one) $d-n+1$
zero eigenvalues (the $+1$ is due to object mean centering\index{centering}).
Hence the limiting spectral distribution of the eigenvalues is a mixture
of a discrete point mass at $0$ (with weight $1-\frac{1}{c}$) together
with a continuous density which is a rescaling of (\ref{eqM.I:MarcPastLSD})
(with weight $\frac{1}{c}$). Using Singular Value Decomposition\index{Singular Value Decomposition}
ideas as in Section \ref{subsecP-A-B:SVD}, the nonzero eigenvalues
of $\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{X}}^{t}$ are
the same as the eigenvalues of the corresponding \emph{Gram Matrix}
$\widetilde{\boldsymbol{X}}^{t}\widetilde{\boldsymbol{X}}$, i.e,
the matrix of inner products of the columns of $\widetilde{\boldsymbol{X}}$.
The normalization\index{normalization} factor of $\frac{1}{n}$ in
the sample covariance matrix results in diagonal entries which are
sample variance estimates (of the underlying theoretical variances
which are all 1) computed as sums of squares of $n$ Gaussians. In
the Gram Matrix case, the $n$ diagonal entries are now sums of squares
of $d$ Gaussians, so a normalization of $\frac{1}{d}$ is needed
for these $n$ diagonal entries to estimate 1, i.e. the sample variances
essentially estimate $c=\frac{d}{n}$. This calculation reflects the
fact that for $d\gg n$ there is far more variation being summarized
by each of the $n$ nonzero eigenvalues in the $d\times n$ matrix.
Hence, the nonzero eigenvalues are roughly $c=\frac{d}{n}$ times
those following $LSD_{\frac{1}{c},\sigma^{2}}$. Several elegantly
symmetric facts follow from this. First considering the mixture structure,
the expected value of the full mixture LSD is still $\sigma^{2}$
(this comes from $\left(1-\frac{1}{c}\right)0+\frac{1}{c}(c\sigma^{2}$))
and similarly the variance remains at $c\sigma^{4}$. Second the above
formulas for support points $a=\sigma^{2}\left(1-\sqrt{c}\right)^{2}$
and $b=\sigma^{2}\left(1+\sqrt{c}\right)^{2}$ require no modification
in this new case of $c>1$. Finally the continuous density part of
the LSD has exactly the same form (\ref{eqM.I:MarcPastLSD}) scaled
by $\frac{1}{c}$. 

This theory can be used in the analysis of hierarchical\index{hierarchical clustering}
\index{clustering!asymptotic analysis}clustering in Section \ref{secK-B:HierarchicalClust}
to understand the standard deviation of PC1 of the Overlapping Classes
data\index{Overlapping Classes data} in Figure \ref{fig:HiDimToyPCA }.
As the largest principal component, the PC1 variance should be at
the upper support point of the Mar\v{c}enko-Pastur LSD. In particular,
that is $b=\sigma^{2}\left(1+\sqrt{c}\right)^{2}=1^{2}\left(1+\sqrt{\frac{1000}{100}}\right)^{2}\approx17.3$.
Taking the square root (to get back to the standard deviation) gives
4.16 which fits well with the spread of points in the upper left panel
of Figure \ref{fig:HiDimToyPCA } and the center panel of Figure \ref{fig:HiDimToyDWD}. 

A novel data based approach to choosing the number of non-noise principal
components, which is based on Mar\v{c}enko-Pastur mathematics, has
been proposed by \citet{ChoiMarron2018NumSpikes}. 

Random matrix theory in general was brought to the statistics literature
by \citet{Johnstone2001LargestEigVal}. A particularly useful tool
has been the null distribution of the largest eigenvalue, derived
by \citet{TracyWidom1994}. See \citet{ZhouMarronWirhgt2018EigenvalueTesting,ZhouMarronWright2018FamiliesIndiv}
for a deep application of the Tracy-Widom distribution in genetics.

An interesting issue is the extent to which the assumption of Gaussianity
can be weakened and still obtain the Mar\v{c}enko-Pastur distribution
for the PCA\index{Principal Component Analysis!eigenvalue distribution}
eigenvalues. Figure \ref{figM.F:MPDiffDist} explores this using scree
plots based on three data sets with $d=100$ and $n=1000$. Each data
set consists of a matrix of independent random variables all with
expected value $0$ and variance $1$. But the data sets differ in
terms of the marginal, with a Gaussian marginal scree plot shown as
black circles, an exponentially distributed plot shown as magenta
plus signs, and the Bernoulli \index{scree plot|)}scree plot in green
x signs. Note that all three are very similar, except for some natural
sampling variation. In particular the Mar\v{c}enko-Pastur LSD is not
affected by either the strong skewness\index{skewness} of the exponential
distribution or the strong discreteness of the Bernoulli. Except in
the case of heavy tails, marginal distributions matter very little,
but essentially independence of all entries of the matrix (thus in
particular independence of the variables forming the data vectors)
is crucial.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterM/OODAbookChpMFigFDiffDistMP}\caption{Scree plots indicating LDS is the same for different marginal distributions.
Gaussian is black circles, Exponential shown as magenta plus signs,
Bernoulli using smaller green x signs. \label{figM.F:MPDiffDist} }
\end{figure}
However there are important situations where classical Mar\v{c}enko-Pastur
theory does not hold, such as scale mixtures of Gaussians. The simplest
of these are of the form 
\begin{equation}
f\left(\boldsymbol{x}\right)=\omega\phi_{d}\left(\boldsymbol{x},\sigma_{1}^{2}\right)+\left(1-\omega\right)\phi_{d}\left(\boldsymbol{x},\sigma_{2}^{2}\right)\label{eqnM.D:ScaleMixNorm}
\end{equation}
where the weight $\omega\in\left(0,1\right)$ and where $\phi_{d}\left(\boldsymbol{x},\sigma^{2}\right)$
denotes the multivariate Gaussian density with expected value $\boldsymbol{0}_{d}$
and covariance matrix $\sigma^{2}\boldsymbol{I}_{d}$. Recall from
the discussion around Figure \ref{figJ.J:BatchAdjDWD} that mixtures
of Gaussians are quite natural in applications where data heterogeneity
is an important issue. Note that scale mixtures of normals also provide
good illustration of the important point that 0 covariance does not
imply independence. In particular the covariance matrix of (\ref{eqnM.D:ScaleMixNorm})
can be shown to be a multiple of the identity matrix\index{identity matrix}
(so all pairwise covariances are 0), yet entries of the vector can
be strongly dependent. E.g. when $\sigma_{1}^{2}\gg\sigma_{2}^{2}$
a large entry suggests the vector came from the first component, so
all other entries tend to be relatively large as well. \citet{LiYao2018ScaleMixRandMatrix}
have discovered an analog of the Mar\v{c}enko-Pastur distribution\index{Marchenko-Pastur distribution|)}
for more general scale mixtures of normals. In simple cases, these
can be computed using the numerical solution of integral equations
of \citet{Dobriban2015CompRandMat} (see \citet{ChoiMarron2018NumSpikes}
for a gentle introduction).

\section{High Dimension Low Sample Size\label{secM-B:HDLSS}}

As noted in the survey paper \citet{Aoshima2018HDLSSsurvey}, the
study of \emph{High Dimension Low Sample Size\index{High Dimension Low Sample Size|(}}
(HDLSS) asymptotics, where $d\rightarrow\infty$ for fixed $n$, was
popularized by \citet{HallMarronNeeman2005HDLSS}. However, earlier
work on that type of asymptotics was done in \citet{CasellaHwang1982HDLSS},
in the context of Stein estimation.

Interesting HDLSS insights arise even in the apparently straightforward
case of a single standard Gaussian random vector 
\[
\widetilde{\boldsymbol{z}}=\left(\begin{array}{c}
\widetilde{z}_{1}\\
\vdots\\
\widetilde{z}_{d}
\end{array}\right)\sim N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{I}_{d}\right).
\]
It is natural to think that observations drawn from this distribution
lie relatively near the origin, because the probability density\index{probability density}
\begin{equation}
f_{Z}\left(\widetilde{\boldsymbol{z}}\right)=\left(2\pi\right)^{-d/2}e^{-\widetilde{\boldsymbol{z}}^{t}\widetilde{\boldsymbol{z}}/2}\label{eqM.A:StandardNormalDensity}
\end{equation}
 has a peak there and decreases very rapidly in all directions as
$\left\Vert \underrightarrow{z}\right\Vert \rightarrow\infty$. However,
using the notation (\ref{eqCHp5EqA:EuclideanDistanceL2}), the Euclidean
distance\index{distance!Euclidean} from $\widetilde{\boldsymbol{z}}$
to the origin is 
\[
\delta_{2}\left(\widetilde{\boldsymbol{z}},\boldsymbol{0}_{d,1}\right)=\left(\sum_{j=1}^{d}\widetilde{z}_{j}^{2}\right)^{1/2},
\]
whose distribution is easy to understand, since $\sum_{j=1}^{d}\widetilde{z}_{j}^{2}$
has a Chi Squared distribution with $d$ degrees of freedom. That
distribution has expected value $d$ and variance $2d$ (i.e. standard
deviation $\sqrt{2d}$). A simple Taylor series (sometimes called
the \emph{delta method}\index{delta method}) calculation shows that
in the limit as $d\rightarrow\infty$ (for fixed $n$)
\begin{equation}
\delta_{2}\left(\widetilde{\boldsymbol{z}},\boldsymbol{0}_{d,1}\right)=\sqrt{d}+O_{p}\left(1\right).\label{eqM.B:NormDistTo0}
\end{equation}
This shows that most of the distribution of $\widetilde{\boldsymbol{z}}$
lies near the surface of a growing sphere of radius $\sqrt{d}$ in
$\mathbb{R}^{d}$, and perhaps paradoxically there is essentially
no chance of finding $\widetilde{\boldsymbol{z}}$ anywhere near the
origin, $\boldsymbol{0}_{d,1}$. This apparent contradiction to the
density (\ref{eqM.A:StandardNormalDensity}) being much larger near
the origin is resolved by recalling that this is density \emph{with
respect to Lebesgue measure}. To see how Lebesgue measure can give
perhaps non-intuitive results such as (\ref{eqM.B:NormDistTo0}) consider
the calculation of the volume of the (solid) unit ball in $\mathbb{R}^{d}$.
This is most easily calculated as an integral using polar coordinates.
Writing the integral with respect to radius last, that integrand goes
up extremely rapidly near the end of the interval $[0,1]$ for large
$d$. This shows that high dimensional Lebesgue measure tends to strongly
push probability mass outwards. The balance point between this effect
and the rapidly diminishing tails of the Gaussian density (\ref{eqM.A:StandardNormalDensity})
to keep mass near the origin is seen to be at $\sqrt{d}$ in (\ref{eqM.B:NormDistTo0}).
This has interesting implications for individuals in that when enough
measurements are made (i.e. $d$ is sufficiently large) it is essentially
impossible for any individual to be ``average'', i.e. near the population
expected value vector. The result (\ref{eqM.B:NormDistTo0}) also
explains the perhaps surprising observation in the discussion of Figure
\ref{figK.G:HierHiDSingView} that two standard Gaussian data points
are so far from the rest of the data. Recall that $d=500$ and $n=30$
for that High Dimensional Gaussian data\index{High Dimensional Gaussian data|(}
set. Thus $\sqrt{500}\approx22.4$ well explains the \index{distance|(}distances
from the two ``outlying'' points to the origin in the first diagonal
panels of Figure \ref{figK.G:HierHiDSingView}.

A result in the same spirit as (\ref{eqM.B:NormDistTo0}) holds for
independent $\widetilde{\boldsymbol{z}}_{1},\widetilde{\boldsymbol{z}}_{2}\sim N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{I}_{d}\right)$,
\begin{equation}
\delta_{2}\left(\widetilde{\boldsymbol{z}}_{1},\widetilde{\boldsymbol{z}}_{2}\right)=\sqrt{2d}+O_{p}\left(1\right).\label{eqM.C:PairNormDist}
\end{equation}
This is not surprising since the difference of independent standard
normals is again normal with standard deviation $\sqrt{2}$. Applying
this idea to a sample $\widetilde{\boldsymbol{z}}_{1},\cdots,\widetilde{\boldsymbol{z}}_{n}$
results in the idea that all of these points are equidistant from
each other. Thus for $n=3,$ they lie at vertices of an equilateral
triangle (a similar phenomenon was discovered by \citet{Kendall1988Triangles}
in the context of shape diffusion), and for $n=4$ at vertices of
a regular tetrahedron. Even conceptualizing this for $n>4$ is challenging
because the point configuration lies in a space with dimension higher
than 3. Humans tend to be good at perception in three dimensions,
but quite poor at higher dimensional perception. This may be because
our perceptual systems come from our ancestors for whom a major task
was finding food. Food exists in three dimensions which seems to be
why we are so good at understanding three dimensional space.

One more related result in this spirit comes from studying the angle
(with vertex at the origin) between independent vectors $\widetilde{\boldsymbol{z}}_{1},\widetilde{\boldsymbol{z}}_{2}\sim N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{I}_{d}\right)$,
\begin{equation}
\angle\left(\widetilde{\boldsymbol{z}}_{1},\widetilde{\boldsymbol{z}}_{2}\right)=90^{\circ}+O_{p}\left(d^{-1/2}\right).\label{eqM.D:NormAngles}
\end{equation}
Note that this right angle completes the Pythagorean relationship
in terms of the right triangle with both leg lengths $\sqrt{d}$ from
(\ref{eqM.B:NormDistTo0}) and hypotenuse $\sqrt{2d}$ from (\ref{eqM.C:PairNormDist}).
The orthogonality of data directions observed in the High Dimensional
Gaussian data in the top center panel of Figure \ref{figK.G:HierHiDSingView}
is also explained by (\ref{eqM.D:NormAngles}).

Equations (\ref{eqM.B:NormDistTo0}), (\ref{eqM.C:PairNormDist})
and (\ref{eqM.D:NormAngles}) have some interesting implications for
the very nature of high dimensional variation. While letting $d\rightarrow\infty$
(for a fixed set of $n$ data objects) clearly adds more and more
noise to a Gaussian data set, that noise manifests itself in a rather
special way. In particular, given a sample $\widetilde{\boldsymbol{z}}_{1},\cdots,\widetilde{\boldsymbol{z}}_{n}$
consider the subspace of dimension $n$ generated by the data. This
is a hyperplane of dimension $n$ passing through the origin. This
plane can be rotated in $\mathbb{R}^{d}$ to the subspace of the first
$n$ coordinates. Within that subspace, all data points are essentially
distance $\sqrt{d}$ from the origin, and $\sqrt{2d}$ from each other
and can be again rotated so that each data point is close to a coordinate
axis. This results in the data lying near the vertices of the $\sqrt{d}$
rescaling of the $n$-dimensional \emph{unit simplex} in $\mathbb{R}^{d}$.
Hence while increasing $d$ creates increasing variation, that goes
mostly into very random rotation. Modulo that rotation and scaling,
the data actually converge to a rigid deterministic structure. This
HDLSS phenomenon has been called \emph{geometric representation\index{geometric representation|(}}
by \citet{HallMarronNeeman2005HDLSS}.

A simple illustration of geometric representation is shown in Figure
\ref{figM.B:ViewGeomRep}. The top left panel shows a conceptual model
of $n=2$ points in $\mathbb{R}^{3}$ as red circles. The top right
panel shows the subspace generated by these points (i.e. the 2-d plane
containing them and the origin) in cyan. Red solid lines are also
used to depict the concept from (\ref{eqM.B:NormDistTo0}) that each
point is approximately distance $\sqrt{d}$ from the origin. Similarly
a dashed line illustrates the roughly $\sqrt{2d}$ distance between
points from (\ref{eqM.C:PairNormDist}). The first step of geometric
representation is to rotate this cyan subspace into the coordinate
system of the first $n=2$ axes in $\mathbb{R}^{3}$ as shown in the
lower left panel of Figure \ref{figM.B:ViewGeomRep}. The second rotation
happens within this $\mathbb{R}^{2}$, to move the points close to
the axes, resulting in the unit simplex (in this case just a line
segment) approximation. 

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterM/OODAbookChpMFigBViewGeomRep}

\caption{Illustration of HDLSS geometric representation using just two data
points. Shows how added noise from increasing dimensionality goes
mostly into random rotation, while modulo rotation data points (red
circles) converge to vertices of a rescaled unit simplex (black circles).
\label{figM.B:ViewGeomRep}}
\end{figure}
Another illustration of geometric representation is shown using a
small simulation study in Figure \ref{eqM.C:PairNormDist}. Each panel
shows seven (the number of colors in Matlab's default color palette)
$N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{I}_{d}\right)$ data
sets of size $n=3$, in the dimensions $d$ shown above each panel.
Each data set determines a hyperplane of dimension 2, which is rotated
into the $\mathbb{R}^{2}$ shown in each panel. Next within that $\mathbb{R}^{2}$
each data set is rotated so the first two data points are equidistant
from the vertical axis, with the minimum of those two below the origin.
These two operations result in modding out the increasing amount of
random rotation, so the points coalesce to vertices of the 3 dimensional
unit simplex, which in this view are the vertices of the shown equilateral
triangle in $\mathbb{R}^{2}$ with side length $\sqrt{2d}$. For $d=2$
shown in the upper left panel, the points appear to be rather random.
Yet at only $d=20$ in the upper right panel the points are already
beginning to coalesce. In the lower left panel for $d=200$ the points
are now quite close to the vertices of the triangle. The much larger
$d=20,000$ is used in the lower right panel to show how the points
are now so close to the vertices that the seven different data sets
are barely distinguishable due to over-plotting.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterM/OODAbookChpMFigCGeomRepSim}

\caption{Simulations demonstrating geometric representation for several dimensions
$d$, each with seven data sets (shown in different colors) of size
$n=3$. Shows for increasing dimension that modulo rotation and rescaling
data points converge to vertices of a rescaled unit simplex. \label{figM.C:GeomRepSim}}
\end{figure}
Geometric representation calculations can explain a number of the
apparently strange phenomena observed earlier. For example in the
discussion of the Overlapping Classes data\index{Overlapping Classes data}
analysis in Figure \ref{fig:HiDimToyDWD}, it was stated (and verified
by the DiProPerm test) that the perhaps surprisingly distinct DWD
difference observed between the two classes (which were generated
from the same probability distribution) in the top left panel was
simply the result of natural variation in Gaussian data. The first
observation is that the $n_{1}$ red points determine a hyperplane
of dimension $(n_{1}-1)$, and similarly for the $n_{2}$ blue points.
The pairwise distance relationship (\ref{eqM.C:PairNormDist}) indicates
that all pairs consisting of one red and one blue point are approximately
$\sqrt{2d}$ apart. This shows that the red and blue hyperplanes are
essentially parallel. Applying DWD to such data will result in a direction
vector that is essentially in the direction of the common normal vectors
of the planes, explaining the two very distinct clusters in the top
left panel of Figure \ref{fig:HiDimToyDWD}. Furthermore the distance
between the clusters can be approximated by the distance between the
hyperplanes. The first step of that computation is to find the distance
from the origin to each plane. It follows from (\ref{eqM.B:NormDistTo0})
that each point is distance $\sqrt{d}$ from the origin. To understand
how this can be used to calculate the distance from the origin to
the plane, rotate say the red plane so that the data points essentially
lie in the first $n_{1}$ positive coordinate directions (which is
possible because of the approximate orthogonality result (\ref{eqM.D:NormAngles})).
Thus the red data points can be thought of as lying near the vertices
of the $\sqrt{d}$ rescaled, $n_{1}$ dimensional unit simplex 
\[
\left\{ \left(\begin{array}{cccccc}
x_{1} & \cdots & x_{n_{1}} & 0 & \cdots & 0\end{array}\right)^{t}\in\mathbb{R}^{d}:\sum_{j=1}^{n_{1}}x_{j}=\sqrt{d},\;x_{j}\geq0\right\} .
\]
A good approximation of the mean of the red points can be found from
the value of $c$ for which the vector $\left(\begin{array}{cccccc}
c & \cdots & c & 0 & \cdots & 0\end{array}\right)^{t}$ lies on the rescaled simplex, which entails $\sum_{j=1}^{n_{1}}c=n_{1}c=\sqrt{d}$,
i.e. $c=\frac{\sqrt{d}}{n_{1}}$. Note that the distance from this
approximate mean to the origin is $\sqrt{\sum_{j=1}^{n_{1}}\left(\frac{\sqrt{d}}{n_{1}}\right)^{2}}=\sqrt{\frac{d}{n_{1}}}$.
A similar idea applies to the approximate simplex of $n_{2}$ blue
points in Figure \ref{fig:HiDimToyDWD}. In particular the distance
from the mean of the blue points to the origin is approximately $\sqrt{\frac{d}{n_{2}}}$.
Note that this is very consistent with one of the key concepts from
the Central Limit Theorem\index{Central Limit Theorem} that averaging
$n$ independent random variables reduces the standard deviation by
a factor of $n^{-1/2}$. Now the direction separating the red and
blue clusters in the upper right panel of Figure \ref{fig:HiDimToyDWD}
is essentially driven by the difference between these two means. Since
the points in the two clusters are all orthogonal, by (\ref{eqM.D:NormAngles}),
their respective subspaces are orthogonal as well. It follows that
the two mean vectors are orthogonal to each other, and thus the distance
between them can be computed using the Pythagorean theorem\index{Pythagorean theorem},
$\sqrt{\sqrt{\frac{d}{n_{1}}}^{2}+\sqrt{\frac{d}{n_{2}}}^{2}}=\sqrt{\frac{d}{n_{1}}+\frac{d}{n_{2}}}$.
Plugging in the values $d=1000$ and $n_{1}=n_{2}=50$ results in
approximate distance between the means of $\sqrt{\frac{d}{n_{1}}+\frac{d}{n_{2}}}\approx6.3$.
This fits well with the observed cluster means shown in the upper
left panel of Figure \ref{fig:HiDimToyDWD}. A related HDLSS analysis
of clustering\index{clustering!asymptotic analysis}, as defined in
Chapter \ref{chapK:Clustering}, can be found in \citet{Borysov2014Clusters}.

A parallel calculation similarly predicts the apparently very large
distances between the colored clusters, for the $d=500$, $n=30$
High Dimensional Gaussian data\index{High Dimensional Gaussian data|)}
shown in Figure \ref{figK.H:HierHiDSWardView}. For the cyan - magenta
clusters the mean difference (left column) prediction is $\sqrt{\frac{500}{8}+\frac{500}{5}}\approx12.7$,
while for the green - red (center column) $\sqrt{\frac{500}{7}+\frac{500}{10}}\approx11.0$. 

A similar explanation of the behavior observed (for MD using a simple
Gaussian toy data set) in the top row of Figure \ref{figJ.C:HighDClass}
is a little harder because the underlying theoretical means are not
0, but instead are 4.4 units apart in the first coordinate direction.
Again reflecting the approximate orthogonality of the first coordinate
with both means the Pythagorean relationship is $\sqrt{\frac{d}{20}+4.4^{2}+\frac{d}{20}}$.
Plugging in $d=10$ (top left panel of Figure \ref{figJ.C:HighDClass})
gives about $4.5$, while $d=40$ (top center) yields around $4.8$
and $d=200$ (top right) gives about $6.3$. An interesting open problem
is similar explanation of the differences (for the more complicated
classifiers LDA and MDP\index{Maximal Data Piling}) in the other
panels of Figure \ref{figJ.C:HighDClass}. 

The above calculations all assumed a standard Gaussian $N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{I}_{d}\right)$
distribution for simplicity of illustration. \citet{HallMarronNeeman2005HDLSS}
extended this to non-Gaussian distributions by allowing dependence
in the form of a mixing condition along the data vector. Mixing conditions
are a natural way to allow some dependence in the laws of large numbers
that lie at the heart of geometric representation. See \citet{Bradley2005MixCond}
and \citet{White2014Prob4Econ} for good overviews of the large literature
on mixing conditions. The first submission of \citet{HallMarronNeeman2005HDLSS}
was rejected by a top statistical journal on the grounds that the
reviewers were unaware of applications where such a dependence assumption
was sensible. That later turned out to be an ironic decision, as that
same year of 2005 was the date of the first publication on Genome
Wide Association Studies (GWAS) by \citet{Klein2005FirstGWAS}. As
discussed in \citet{ZhouMarronWirhgt2018EigenvalueTesting,ZhouMarronWright2018FamiliesIndiv}
and in Section \ref{secO.B.B:GWAS} GWAS is an excellent example of
a data analysis context where geometric representation ideas have
given clear insights, and the mixing condition assumption of \citet{HallMarronNeeman2005HDLSS}
is a very natural consequence of the splitting and merging of chromosomes
that constitutes the reproductive process.

A number of other assumptions beyond normality for geometric representation
have been developed by \citet{AhmMarronMuller2007HDLSS} and \citet{YataAoshima2012GeomRep}.
It was pointed out by John Kent that when exploring sufficient conditions
for geometric representation a useful example to keep in mind is the
scale mixture of Gaussians in (\ref{eqnM.D:ScaleMixNorm}). In that
case the sharp $\sqrt{d}$ in (\ref{eqM.B:NormDistTo0}) becomes just
$O_{p}\left(\sqrt{d}\right)$ because with probability $\omega$ the
observation is $\sigma_{1}\sqrt{d}$ from the origin, and at distance\index{distance|)}
$\sigma_{2}\sqrt{d}$ otherwise. However it is important to note that
geometric representation is more a property of the \emph{dependence
structure} than of the marginal distribution (consistent with the
ideas of Figure \ref{figM.F:MPDiffDist}). As noted in the discussion
of (\ref{eqL.M:IsoVar}), an important consequence of assuming independent
variables is the joint distribution has common variance of projections
in all directions. Another consequence is that in the case of GWAS
discussed above, the marginal distributions are discrete (often binary),
yet geometric representation\index{geometric representation|)} ideas
still give useful insights.

Because visualization\index{visualization} of PCA\index{Principal Component Analysis!HDLSS|(}
\index{scores|(}scores has been so useful at many points above, the
HDLSS mathematical underpinnings of these are considered here. In
view of the tendency of angles to converge to $90^{\circ}$ (as seen
in (\ref{eqM.D:NormAngles})) a natural starting point is the behavior
of the direction vectors used for projection in PCA, i.e. the eigenvectors
of the sample covariance matrix as detailed in (\ref{eqP-O:EigenAnalofCov}).
For this assume the random sample\index{random sample} $\widetilde{\boldsymbol{x}}_{1},\cdots,\widetilde{\boldsymbol{x}}_{n}$
is drawn (independently) from the Gaussian distribution $N_{d}\left(\boldsymbol{0}_{d,1},\boldsymbol{\varSigma}_{d}\right)$,
where the $d$ dimensional theoretical covariance matrix $\boldsymbol{\varSigma}_{d}$
has eigenvalues $\lambda_{1,d},\cdots,\lambda_{d,d}$ and first eigenvector
$\boldsymbol{u}_{1,d}\in\mathbb{R}^{d}$. Further assume a \emph{\index{spike covariance model}spike
covariance model}, in the spirit of \citet{Paul2007SpikeCov}, where
most eigenvalues are fixed, but the first eigenvalue grows with the
dimension. In particular assume 
\begin{equation}
\lambda_{1,d}=d^{\alpha},\;\lambda_{2,d}=\cdots=\lambda_{d,d}=1,\label{eqnM.E:SpikeCovAssumption}
\end{equation}
where $\alpha>0$ controls the growth rate of the first eigenvalue.
This distribution has a hot dog shape, with a single long axis in
the direction $\boldsymbol{u}_{1,d}$, and smaller spherically symmetric
variation in all orthogonal directions. Next let $\widehat{\boldsymbol{\varSigma}}_{d}$
denote the sample covariance matrix defined at (\ref{eqB.F:SampleCovariance}),
which is a natural estimate of $\boldsymbol{\varSigma}_{d}$ and similarly
define estimators $\widehat{\boldsymbol{u}}_{1,d},\hat{\lambda}_{1,d},\cdots,\hat{\lambda}_{d,d}$
using the corresponding sample eigenanalysis of $\widehat{\boldsymbol{\varSigma}}_{d}$
(again see (\ref{eqP-O:EigenAnalofCov})). Because $\boldsymbol{u}_{1,d}$
and $\widehat{\boldsymbol{u}}_{1,d}$ are unit vectors, it is natural
to study their relationship using angles. \citet{Jung2009Marron}
showed that the power $\alpha$ is critical to this angle. In particular,
in the limit as $d\rightarrow\infty$ with fixed $n$,
\begin{equation}
\angle\left(\widehat{\boldsymbol{u}}_{1,d},\boldsymbol{u}_{1,d}\right)\rightarrow\begin{cases}
\begin{array}{c}
0^{\circ}\\
90^{\circ}
\end{array} & \begin{array}{c}
\alpha>1\\
\alpha<1
\end{array}.\end{cases}\label{eqM.F:ConsistStrongInconsist}
\end{equation}
The first part of (\ref{eqM.F:ConsistStrongInconsist}) has the form
of standard \emph{consistency\index{consistency|(}}, in the sense
that when the first eigenvalue is large enough the empirical eigenvector
$\widehat{\boldsymbol{u}}_{1,d}$ appropriately converges to the underlying
true eigenvector $\boldsymbol{u}_{1,d}$. The second part of (\ref{eqM.F:ConsistStrongInconsist})
is much less standard in mathematical statistics. It says when the
eigenvalue is too small, not only do these directions not converge,
but in fact they diverge to the greatest extent possible (recall eigenvectors
are determined only up to $\pm$ sign flips, so angles $>90^{\circ}$
should be replaced by their supplements). Hence this condition is
called \emph{strong inconsistency}\index{inconsistency!strong|(}.
 To understand why $\alpha=1$ is critical here, recall from (\ref{eqM.B:NormDistTo0})
that Standard Gaussian data tend to lie near the surface of the unit
sphere with radius $\sqrt{d}$. When $\alpha>1$ the stretched Gaussian
distribution has standard deviation $d^{\alpha/2}\gg\sqrt{d}$, where
the exponent is $\alpha/2$ because eigenvalues are on the scale of
variance. Thus most of the probability mass of the stretched Gaussian
is outside the sphere of radius $\sqrt{d}$ pulling the empirical
eigenvector $\widehat{\boldsymbol{u}}_{1,d}$ in the direction of
$\boldsymbol{u}_{1,d}$. Furthermore when $\alpha<1$, it follows
that $d^{\alpha/2}\ll\sqrt{d}$ so essentially all of the mass due
to the spike lies within the standard sphere.  In this case $\widehat{\boldsymbol{u}}_{1,d}$
will take some random direction, but according to (\ref{eqM.D:NormAngles})
random directions tend to be orthogonal, resulting in the $90^{\circ}$
part of (\ref{eqM.F:ConsistStrongInconsist}). Important extensions
of these results, including to more than one spike eigenvalue and
also to a broader range of asymptotics can be found in \citet{Shen2016GenFrame,Shen2016HDLSS}. 

Note that (\ref{eqM.F:ConsistStrongInconsist}) leaves open the boundary
case of $\alpha=1$. That was resolved in \citet{JungSenMarron2012Boundary}
with the establishment of appropriate limit distributions that appropriately
meld the two results in (\ref{eqM.F:ConsistStrongInconsist}). That
set of mathematics provides the basis for an important proposal for
selection of the number of eigenvalues in PCA by \citet{JungLeeAhn2018HDLSSRankSel}.

\citet{Jung2009Marron} also studied the behavior of the first eigenvalue
and showed that under the eigenvalue assumption (\ref{eqnM.E:SpikeCovAssumption}),
again in the limit as $d\rightarrow\infty$, 
\begin{equation}
\frac{\hat{\lambda}_{1,d}}{\lambda_{1,d}}\stackrel{d}{\rightarrow}\frac{\chi_{n}^{2}}{n}.\label{eqM.G:InconsistEigenvalue}
\end{equation}
Thus the first spike sample eigenvalue is inconsistent, although in
a rather controlled and interpretable way. For example $E\chi_{n}^{2}=n$,
so $\hat{\lambda}_{1,d}$ is unbiased in a crude sense. This issue,
together with a broad set of different conditions for consistency
and strong inconsistency along the lines of (\ref{eqM.F:ConsistStrongInconsist})
has been explored in a series of paper by \citet{YataAoshima2009PCAconsist,YataAoshima2010aEffectivePCA,YataAoshima2010bIntrinsicDim,YataAoshima2012GeomRep,YataAoshima2013PowerSpike}.
Some landmarks of those results are proposed variations of PCA that
features a consistent first eigenvalue, thus overcoming the perhaps
disappointing result (\ref{eqM.G:InconsistEigenvalue}), as well as
obtaining consistency for some values of $\alpha<1$ in (\ref{eqnM.E:SpikeCovAssumption}).
A less elegant way to overcome the inconsistency of (\ref{eqM.G:InconsistEigenvalue})
is to follow the $\underset{d\rightarrow\infty}{\lim}$ with a $\underset{n\rightarrow\infty}{\lim}$
as discussed in Section \ref{secM-C:HDMSS}. 

Given the popularity of \emph{sparse} methods\index{sparse methods}
(those that approach high dimensional challenges by assuming large
numbers of 0s) in the world of mathematical statistics, it is natural
to investigate parallels of (\ref{eqM.F:ConsistStrongInconsist})
in sparse contexts. \citet{Shen2013SparsePCA} developed a complementary
theory based on a sparsity\index{sparsity} \index{sparsity index}index
$\beta$ (in addition to the signal strength parameter $\alpha$),
and laid out interesting regions for consistency and strong \index{inconsistency!strong|)}inconsistency
in the $\alpha,\beta$ parameter space. 

Results of the form (\ref{eqM.F:ConsistStrongInconsist}) and (\ref{eqM.G:InconsistEigenvalue})
have generated some considerable skepticism. An interesting objection
is that (\ref{eqM.F:ConsistStrongInconsist}) holds even for sample
size $n=1$. This may be surprising, however in that case the direction
$\widehat{\boldsymbol{u}}_{1,d}$ points in the direction of the single
data vector $\widetilde{\boldsymbol{x}}_{1}$. The stretch of the
spiked Gaussian density described above in the direction $\boldsymbol{u}_{1,d}$
then explains why (\ref{eqM.F:ConsistStrongInconsist}) holds. Discussing
consistency in the case $n=1$ is certainly not typical, and one might
interpret this as an indication that the spike assumption (\ref{eqnM.E:SpikeCovAssumption})
is too strong to have any practical relevance. However the two sides
of (\ref{eqM.F:ConsistStrongInconsist}) reveal that when no assumption
in the spirit of $\alpha>1$ holds, PCA cannot give any useful results.
Yet at many points above, e.g. Figures \ref{fig:LungCancerLogCurves}
- \ref{fig:LungCancerLogCurvesBrushed}, PCA in very high dimensions
($d=1709$ in that case) clearly gives nonrandom, useful results.
Hence while (\ref{eqnM.E:SpikeCovAssumption}) may appear strong to
those theoreticians whose thoughts are unencumbered by real data,
in fact it evidently is reasonable in interesting practical contexts. 

Another potential objection to this theory is that (\ref{eqM.G:InconsistEigenvalue})
hints at potential inconsistency of PC scores in high dimensions.
To investigate this, for $i=1,\cdots,n$ let $\widehat{s}_{i,j}$
denote the $i$th PC score for the $j$th component, i.e. the projection
of $\widetilde{\boldsymbol{x}}_{i}$ onto the eigenvector $\widehat{\boldsymbol{u}}_{1,d}$.
Now $\widehat{s}_{i,j}$ can be viewed as an estimate of the theoretical
score $s_{i,j}$, which is the projection of $\widetilde{\boldsymbol{x}}_{i}$
onto the theoretical eigenvector $\boldsymbol{u}_{1,d}$. Using results
from \citet{Shen2016HDLSS}, it can be shown that under the appropriate
analogs of the eigenvalue assumption (\ref{eqnM.E:SpikeCovAssumption})
\begin{equation}
\frac{\widehat{s}_{i,j}}{s_{i,j}}\stackrel{d}{\rightarrow}R_{j},\label{eqM.H:InconsistScores}
\end{equation}
as $d\rightarrow\infty$, where $R_{j}$ is a non-degenerate random
variable. Consistency of PC scores would correspond to a degenerate
$R_{j}=1$, thus revealing general \index{inconsistency|)}inconsistency
of PC scores in HDLSS settings. This seems like an apparent contradiction
to the clear usefulness of PCA in finding interesting structure in
data, e.g. as in Figures \ref{fig:LungCancerPCA} and \ref{fig:LungCancerPCABrushed}.
That is resolved by noting that while the limiting random variable
$R_{j}$ depends on the component index $j$, it does not depend on
the data index $i=1,\cdots,n$. Thus for each data point, it has the
same realization, so while the scores are indeed off by a random factor
of $R_{j}$, it is the \emph{same factor} for each $\widetilde{\boldsymbol{x}}_{i}$.
Hence in PCA scatterplot matrices\index{scores!scatterplot} in high
dimensional situations, such as Figures \ref{fig:LungCancerPCA} and
\ref{fig:LungCancerPCABrushed}, the numbers on the axes do not reflect
the underlying probability distribution. However the \emph{relative
positioning} of the data points is correct which is why PCA is very
useful in high dimensional data analysis. In addition to the results
mentioned above, \citet{YataAoshima2010aEffectivePCA,YataAoshima2012GeomRep}
have proposed modifications of PCA which \emph{are} consistent in
the sense that $R_{j}=1$. 

Another point of skepticism that has been raised is that it is unclear
how to check assumptions such as (\ref{eqnM.E:SpikeCovAssumption})
``in practice''. But this loses sight of the fact that these results
are merely mathematical statistics, and such issues exist almost everywhere
in that area. For example, the Central Limit Theorem\index{Central Limit Theorem}
also relies on assumptions such as independence and second moments
which similarly cannot be proven from any data set, yet its value
in explaining naturally observed phenomena is unquestioned. This same
principle applies as well to the asymptotics described here, as indeed
it applies to all of mathematical statistics. 

A particularly interesting HDLSS characterization of the usefulness
of PCA\index{Principal Component Analysis!HDLSS|)} scores can be
found in \citet{HelltonThoresen2017PCA}, which explores the impact
of measurement error on HDLSS scores\index{scores|)}. Interesting
work generalizing rank based nonparametric\index{nonparametric} statistics
to HDLSS contexts can be found in \citet{Biswas2014TwoSample,Biswas20165OneSample,Biswas2016TestIndep}
and \citet{Ghosh2016Hyperplanes}. 

As noted in \citet{Aoshima2018HDLSSsurvey}, the mathematics behind
many of these results can be established using the theory called \emph{concentration
of measure}\index{concentration of measure}, which has a large literature
started by \citet{Talagrand1991ConcMeas,Talagrand1995ConcMeas}. See
\citet{Ledoux2001ConcMeas} for good overview, and \citet{Kolchinskii2017ConcMeas,Koltchinskii2016ConcMeas,Koltchinskii2017PCA}
for useful directly related results. This is the theory that seems
to have generated the common misunderstanding, typified by the quote
from \citet{Wainwright2019HiDimStat} at the beginning of this chapter.

Closely related asymptotics also based on fixed sample size $n$,
but following the tangential path of hypothesis testing for repeated
measures, have have been explored by \citet{Pesarin2010HDLSSRepMeas},
who propose a notion called \emph{finite sample consistency\index{consistency|)}}.

\section{High Dimension Medium Sample Size\label{secM-C:HDMSS} }

An asymptotic domain that lives between the random matrix theory\index{random matrix theory|)}
discussed in Section \ref{secM-A:RandomMatrix} and the HDLSS context
of Section \ref{secM-B:HDLSS} is where both $d,n\rightarrow\infty$
with $d\gg n$. This has been called \emph{High Dimension Medium Sample
Size\index{High Dimension Medium Sample Size}} (HDMSS) by \citet{YataAoshima2013PowerSpike}
and \emph{ultra high dimensionality\index{ultra high dimensionality}}
by \citet{FanLv2008UltraHiDim}. The added $\underset{n\rightarrow\infty}{\lim}$
results in a theory that corresponds somewhat more with conventional
mathematical statistics. In particular some of the inconsistency results
such as (\ref{eqM.G:InconsistEigenvalue}) and (\ref{eqM.H:InconsistScores})
turn into more conventional consistency statements. There is still
a sharp boundary between consistency and strong inconsistency of the
sample eigenvectors of the type in (\ref{eqM.F:ConsistStrongInconsist}),
but as clearly expressed in Figure 1 of \citet{Shen2016GenFrame},
the $\alpha=1$ threshold that was critical to the boundary is now
replaced by one driven by the relative rates at which $d$ and $n$
tend to $\infty$.

A conceptual issue is the parameterization of the asymptotics in the
HDMSS domain. \citet{FanLv2008UltraHiDim} took a conventional approach
working with the limit as $n\rightarrow\infty$ and taking for example
$d\sim n^{\theta}$. On the other hand, \citet{YataAoshima2013PowerSpike}
base their analysis on the limit as $d\rightarrow\infty$ and working
with $n=d^{\rho}$. Of course they are essentially the same domain
but note that the Yata-Aoshima formulation gives a much more natural
interface between the HDMSS and HDLSS\index{High Dimension Low Sample Size|)}
worlds.
