
\chapter{Classification - Supervised Learning\label{chapJ:Classification}}

Statistical \emph{classification\index{classification|(}}, also called
\emph{discrimination\index{discrimination}} is easily understood
in the context of data-driven automatic disease diagnosis. The process
is based on a \emph{\index{training data|(}training set} consisting
of patients with known labels such as the presence or absence of the
disease under consideration. The training data also includes a \index{trait}trait
(i.e. feature\index{feature}) vector of relevant data for each patient.
The goal of classification is a data-based rule for assigning new
cases to either the disease present or absent classes on the basis
of their data vector. 

As noted in Section \ref{secB4:Further-Statistical-Tasks} and even
this chapter title, classification is also called \emph{supervised
learning\index{supervised learning}} in machine learning\index{machine learning},
see \citet{James2013RStatLearn} and \citet{HastiTibsFried2009ElStatLearn},
where the given class labels constitute the supervision. This provides
a useful contrast with \emph{\index{unsupervised learning}unsupervised
learning}, which is called by the more classical name of \index{cluster analysis}cluster
analysis in Chapter \ref{chapK:Clustering}. Both focus on subsets
of the data, that are considered known in the former, and to be discovered
in the latter.

Early work in this area, e.g. \citet{Fisher1936LDA}, had strong connections
with biology, where the term ``classification'' refers to development
and refinement of the taxonomic organizational structure of kingdom-phylum-...-genus-species
that has been at the heart of that field for many years. Hence as
the term ``classification'' was already used, discrimination was
once preferred terminology for the statistical task studied here.
However, that term has since become so socially charged that classification
has become more commonly used in statistical contexts. This may be
connected to ``discrimination'' having been coined by advocates
of the widely discredited area of eugenics. Another name, popular
in the machine learning literature, is \emph{pattern recognition\index{pattern recognition}}.

The general classification problem involves $K$ classes, but for
ease of notation typically the \emph{\index{binary classification|(}binary
classification} case of $K=2$ will usually be considered here. There
are many notational systems in use for discussing classification,
but a particularly useful one (especially for linear classifier\index{linear classifier|(}s)
is to assume \emph{training data} 
\begin{equation}
\left(\widetilde{\boldsymbol{x}}_{1},\widetilde{y}_{1}\right),\cdots,\left(\widetilde{\boldsymbol{x}}_{n},\widetilde{y}_{n}\right),\label{eqJ.H:TrainingData}
\end{equation}
where $\widetilde{\boldsymbol{x}}_{i}\in\mathbb{R}^{d}$ and $\widetilde{y}_{i}\in\left\{ -1,1\right\} $
for $i=1,\cdots,n$. As discussed around Table \ref{tab7-1:TildeNotation}
the tilde on top of scalars and vectors are used to indicate that
these are random quantities. . Let $\left(\widetilde{\boldsymbol{x}}_{0},\widetilde{y}_{0}\right)$
denote a separate test case. The goal is to find a \emph{classification
rule\index{classification rule|(}} for using $\widetilde{\boldsymbol{x}}_{0}$
to predict the corresponding class label $\widetilde{y}_{0}$, i.e.
a function $c\left(\boldsymbol{x}\right):\mathbb{R}^{d}\rightarrow\left\{ -1,1\right\} $.

While not everyone in the machine learning literature does this, it
is common to assume that the training data are drawn independently
from an underlying joint probability distribution $f\left(\boldsymbol{x},y\right)$.
It is useful to think in a Bayesian way, starting with a loss function
of $0$ for a correct classification and $1$ for incorrect. Next
define \emph{\index{prior probabilities}prior probabilities} to be
the marginals $p_{+}=P\left\{ \widetilde{y}=+1\right\} $ and $p_{-}=P\left\{ \widetilde{y}=-1\right\} =1-p_{+}$.
Performance of a classifier, $c$, is assessed in a natural way using
the \index{Bayes risk}Bayes risk (expected loss), $P\left\{ c(\widetilde{\boldsymbol{x}})=\widetilde{y}\right\} $.
The \index{likelihood ratio rule}likelihood ratio rule, 
\[
c_{LR}\left(\boldsymbol{x}\right)=\left\{ \begin{array}{cc}
+1, & \mathrm{when}\;\frac{f\left(\boldsymbol{x},+1\right)}{f\left(\boldsymbol{x},-1\right)}\geq1\\
-1, & \mathrm{when}\;\frac{f\left(\boldsymbol{x},+1\right)}{f\left(\boldsymbol{x},-1\right)}<1
\end{array}\right.
\]
is Bayes risk optimal in the sense of minimizing the Bayes risk over
all choices of $c$. Of course $f$ is typically unknown, but the
Bayes rule still provides an ideal to keep in mind and its Bayes risk
is an upper bound on the classification performance of all possible
methods. 

An important general class of methods are the \emph{linear classifiers}.
Given a \emph{\index{classification direction vector}direction vector}
$\boldsymbol{w}\in\mathbb{R}^{d}$ (i. e. $\left\Vert \boldsymbol{w}\right\Vert =1$)
and an intercept $\beta\in\mathbb{R}$, define 
\begin{equation}
c_{\boldsymbol{w},\beta}\left(\boldsymbol{x}\right)=1-2\cdot1_{\left\{ \boldsymbol{w}^{t}\boldsymbol{x}<\beta\right\} },\label{eqJ.A:LinearClasser}
\end{equation}
where $1_{\left\{ \cdot\right\} }$ denotes the \emph{\index{indicator function}indicator
function} (in contrast to the boldface $\boldsymbol{1}_{d,n}$ notation
for the matrix of ones in (\ref{eq:ones matrix})) that is 1 when
the condition $\left\{ \cdot\right\} $ is true, and is 0 otherwise.
Note that $c_{\boldsymbol{w},\beta}$ assigns $\boldsymbol{x}$ to
$+1$ exactly when it is on the positive (with respect to the direction
of $\boldsymbol{w}$) side of the hyperplane whose normal vector is
$\boldsymbol{w}$ and intercept is $\beta$. I.e. there is a \emph{\index{separating surface|(}separating
surface} between the two class regions that is a hyperplane.

An attractive aspect of linear classifiers, relative to some others
discussed below is that they tend to be quite interpretable. This
is because the entries of the vector $\boldsymbol{w}$ can be called
\emph{classification loadings\index{loadings!classification}}, in
the spirit of PCA loadings as discussed around Figures \ref{fig:2dToy-PC1}
and \ref{fig:2dToy-PC2}. Thus the drivers of the classification can
be understood using loadings plots\index{loadings!plots|(} of the
type shown in Figure \ref{figB-X:PanCanDWDLoadings}.

\emph{\index{Fuzzy classification}Fuzzy classification} is a sometimes
useful parallel to classification, where instead of trying to predict
the class, the goal is to instead estimate an analog of the posterior
probabilities $P\left\{ \widetilde{y}=y|\boldsymbol{x}\right\} $.
This gives not only the class prediction, but also a measure of confidence
in that prediction. The terminology ``fuzzy'' in this context seems
to go back to \citet{Zadeh1965FuzzySets}.

A comprehensive reference to many important ideas and methods in classification
can be found in \citet{DudaHartStork2001PattRecog}. Another good
overview can be found in \citet{Maclachlan2004DiscAnal}.

Section \ref{secJ-A:ClassicalMethods} studies some classical statistical
approaches to classification, but some non-standard ideas are revealed.
Sections \ref{secJ-B:KernelMethods} and \ref{secJ-C:SVM} provide
a graphical introduction to the kernel trick\index{kernel trick}
and Support Vector Machines, which were at the heart of the early
days of machine learning. An improvement of the Support Vector Machine
called Distance Weighted Discrimination is introduced and discussed
in Section \ref{secJ-D:DWD}. 

Classification is an area where the ``every dog has its day'' quote
in the preface is particularly applicable, as a large number of methods
(each with relative strengths and weaknesses) are available.

\section{Classical Methods\label{secJ-A:ClassicalMethods} }

Perhaps the simplest classification method is called \emph{one nearest
neighbor}\index{classifier!nearest neighbor}, where $c\left(\boldsymbol{x}_{0}\right)$
is the $\widetilde{y}_{i}$ of the closest training point to $\boldsymbol{x}_{0}$.
An extension that generally has better Bayes risk is the $k$ nearest
neighbors rule, where a vote is taken among the $k$ nearest neighbors
to $\boldsymbol{x}_{0}$. Besides the appealing property of simplicity
a convenient aspect of these methods is that they are \emph{distance\index{distance|(}
based} as discussed in Chapter 5. Thus they are easily constructed
in any \index{metric}metric space. A serious drawback of the $k$
nearest neighbor rule is that it requires choice of $k$, which is
generally as hard as bandwidth\index{bandwidth!selection} selection
as discussed in Section \ref{secN-C:SmooParSel}.

As discussed in Section \ref{secD2:Preprocess-Trans}, another very
simple classifier is the \emph{\index{Mean Difference|(}Mean Difference}
(MD)\index{classifier!Mean Difference} or \emph{Centroid\index{classifier!centroid}}
approach. One way of thinking of the MD is to simply take the class
whose mean is closest. But MD can also be thought of as a linear classifier
of the form (\ref{eqJ.A:LinearClasser}). In particular the separating
hyperplane has the line between the class means as the normal direction,
and intercept as the midpoint between. A notational benefit of the
$\pm1$ class labels is that $\widehat{\boldsymbol{w}}=\frac{\sum_{i=1}^{n}\widetilde{y}_{i}\widetilde{\boldsymbol{x}}_{i}}{\left\Vert \sum_{i=1}^{n}\widetilde{y}_{i}\widetilde{\boldsymbol{x}}_{i}\right\Vert }$
and $\widehat{\beta}=n^{-1}\sum_{i=1}^{n}\widetilde{\boldsymbol{x}}_{i}$,
where here and in the following the hat symbol is used to indicate
an estimated quantity. When both classes are balanced and Gaussian
with identity covariance, the MD is the likelihood ratio rule so it
is Bayes risk optimal. 

In situations where the variables may have arbitrarily different scalings
(e.g. are measured in different units) it is sensible to first rescale\index{scale}
each variable in the full training data set by subtracting its mean
and dividing by its standard deviation (and applying the same standardization\index{standardization}
to $\widetilde{\boldsymbol{x}}_{0}$). This variant of MD appears
to have been named \emph{\index{classifier!Naive Bayes}Naive Bayes}
by \citet{Domingos1997NaiveBayes}. This is Bayes risk optimal in
balanced Gaussian settings where the common covariance matrix is diagonal.
The same issues (about the strengths and weaknesses of accentuating
variables by standardization) raised during the analysis of the Two
Scale Curves data\index{Two Scale Curves data} in Figures \ref{fig:ToyPCAnScaling}
and \ref{fig: ToyPCAnScalingCorr} apply in this context as well.

The MD and Naive Bayes methods work well when there is no correlation
among the variables, but otherwise they can often be substantially
improved as demonstrated in using the \index{Shifted Correlated Gaussians data|(}Shifted
Correlated Gaussians example in Figure \ref{figJ.A:WhyNotMD}. The
left panel displays a $d=2$ dimensional training data set with $n_{+}=20$
class +1 data points colored red and $n_{-}=20$ labeled $-1$ colored
blue, for a total of $n=n_{+}+n_{-}=40$ points. Green circled x signs
denote the two sample means and the circled plus is the overall mean.
The green direction vector from the overall mean shows the MD direction.
The MD separating hyperplane (normal to the MD direction) is shown
as the dashed green line. Projections of the training data onto the
subspace (line) generated by the MD direction are shown in the right
panel. Substantial overlap of the projected data indicates this is
not an outstanding direction for linear classification. The left panel
suggests that a direction that is more orthogonal to the point clouds
could give much better performance.

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigAWhyNotMD}

\caption{Left panel shows the Shifted Correlated Gaussians training data set
with color indicating classes, green direction vector (based at the
circled plus sign) shows the MD direction, and separating hyperplane
is the green dashed line. Right panel is projections in the MD direction.
Shows poor class separation performance of MD when training variables
are correlated.\label{figJ.A:WhyNotMD}}
\end{figure}
The above ideas can be used to give a graphical derivation of the
widely used and studied \emph{\index{Linear Discriminant Analysis|(}Linear
Discriminant Analysis} (LDA). In the case of binary classification\index{binary classification|)},
another common name for LDA is \emph{\index{Fisher Linear Discriminant}Fisher
Linear Discriminant}. Figure \ref{figJ.A:WhyNotMD} suggests that
MD failed because it only uses the class means and ignores the clearly
important covariance information. LDA was proposed by \citet{Fisher1936LDA}
who made the key observation that when both classes are assumed to
be Gaussian with \emph{common} covariance matrices the likelihood
ratio classifier is linear, i.e. of the form (\ref{eqJ.A:LinearClasser}).
For this reason almost all (2 class) derivations of LDA are done from
a Gaussian likelihood viewpoint. As in other situations, this can
lead to the mistaken impression that LDA somehow requires Gaussianity
to be effective. To counter this, a fully nonparametric\index{nonparametric}
visual introduction to LDA is given in Figure \ref{figJ.B:FLDIntro}.

The same Shifted Correlated Gaussians data appears in the upper left
panel. For simplicity of exposition, only the balanced case where
the number of points in each class are the same is treated here. It
is straightforward but notationally heavier to treat the general case
by using appropriate class sample size weighting everywhere. Critical
components of LDA are the class means $\overline{\boldsymbol{x}}^{+}$
and $\overline{\boldsymbol{x}}^{-}$, shown as circled x signs in
both the left panel of Figure \ref{figJ.A:WhyNotMD}, and the lower
left of Figure \ref{figJ.B:FLDIntro}. Also important are the two
within class sample covariance matrices, $\widehat{\boldsymbol{\varSigma}}^{+}$
and $\widehat{\boldsymbol{\varSigma}}^{-}$ (recall the definition
(\ref{eqB.F:SampleCovariance})). Under the assumption that the underlying
population covariance matrices are the same, it make sense to estimate
that common matrix by \emph{pooling} the within class sample covariance
matrices to get $\widehat{\boldsymbol{\varSigma}}^{w}$, by taking
a weighted (according to class size) average of $\widehat{\boldsymbol{\varSigma}}^{+}$
and $\widehat{\boldsymbol{\varSigma}}^{-}$. Some algebra shows that
$\widehat{\boldsymbol{\varSigma}}^{w}$ can be written in a form that
is quite similar to the overall sample covariance matrix $\widehat{\boldsymbol{\varSigma}}$,
with the important exception that the overall mean is replaced by
the respective class means. For this reason $\widehat{\boldsymbol{\varSigma}}^{w}$
is called the sample \emph{\index{within class covariance}within
class covariance}. Note that $\widehat{\boldsymbol{\varSigma}}^{w}$
effectively captures the covariance structure (with strong positive
correlation) that causes problems for MD in Figure \ref{figJ.A:WhyNotMD}.
Hence it is natural to undo that by a \emph{\index{sphering|(}sphering}
operation, i.e. transforming each data point to 
\begin{equation}
\widetilde{\boldsymbol{z}}_{i}=\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2}\widetilde{\boldsymbol{x}}_{i}.\label{eqJ.B:SphereTrans}
\end{equation}
The success of this is seen in the plot of the $\widetilde{\boldsymbol{z}}_{i}$
in the upper right panel of Figure \ref{figJ.B:FLDIntro}, where each
class indeed looks spherically distributed. That is the ideal situation
for the MD classifier, which is applied in the lower right panel.
The corresponding transformed class means $\overline{\boldsymbol{z}}^{+}$
and $\overline{\boldsymbol{z}}^{-}$ are shown as magenta circled
x signs, with the overall mean (just the average of the two class
means) shown as the circled plus. The MD direction in this transformed
space is shown as the magenta direction vector based at the overall
mean, and the orthogonal separating plane appears as the dashed magenta
line. Classification of the data in this space could be performed
by applying the \index{transformation!sphering|(}transformation to
new data, but it is more insightful to transform the linear classification
algorithm back to the original space. Note that the circled plus center
point is transformed back to $\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2}\left(\frac{1}{2}\overline{\boldsymbol{z}}^{+}+\frac{1}{2}\overline{\boldsymbol{z}}^{-}\right)$,
while the back transformed normal vector is proportional to $\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2}\left(\overline{\boldsymbol{z}}^{+}-\overline{\boldsymbol{z}}^{-}\right)$.
Now inverting the sphering transformation (\ref{eqJ.B:SphereTrans})
gives the LDA centerpoint $\widehat{\boldsymbol{\mu}}_{LDA}=\frac{1}{2}\overline{\boldsymbol{x}}^{+}+\frac{1}{2}\overline{\boldsymbol{x}}^{-}$
and the direction vector 
\begin{equation}
\widehat{\boldsymbol{w}}_{LDA}=\frac{\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1}\left(\overline{\boldsymbol{x}}^{+}-\overline{\boldsymbol{x}}^{-}\right)}{\left\Vert \left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1}\left(\overline{\boldsymbol{x}}^{+}-\overline{\boldsymbol{x}}^{-}\right)\right\Vert }.\label{eqJ.C:FLDdirn}
\end{equation}
Thus this is a linear classifier (\ref{eqJ.A:LinearClasser}) where
the intercept is the inner product $\widehat{\beta}_{LDA}=\widehat{\boldsymbol{w}}_{LDA}^{t}\widehat{\boldsymbol{\mu}}_{LDA}$.
The corresponding separating hyperplane is shown as the magenta dashed
line in the lower left panel, which is clearly sensible for this data
set. 

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigBFLDIntro}

\caption{Graphical introduction to LDA. Top left panel shows the same Shifted
Correlated Gaussians data. Top right panel shows the result of transforming
using the within class sphering operation. Bottom right panel shows
the result of applying the MD classifier to the transformed data.
Bottom left panel shows the effectiveness of the MD classifier, which
is essentially an inversion of the sphering transformation\index{transformation!sphering|)}.
\label{figJ.B:FLDIntro}}
\end{figure}
This derivation of LDA is also interpretable in terms of the important
idea of Mahalanobis distance\emph{\index{distance!Mahalanobis}},
proposed by \citet{Mahalanobis1936}. Given independent random vectors
$\widetilde{\boldsymbol{x}}_{1},\;\widetilde{\boldsymbol{x}}_{2}\sim N_{d}(\boldsymbol{\mu},\boldsymbol{\Sigma})$,
a natural distance measure is 
\[
d_{M}=\left[\left(\widetilde{\boldsymbol{x}}_{1}-\widetilde{\boldsymbol{x}}_{2}\right)^{t}\boldsymbol{\varSigma}^{-1}\left(\widetilde{\boldsymbol{x}}_{1}-\widetilde{\boldsymbol{x}}_{2}\right)\right]^{1/2}.
\]
This distance is interpretable along the familiar lines of ``standard
deviations from the mean''. Note that Mahalanobis distance is computable
as the standard Euclidean distance applied to the sphered variables
$\boldsymbol{\varSigma}^{-1/2}\widetilde{\boldsymbol{x}}_{1}$ and
$\boldsymbol{\varSigma}^{-1/2}\widetilde{\boldsymbol{x}}_{2}$. Thus
it is essentially Euclidean distance in the standardized space shown
in the right column of Figure \ref{figJ.B:FLDIntro}.

LDA has received far more attention than MD in the statistical literature,
perhaps because it was perceived as doing much better in situations
such as those studied in Figure \ref{figJ.B:FLDIntro}, and being
roughly similar otherwise. However that ignores the increasingly important
high dimensional case where $d>n$. This is a challenge because then
the sample covariance matrix $\hat{\varSigma}^{w}$ is not invertible.
Various domains of high dimensional data are discussed in Chapter
\ref{chapM:HDasy}, where an important principal is that when matrix
inverses are not available, often generalized inverses contain the
needed information. The Moore-Penrose \index{generalized inverse}generalized
inverse (also called pseudo-inverse\index{pseudo-inverse}) is a well
defined version of generalized inverse that essentially consists of
inverting nonzero eigenvalues (while leaving zero eigenvalues as is).
See Section \ref{subsecP-A-B:SVD} for more discussion of eigenvalues.
While use of the Moore-Penrose generalized inverse gives a well defined
version of LDA, performance in high dimensions is rather poor, as
seen in Figure \ref{figJ.C:HighDClass}.

The panels in each column of Figure \ref{figJ.C:HighDClass} are based
on the $n_{+}=n_{-}=20$ data vectors of dimension $d=10,\;40,\;200$
dimensions in rows 1, 2, 3 respectively. The $d=200$ vectors in each
column were simulated independently from Gaussian distributions having
mean vectors $\boldsymbol{\mu}_{+}=\left(+2.2,0,\cdots,0\right)^{t}$
and $\boldsymbol{\mu}_{-}=\left(-2.2,0,\cdots,0\right)^{t}$ for the
$+1$ (red) and the $-1$ (blue) classes respectively. To keep the
realizations as similar as possible, the $d=10$ and $40$ versions
are the subvectors of the first few entries. Thus each data set has
the same signal, with noise level increasing with the dimension from
left to right. Each panel shows the projection of the training dataset
onto a linear classification direction, i.e. a choice of $\widehat{\boldsymbol{w}}$
in (\ref{eqJ.A:LinearClasser}). In the top row $\widehat{\boldsymbol{w}}$
is the MD direction. Note that for $d=10$ in the top left panel,
the projections seem close to what is expected from projection onto
the optimal direction $\boldsymbol{w}_{OPT}=\left(1,0,\cdots,0\right)^{t}$.
This is verified by displaying the angle between this estimated $\widehat{\boldsymbol{w}}_{MD}$
and $\boldsymbol{w}_{OPT}$ which is a very small $16^{\circ}$ in
this relatively easy case. Performance of MD remains good for $d=40$,
as seen in the center top panel, although the angle to $\boldsymbol{w}_{OPT}$
has grown somewhat. Even for the quite challenging case of $d=200$
(top right panel) the MD direction continues to separate the classes,
and perhaps surprisingly the classes seem to have actually moved farther
apart, while the angle to $\boldsymbol{w}_{OPT}$ has also gotten
larger. This phenomenon will be explained using high dimensional concepts
developed in Section \ref{secM-B:HDLSS}. In the middle row, $\widehat{\boldsymbol{w}}$
is the LDA direction $\widehat{\boldsymbol{w}}_{LDA}$. Performance
in the middle left panel ($d=10$) is reasonably good, although as
expected (because of the overhead of estimating $\widehat{\boldsymbol{\varSigma}}^{w}$)
not as good as MD, neither visually nor in terms of the angle to $\boldsymbol{w}_{OPT}$.
At first glance the $d=40$ performance of LDA might be considered
excellent, in the sense of separating the classes very well with respect
to the projected within class variances. But the angle to $\boldsymbol{w}_{OPT}$
of $80^{\circ}$ gives substantial pause as that suggests that this
classifier will have very poor generalizability properties when it
is given new data. At $d=200$, the LDA falls apart completely with
a very poor angle to $\boldsymbol{w}_{OPT}$ as well as big overlap
of the projected classes. 

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigCHighDClass}

\caption{Comparison of classification methods MD (top row), LDA (middle) and
MDP (bottom), over varying dimensions. Angles to the optimal in degrees
for each direction are included. Shows MD very robust, LDA has poor
high dimensional performance, and perhaps surprisingly good behavior
of MDP. \label{figJ.C:HighDClass}}
\end{figure}
The bottom row shows a method (that is not commonly considered in
discussions of classification methods) called the Maximal Data Piling\index{Maximal Data Piling|(}
(MDP) direction. MDP was first defined and studied in detail in \citet{AhnMarron2010MDP}.
In Figure \ref{figJ.C:HighDClass}, MDP has quite different properties
from the above two methods in the higher dimensional cases. The $d=40$
direction is quite spurious with an angle to $\boldsymbol{w}_{OPT}$
of $89^{\circ}$, yet the projections have a very systematic structure
of all the red points piled up at the \emph{same value}, and similarly
for the blue points. This piling of the two classes happens with probability
one (for distributions absolutely continuous with respect to Lebesgue
measure) for $d\geq n$, and is where the name MDP came from. In particular,
as noted in \citet{AhnMarron2010MDP}, out of all directions where
both classes completely pile, MDP has the largest separation of the
class piling points. For higher $d=200$ the piling points spread
substantially, and once again perhaps surprisingly the angle to $\boldsymbol{w}_{OPT}$
becomes quite comparable with that of MD. This suggests that in high
dimensions the generalizability of MDP is fairly comparable to that
of the optimal likelihood ratio rule, MD. Another high dimensional
quirk of MDP, shown by \citet{AhnMarron2010MDP}, and explained in
the PhD dissertation \citet{Miao2015Dissertation}, is that for auto-correlated
errors, MDP actually performs \emph{better} than MD\index{Mean Difference|)}.

It is also seen in \citet{AhnMarron2010MDP} that the form of MDP
is rather surprisingly similar to LDA. In particular the only difference
is that the within class covariance matrix $\widehat{\boldsymbol{\varSigma}}^{w}$
is replaced by the overall covariance matrix $\widehat{\boldsymbol{\varSigma}}$,
so that 
\begin{equation}
\widehat{\boldsymbol{w}}_{MDP}=\frac{\left(\widehat{\boldsymbol{\varSigma}}\right)^{-1}\left(\overline{\boldsymbol{x}}^{+}-\overline{\boldsymbol{x}}^{-}\right)}{\left\Vert \left(\widehat{\boldsymbol{\varSigma}}\right)^{-1}\left(\overline{\boldsymbol{x}}^{+}-\overline{\boldsymbol{x}}^{-}\right)\right\Vert }.\label{eqJ.D:MDPdirn}
\end{equation}
 In addition, it is seen that for $d\leq(n-2)$ the direction vectors
$\widehat{\boldsymbol{w}}_{LDA}$ and $\widehat{\boldsymbol{w}}_{MDP}$
are the \emph{same}, despite the different numerators in (\ref{eqJ.C:FLDdirn})
and (\ref{eqJ.D:MDPdirn}). This can be seen by a careful look at
the $d=10$ left center and bottom panels of Figure \ref{figJ.C:HighDClass}.
Note that both the horizontal and vertical coordinates of each symbol
are exactly the same. This is an example of when using order in the
data as heights in jitter plot\index{jitter plot}s conveys more information
than using a random height, as discussed in Section \ref{secB2:Data-Visualization}.
Daniel Peņa has explained this phenomenon as illustrated in Figure
\ref{figJ.D:MDPIntro} using the above viewpoint. This uses the same
Shifted Correlated Gaussians data set as Figures \ref{figJ.A:WhyNotMD}
and \ref{figJ.B:FLDIntro} and this time only the corresponding bottom
row is shown. The sphering that occurs in the right panel is done
by $\left(\widehat{\boldsymbol{\varSigma}}\right)^{-1/2}$ (as opposed
to the $\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2}$
used to get LDA above). Instead of individually sphering each class
separately (as done in the right column of Figure \ref{figJ.B:FLDIntro}),
now the full data set is sphered. The two point clouds are still apparent,
but they have been adjusted so that the overall sample covariance
matrix is the identity. Note that this transformation also nicely
handles the offset of the class means that was the original problem
with MD highlighted in Figure \ref{figJ.A:WhyNotMD}. In particular,
the sample means following this transformation are shown as magenta
circled x signs, and the corresponding separating hyperplane is very
sensible. The result of back transformation, following the same lines
as in Figure \ref{figJ.B:FLDIntro} result in the linear classification
rule in the lower left panel, which again is the \emph{same} as in
the LDA case.

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigDMDPIntro}

\caption{Visual introduction to MDP, using the same data and format as Figure
\ref{figJ.B:FLDIntro}. Shows how (for $d<n$) overall sphering\index{sphering|)}
results in the same classification rule as LDA. \label{figJ.D:MDPIntro}}
\end{figure}
A natural conclusion of the above is that the ubiquitous LDA should
be \emph{replaced} by MDP as the simple classifier of choice. This
is because they are the \emph{same} in the case $n<d$, and when they
are different MDP is superior. Furthermore MDP is (very slightly)
simpler to compute.

An unfortunately misleading presentation of these issues can be found
in Section 3.8.2 of \citet{DudaHartStork2001PattRecog}, where it
is stated that LDA optimizes the ratio of the difference between the
means and the sum of the standard deviations of the projections onto
the normal vector. While the standard LDA formula is given, it is
in fact MDP that has the stated properties. That presentation is not
mathematically incorrect, because it is assumed in Section 3.7.2 there
that $d<n$, in which case all quantities are the same. But it is
important to be aware that several statements in Section 3.8.2 are
incorrect when $d\geq n$. 

In situations where the within class covariance matrices are not similar,
one can still appeal to the \index{classifier!Gaussian Likelihood Ratio}Gaussian
Likelihood Ratio to construct a classifier. While it is straightforward
to compute the classification rule $c_{GLR}$ the classification boundary
(\index{separating surface|)}separating surface) is determined by
quadratic equations and explicit representation is rather complicated
and case-wise. This method is not treated further here, as it has
been the subject of many classical texts, and because its high dimensional
performance is even worse than that of LDA as illustrated in Figure
\ref{figJ.C:HighDClass}.

Another branch of classification methods reviewed in Chapter 4 of
\citet{DudaHartStork2001PattRecog}, assumes that each class has a
smooth probability density, and uses smoothing\index{smoothing} methods
as discussed in Chapter \ref{chapN:Smoothing} to construct classifiers.
Early results on optimal rates of convergence (of the type discussed
in Section \ref{secN-B:SmoothBasics}) were established in \citet{Marron1983}.
For much more on optimality in related contexts see \citet{DevroyeGyorfi1996PatRecog}.

A point that will become relevant in Section \ref{secJ-C:SVM} is
that all of the methods developed in this section can be thought of
as being based on probability distributions.

While the above discussion has been about binary classifiers for the
sake of simplicity, an extension to $K$ class methods of LDA is well
worth describing. In the case of $K>2$ the likelihood approach is
no longer useful, but a direction based approach is available. From
Figure \ref{figJ.B:FLDIntro} and (\ref{eqJ.C:FLDdirn}) it is clear
that the LDA direction vector points in the direction of the difference
of class means after appropriately adjusting by the pooled within
class covariance structure. For multiple means, the single direction
between means can be usefully replaced by a PCA of the set of mean
vectors. That would be computed\index{computed} as an eigen analysis
of the covariance matrix of the set of class means, which is sometimes
called the \emph{between class covariance} matrix $\widehat{\boldsymbol{\varSigma}}^{b}$.
As LDA got some benefit from adjusting the MD using the within class
covariance matrix $\widehat{\boldsymbol{\varSigma}}^{w}$, the same
approach makes sense here. This results in a set of direction vectors
(analogs of the PCA loadings vectors detailed in Figures \ref{fig:2dToy-PC1}
\& \ref{fig:2dToy-PC2} of Section \ref{secBA.1:Foundations_Terminology}),
called \emph{Canonical Variate Analysis }(CVA)\index{Canonical Variate Analysis|(}
in \citet{MardiaKentBibby1979MultivarAnal} and\emph{ \index{Multiple Discriminant Analysis}Multiple
Discriminant Analysis} in \citet{DudaHartStork2001PattRecog}, as
previously mentioned in Section \ref{sec F.D:AlternativeDirns}. Potential
differences between the Gaussian likelihood and algebraic derivation
of these methods have been resolved by \citet{vanMeegen2019LDA}

The canonical variates can be calculated as an eigen analysis of $\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1}\widehat{\boldsymbol{\varSigma}}^{b}$,
by solving $$\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1}\widehat{\boldsymbol{\varSigma}}^{b} \boldsymbol{w} = \lambda \boldsymbol{w}$$
In practice we use the eigen-decomposition of the matrix $\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2}\widehat{\boldsymbol{\varSigma}}^{b}\left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2}$
which is symmetric (hence more amenable to standard eigenanalysis
algorithms) and has the same eigenvalues but with corresponding eigenvectors
$\boldsymbol{v}$ related to the canonical variate vectors $\boldsymbol{w}$
by $\boldsymbol{w} = \left(\widehat{\boldsymbol{\varSigma}}^{w}\right)^{-1/2} \boldsymbol{v}$.
This computation is a type of generalized eigen analysis\index{generalized eigen analysis}\index{eigen analysis|)},
where eigenvectors of one matrix are computed with respect to another
matrix (see \citet{DrydenMardia2016}, page 298). Both the high and
low dimensional issues of the relationship between LDA and MDP based
on replacing $\widehat{\boldsymbol{\varSigma}}^{w}$ with $\widehat{\boldsymbol{\varSigma}}$
are likely to apply here, although this does not seem to have been
explored as yet. 

An example of using Canonical Variate Analysis is seen in Figure \ref{figJ.O:CVA}.
In this dataset there are six groups of great ape skull landmarks:
males and females for each of Gorilla (gor), Chimpanzee (pan) and
Orangutan (pongo). The dataset is available in the {\tt shapes} package
in R, \citet{Dryden2019RShapes}. Generalized Procrustes analysis
(as discussed in Section \ref{subsec:Generalized-Procrustes-analysis})
is carried out to remove location, rotation and scale\index{scale}
information to retain the shapes, which lie in the non-Euclidean shape
space. Tangent space coordinates are obtained at the overall mean
and both PCA and CVA are carried out. We plot the first two PC scores
in the left hand figure (which summarize 37.5\% and 28.1\% of the
shape variability), and the first two canonical variate scores\index{scores!canonical}
are shown on the right. The CVA gives much better visual separation
of the groups than PCA. For CVA it is clear that there is good separation
in shape for all six groups, except there is some overlap between
the male and female chimpanzees. Careful inspection of the axis labels
reveals how PCA directions maximize projected variation which is quite
different from maximizing class differences as done by CVA. Note that
Canonical Variate Analysis\index{Canonical Variate Analysis|)} is
equivalent to Canonical Correlation Analysis\index{Canonical Correlation Analysis}
(Section \ref{subsecP-B-B:CCA}) between the data matrix $\boldsymbol{X}$
and a new dummy indicator matrix $\boldsymbol{Y}$ of zeros and ones
to indicate which group the observation belongs to (see \citet{MardiaKentBibby1979MultivarAnal},
Exercise 11.5.4, page 330). 

\begin{figure}
\includegraphics[scale=0.3]{ChapterJ/ApePCA}\includegraphics[scale=0.3]{ChapterJ/CVAplot}

\caption{PCA and CVA\index{Canonical Variate Analysis|)} of a dataset of six
groups of great apes: males (m) and females (f) for each of Gorilla
(gor), Chimpanzee (pan) and Orangutan (pongo). The left figure shows
the first two PC scores and the right figure shows the canonical variate
scores. There is much better separation of the groups using the canonical
variates compared to PCA. \label{figJ.O:CVA}}
\end{figure}

\section{Kernel Methods\label{secJ-B:KernelMethods}\index{}}

Kernel methods, as developed in the machine learning\index{machine learning}
community are usefully viewed as a type of transformation, as discussed
in Section \ref{secD4:Standardization}, with the goal of making data
more accessible to linear analysis. The fundamental idea goes back
to \citet{AizermanBravermanRozoner1964}. A canonical classification
example called the Donut 2-d data\index{Donut 2-d data|(} is shown
in Figure \ref{figJ.E:DonutPolyEmbed}. Each panel is based on the
same data set, with $n_{+}=100$ points in the positive class (donut
hole) shown as red plus signs from 0.6 times a standard Gaussian distribution
(i.e. $N\left(0,(0.6)^{2}\right)$, and $n_{-}=100$ negative class
points (forming the donut) appearing as blue circles generated using
polar coordinates that are 2.5 + 0.5 times a standard Gaussian (i.e.
$N\left(5,(0.5)^{2}\right)$) in the radial component and uniformly
over the angular component. Training data points $\left(x_{1,i},x_{2,i}\right)^{t}$
for $i=1,\cdots,200$ are plotted on the $x_{1}$ (horizontal) and
$x_{2}$ (vertical) axes in each panel. Performance of the classifier
based on each polynomial embedding is shown by treating each background
pixel\index{pixel|(} (recall this terminology from Section \ref{secF.A:Heatmap-Views})
as a new point and classifying it with yellow denoting a $+$ (red)
classification and cyan indicating $-$ (blue). Both colors are scaled
towards white for pixels near the boundary. In the upper left panel,
LDA classification is done in the original $\mathbb{R}^{2}$, where
LDA has very poor performance because it is constrained to be a linear
classifier of the form (\ref{eqJ.A:LinearClasser}), none of which
can separate this data well. 

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigEDonutPolyEmbed}

\caption{Donut 2-d toy example, classified using LDA in various polynomial
embedding spaces. Background reflects classification of each pixel
location. Shows increasing benefits of higher order embedding. \label{figJ.E:DonutPolyEmbed}}
\end{figure}
The remaining panels of Figure \ref{figJ.E:DonutPolyEmbed} show the
results of \emph{polynomial kernel embedding}\index{kernel embedding|(}\index{kernel embedding!polynomial},
where data are nonlinearly transformed, i.e. embedded, in a higher
dimensional space where the linear method LDA is then applied. Again
impact of the embedding process on the classification is shown by
classifying the background pixels in the same way. The upper right
panel shows the results of embedding into $\mathbb{R}^{3}$ by mapping
the $i$-th data point to $\left(x_{1,i},x_{2,i},x_{1,i}^{2}\right)^{t}$.
That mapping puts the data onto a sheet that is warped in the $x_{1}$
direction into a parabolic cylinder. LDA slices that surface in $\mathbb{R}^{3}$
in a horizontal way so that pixels with $x_{1}$ coordinate near 0
are assigned to the $+$ (red) class and the rest are $-$ (blue),
which gives a marked improvement over the $\mathbb{R}^{2}$ version
of LDA in the upper left panel. The lower left panel is similar to
the upper right, except now the embedding is to $\left(x_{1,i},x_{2,i},x_{2,i}^{2}\right)^{t}$,
so the parabolic cylinder in $\mathbb{R}^{3}$ is now folded in the
$x_{2}$ direction, resulting in a rotation of the classification
results. The lower right panel studies the embedding to $\mathbb{R}^{4}$,
using $\left(x_{1,i},x_{2,i},x_{1,i}^{2},x_{2,i}^{2}\right)^{t}$.
It is hard to conceive of the 2-d sheet containing the data in $\mathbb{R}^{4}$,
so interpretation is not so straightforward. However, the yellow-cyan
coloring makes it clear that classification is happening mainly in
terms of $x_{1,i}^{2}+x_{2,i}^{2}$ (i.e. squared distance to the
origin), giving a very good result in this case. In particular, this
nonlinear transformation has taken a data set that is clearly intractable
to linear methods and renders it to be quite amenable to the simple
LDA.

While polynomial embedding is clearly ideal for the Donut 2-d toy
data in Figure \ref{figJ.E:DonutPolyEmbed} it also has some drawbacks
that include potential lack of flexibility as well as poor interpretability
and non-obvious choice of polynomial degree. This is illustrated in
Figure \ref{figJ.F:Checkerboard}, which again starts with a data
set in $\mathbb{R}^{2}$. But this time the data follow a checkerboard
pattern that is even less amenable to linear classification. The format
of Figure \ref{figJ.F:Checkerboard} is the same as for Figure \ref{figJ.E:DonutPolyEmbed}.
Each class consists of $8$ standard Gaussian clusters of $25$ points
each, that are spaced 6 units apart in both directions. The polynomial
embeddings shown in Figure \ref{figJ.E:DonutPolyEmbed} were not sufficiently
flexible to give reasonable classification. Somewhat more flexible
is the cubic embedding to $\left(x_{1,i},x_{2,i},x_{1,i}^{2},x_{2,i}^{2},x_{1,i}^{3},x_{2,i}^{3}\right)^{t}$
shown in the left panel of Figure \ref{figJ.F:Checkerboard}, which
still gives very poor classification performance.

Much better performance appears in the right panel of Figure \ref{figJ.F:Checkerboard}.
This uses a much different embedding, using kernel density estimation\index{density estimation}
ideas as discussed in Section \ref{secN-A:NoHisto}. In particular
the data are embedded in $\mathbb{R}^{49}$ by mapping $\left(x_{1,i},x_{2,i}\right)^{t}$
to 
\[
\left(\phi\left(\left(x_{1,i}-g_{1,1}\right)/h\right)\cdot\phi\left(\left(x_{2,i}-g_{2,1}\right)/h\right),\cdots,\phi\left(\left(x_{1,i}-g_{1,49}\right)/h\right)\cdot\phi\left(\left(x_{2,i}-g_{2,49}\right)/h\right)\right)^{t},
\]
where $\phi$ is the standard Gaussian density, where $\left(g_{1,1},g_{2,1}\right),\cdots,\left(g_{1,49},g_{2,49}\right)$
are a $7\times7$ grid of kernel centerpoints, and where the spread
of each Gaussian kernel is controlled by the bandwidth\index{bandwidth}
$h=1$. This much more complicated embedding results in excellent
classification where each training point is correctly classified.
This shows why Gaussian kernels, a commonly used type of what are
called \emph{radial basis functions\index{radial basis function|(}}
in the kernel machine literature, have been very popular.

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigFCheckerboard}

\caption{Checkerboard 2-d toy example. Shows poor classification performance
by a polynomial embedding, but a much better performance using LDA
on a radial basis function embedding. \label{figJ.F:Checkerboard}}
\end{figure}
A clever variation of the kernel embedding idea is \emph{kernel PCA}\index{kernel PCA|(},
proposed by \citet{Scholkopf1997KernelPCA}. This is illustrated in
Figure \ref{figJ.K:KernelPCA}, based on the Four Cluster \index{Four Cluster data}
toy data set which is the basis of several examples in Chapter \ref{chapK:Clustering}.
This data set has two spherical clusters shown as blue circles and
cyan squares, a stretched cluster of green plus signs, and a small
two point cluster (perhaps outliers?) shown as red x signs. Kernel
PCA starts with conventional PCA performed in the kernel space, with
visualization\index{visualization|(} done in the original object
space\index{object space}. The latter comes from projecting every
point in the object space (essentially each pixel\index{pixel|)}
value) onto each direction vector, and coloring the image with the
score, using gray for 0, and intensity of black (white) for magnitude
of positive (negative respectively) score. 

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigKkernelPCA}

\caption{Four Cluster toy example illustrating kernel PCA. Shows first three
modes of variation. First (left) contrasts green cluster with union
of blue and cyan. Second (center) contrasts blue and cyan clusters.
Third (right) splits the stretched green cluster. \label{figJ.K:KernelPCA}}
\end{figure}
The first component (i. e. \index{modes of variation!kernel PCA}mode
of variation, shown in the left panel) highlights a coarse scale clustering
which provides a contrast between the stretched green cluster, while
treating the union of the blue and cyan as a single coarse scale cluster.
The second PC mode of variation refines the clustering view by separating
the blue and cyan clusters. The third mode seems to split the long
green cluster into two subclusters. The higher components (not shown
here) are similarly interpretable as a further split, into 3, of the
stretched green cluster, followed by a focus on the two point red
cluster. Gaussian kernel PCA and also other kernel methods depend
on the choice of the window width. In Figure \ref{figJ.K:KernelPCA}
this was chosen by trial and error to give good interpretation of
the modes of variation. See \citet{Ahn2010RBFbandwidth} for further
useful ideas on this choice.

While examples like this are clearly compelling, this view of kernel
PCA does not seem to have had much impact on real data analysis. This
appears to be due to the fact that this visualization requires a low
dimensional (really $d=2$) object space. A much more successful visualization
based on kernel PCA is the \index{t-SNE|(}t-SNE (t-distribution Stochastic
Neighbor Embedding) approach of \citet{MaatenHinton2008tSNE}. This
is usefully viewed as a type of inversion of the Gaussian kernel PCA,
using a Cauchy kernel (i.e. very heavy tails) analog of kernel PCA\index{kernel PCA|)}.
Inversion here is done much as in MDS (recall Multi-Dimensional Scaling\index{Multi-Dimensional Scaling|(}
from Section \ref{secE.2:Multi-dimensional-scaling}), where the key
idea was to find a set of representer\index{representer|(}s (of each
data object) whose position best fits the set of distances between
data objects. This approach is called \emph{auto-encoding\index{auto-encoding}}
in the machine learning\index{machine learning} literature. In t-SNE
the Cauchy kernel embedding\index{kernel embedding!Cauchy}s of the
representers are fit to the positions of the data objects in the kernel
space. The heavy tails of the Cauchy kernel, relative to the Gaussian
give a representation which keeps nearby points close to each other,
while somehow pushing away points that are less close, which can visually
accentuate clusters. This can result in impressive looking graphics,
although it does rely on choice of a tuning parameter, and does not
resolve issues such as which clusters represent true underlying structure,
as done e.g. by the SigClust\index{SigClust} method of Section \ref{secL2:SigClust}.
Finally note that the name t-SNE comes from the fact that the standard
Cauchy distribution is also Student's t distribution with one degree
of freedom.

Figure \ref{figJ.L:tSNE} shows three applications of t-SNE to the
same Four Cluster toy data set. The 3 panels show different values
of the \emph{perplexity}\index{t-SNE!perplexity}, which is a tuning
parameter. The default perplexity of 30 is shown in the center panel
and gives a reasonable view of the data, except for the placement
of the two red points in the small isolated cluster. The right panel,
with perplexity 60, may be somewhat better, although the blue and
cyan clusters seem to be uncomfortably close. The perplexity 12 visualization
in the left panel has the red points even farther apart and two of
the blue points separated from the rest. An important aspect is that
the global structure is quite arbitrary, but the local clusters tend
to be reasonably highlighted, which is the goal of t-SNE.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterJ/OODAbookChpJFigLtSNE}

\caption{Example of t-SNE visualization for the same Four Cluster data set,
for different values of the perplexity parameter. Default of 30 is
in the center panel, but the most reasonable choice may be 60 as shown
in the right panel. \label{figJ.L:tSNE}}
\end{figure}
An application of t-SNE to the \index{Osteo-Arthritis data}Osteo-Arthritis
data set of \citet{NelsonFang2019FNIHtSNE} is shown in Figure \ref{figJ.N:tSNEonOA},
which is Figure 4 of that paper. This data set has four distinct subpopulations,
which are first dichotomized as patients whose disease progressed
over time (blue and cyan) versus those who did not (red and green),
and second split on males (blue and red) versus females (cyan and
green). The t-SNE visualization nicely highlights these as four clusters,
suggesting these are useful partition\index{partition}s of this data
set. However, the visual \index{clusters!relationships between}relationships
between clusters appear to be symmetric. The black lines and numbers
show the results of a quantitative investigation of the difference
between these clusters, using the \index{DiProPerm}DiProPerm hypothesis
test (discussed in Section \ref{secL1:DiProPerm}) of the difference
between each pair of subpopulations. The numbers are Z-scores which
reflect statistical significance on the standard Gaussian scale, so
larger than 2 is statistically significant. Larger Z-scores indicate
stronger statistical significance. In Figure \ref{figJ.N:tSNEonOA}
this shows that the male versus female subpopulations are actually
far more distinct than progression versus not. This highlights the
fact that because of its stress on locality t-SNE cannot be counted
on to visually indicate strength of relationships between clusters. 

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterJ/FNIHfig4crop}

\caption{Application of t-SNE to Osteo-Arthritis data. Highlights four subsets
in the data, Red for male non-progressors, Green for female non-progressors,
Cyan for female progressors and Blue for male progressors. DiProPerm
Z-scores show that, contrary to t-SNE visual impression, some subgroups
are actually much more distinct than others. \label{figJ.N:tSNEonOA} }
\end{figure}
At the time of this writing, t-SNE visualization may be fading in
popularity in the bioinformatics community, in favor of UMAP proposed
by \citet{McinnesHealy2018umap}. Like t-SNE\index{t-SNE|)}, UMAP
seeks a set of visual representer\index{representer|)}s of high dimensional
data, but deep ideas from differential geometry and fuzzy topology
provide the motivation for this method.

A fundamental issue is how a kernel embedding relates to the inner
products that lie at the heart of many data analytic methods such
as PCA and LDA. In the examples shown in Figures \ref{figJ.E:DonutPolyEmbed}
and \ref{figJ.F:Checkerboard} an \emph{explicit \index{kernel embedding!explicit}kernel
embedding} is used where the embedding to the kernel space is performed
first, and then LDA\index{Linear Discriminant Analysis} is explicitly
applied in that space. But there are major computational advantages
to first computing the matrix of inner products of the data objects,
and then mapping those to the kernel space, which is usefully labeled
\emph{implicit kernel embedding}\index{kernel embedding!implicit}.
The details of the relationship between explicit and implicit embeddings
have been carefully laid out in the kernel machines literature. A
fundamental tool in that area is Mercer's theorem that relates embedding
types. See \citet{CristianiniShaweTaylor2000BookSVM,ShaweTaylorCristianini2004BookKernels}
for a good introduction, and \citet{ScholkopfSmola2002BookKernelSVM}
for a more detailed investigation. 

One implicit kernel embedding is the use of the matrix of centered
inner products, which are used in the computation of distance matrices
in MDS (Multi-Dimensional Scaling\index{Multi-Dimensional Scaling|)},
Section \ref{secE.2:Multi-dimensional-scaling}). In particular in
classical MDS the ``centered inner product matrix'' is formed to
compute the matrix of pairwise Euclidean distances, see 14.2.1 of
\citet{MardiaKentBibby1979MultivarAnal}). Although sometimes a horseshoe
effect\index{horseshoe effect} may be present in such representations
as mentioned at the end of Section \ref{secE.2:Multi-dimensional-scaling},
the higher dimensional \index{kernel embedding!higher dimensional}kernel
embedding may be useful in classification tasks.

The great flexibility of kernel methods comes with great danger for
overfitting\index{overfitting|(}. The \index{BiDirectional Discrimination}BiDirectional
Discrimination ideas of \citet{Huang2012BDD} overcome this with a
methodology that involves a product of two (or more) linear direction
vectors for a more controllable enhanced flexibility. 

While the above examples make clear the strong data analytical advantages
available from kernel embedding, an interesting limitation has been
pointed out by \citet{ElKaroui2010}. The fundamental result presented
there is that, in the case of very high dimension, kernel methods
give the same classification performance as conventional linear methods.
In particular, the analysis of \citet{ElKaroui2010} uses random matrix
theory\index{random matrix theory} asymptotics, where both the sample
size and the dimension $n,d\rightarrow\infty$, as studied in Section
\ref{secM-A:RandomMatrix}.

\section{Support Vector Machines\label{secJ-C:SVM}}

There is an important contrast to notice between the classical methodologies
reviewed in Section \ref{secJ-A:ClassicalMethods} and the kernel
embedding\index{kernel embedding|)} context of Section \ref{secJ-B:KernelMethods}.
The methods in the former are strongly rooted in probability distributions
(the foundations of the popular likelihood approach to statistics),
which are very well developed and best understood in flat Euclidean
spaces $\mathbb{R}^{d}$. However the key to the improvements demonstrated
in the latter is to map the data onto a curved manifold (as studied
in Chapter \ref{chapG:Manifold-Data-Analysis}). For example the parabolic
cylinders in the upper right and lower left panels of Figure \ref{figJ.E:DonutPolyEmbed},
as well as the harder to understand 2-d curved surface in $\mathbb{R}^{4}$
that underlies the classification result shown in the lower right
panel. Because probability distributions on manifolds are far less
well developed (e.g. recall the discussion about \index{Central Limit Theorem}Central
Limit Theorems on manifolds in Section \ref{subsecG-A-A:CLTon Manifolds})
and intuitively understood, it makes sense to consider non-probabilistic
approaches to the invention of classification methods.

Such considerations appear to have motivated the invention of the
\emph{Support Vector Machine\index{Support Vector Machine|(}} (SVM)
by \citet{Vapnik1982SVM,Vapnik1995SVM}. A very readable short introduction
to the main ideas can be found in \citet{Burges1998SVMtutorial}.
While the SVM is aimed at data lying on curved manifolds, its characteristics
are best understood in a simple Euclidean space such as the 2-d toy
example shown in Figure \ref{figJ.G:SVMIntro}. The data set shown
in both panels is similar in spirit to the Shifted Correlated Gaussians\index{Shifted Correlated Gaussians data|)}
data  in Figures \ref{figJ.A:WhyNotMD}, \ref{figJ.B:FLDIntro} and
\ref{figJ.D:MDPIntro} except there are now only $n_{+}=15$ red circles
in the $+$ class, and $n_{-}=15$ blue circles in the $-$ class.
These data are \emph{separable}\index{separable classification data},
in the sense that there is a hyperplane between the classes with all
of one class on one side, and the other class entirely on the other
side. In fact there are many such separating planes, some of which
are shown using colored dashed lines in the left panel of Figure \ref{figJ.G:SVMIntro}.
A useful view of the linear SVM is that it seeks to find the ``best''
among all such separating hyperplanes. Note this way of thinking is
quite different from the mostly probability distribution based approaches
in Section \ref{secJ-A:ClassicalMethods}. The SVM notion of best
is illustrated in the right panel of Figure \ref{figJ.G:SVMIntro},
where given a candidate separating hyperplane the residuals $\tilde{r}_{i}$
from the projection of each data point $\boldsymbol{x}_{i}$ onto
the plane are shown as thin magenta line segments. SVM chooses the
hyperplane which is farthest from every data point in the sense of
maximizing the smallest of these projected distances, as shown in
the right panel. The minimizing distance is called the \emph{margin}\index{classifier!margin},
and planes shown as black dotted lines parallel to the separating
plane are shown. Note that there are 3 points on these lines, which
are the data points that achieve the smallest \index{distance|)}distances.
These are called \emph{support vectors\index{classifier!support vectors}}
(hence the name of the method) and are highlighted with black boxes. 

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterJ/OODAbookChpJFigGSVMPIntro}

\caption{Graphical introduction to the SVM. Left panel shows separable toy
data, together with several separating hyperplanes, Right panel shows
SVM direction vector (solid thick magenta line) and orthogonal SVM
separating plane (dashed magenta line). Thin magenta lines show projected
residuals. SVM margins are dotted lines with support vectors highlighted
by black boxes. \label{figJ.G:SVMIntro}}
\end{figure}
To formulate the linear SVM as an optimization problem, let the $d\times n$
matrix $\boldsymbol{X}=\left[\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{n}\right]$
and similarly combine the class labels into an $n\times1$ vector
$\boldsymbol{y}=\left[y_{1},\cdots,y_{n}\right]^{t}$. Note that the
tilde notation used above to denote random quantities is deliberately
not used here to reflect the common machine learning\index{machine learning}
viewpoint that the input data are just numbers and not random variables.
Also let $\boldsymbol{Y}$ denote the $n\times n$ diagonal matrix
with $y_{1},\cdots,y_{n}$ on its diagonal. Given a direction vector
$\boldsymbol{w}\in\mathbb{R}^{d}$ (thus $\left\Vert \boldsymbol{w}\right\Vert =1$)
and an intercept $\beta\in\mathbb{R}$, the $n$ vector of magenta
residual lengths in the right panel of Figure \ref{figJ.G:SVMIntro}
can be written in the form $\check{\boldsymbol{r}}=\boldsymbol{Y}\boldsymbol{X}^{t}\boldsymbol{w}+\beta\boldsymbol{y}$.
The simplicity of this representation is due to the choice of $\pm1$
as class labels. In the separable case, the SVM optimization problem
is then written as $\underset{\boldsymbol{w},\beta}{\max}\underset{i}{\min}\:\check{r}_{i}$.
This can be solved using quadratic programming, by introducing a new
variable $\tau$, and maximizing $\tau$ subject to $\check{r}_{i}\geq\tau$
for $i=1,\cdots,n$. Since $\check{\boldsymbol{r}}$ scales with $\boldsymbol{w}$
and $\beta$, the maximization of $\tau$ subject to $\left\Vert \boldsymbol{w}\right\Vert =1$
is equivalent to minimizing $\left\Vert \boldsymbol{w}\right\Vert $
subject to $\tau=1$. When the data cannot be assumed to be separable,
nonnegativity of the residuals is maintained by the use of slack variables
$\boldsymbol{\xi}=\left[\xi_{1},\cdots,\xi_{n}\right]^{t}$ which
are incorporated into the modified matrix of residuals $\boldsymbol{r}=\boldsymbol{Y}\boldsymbol{X}^{t}\boldsymbol{w}+\beta\boldsymbol{y}+\boldsymbol{\xi}$.
Given a tuning parameter $\lambda$, these modified residuals are
then constrained to be nonnegative in the more general version of
the SVM optimization: 
\begin{equation}
\underset{\boldsymbol{w},\beta,\boldsymbol{\xi}}{\min}\left(\left\Vert \boldsymbol{w}\right\Vert ^{2}-\lambda\boldsymbol{1}_{1,n}\boldsymbol{\xi}\right)\label{eqJ.E:SVMopt}
\end{equation}
subject to the constraints $r_{i}\geq0$, $\xi_{i}\geq0$, for $i=1,\cdots,n$,
using the notation from (\ref{eq:ones matrix}) of $\boldsymbol{1}_{d,n}$
for a $d\times n$ matrix of ones. Small values of $\lambda$ give
a solution called the \emph{\index{hard margin SVM}hard margin} SVM,
which attempts to put as many data points as possible outside the
margins (essentially the dotted black lines shown in the right panel
of Figure \ref{figJ.G:SVMIntro}). Larger values of $\lambda$ relaxes
that goal and allows more \emph{violators}, i.e. points not outside
of the margins.

Much more on SVMs from a kernel learning viewpoint can be found in
the books by \citet{CristianiniShaweTaylor2000BookSVM,ShaweTaylorCristianini2004BookKernels}
and \citet{ScholkopfSmola2002BookKernelSVM}. A more statistically
oriented overview can be found in \citet{LinWahba2002SVM}. For approaches
to choice of the SVM tuning parameter, see \citet{Joachims2000SVM}
and \citet{Wahba1999SVM,Wahba2001SVM}.

A detailed analysis of several important properties of the SVM can
be found in \citet{CarmichaelMarron2017SVManal}. One is a study of
the solution path, through the space of linear classifiers (\ref{eqJ.A:LinearClasser}),
as a function of the tuning parameter $\lambda$ for a given data
set. At one end of this path is the hard margin classifier studied
in Figure \ref{figJ.G:SVMIntro}. At the other end of the path is
a trimmed version of the MD\index{Mean Difference}, which is exactly
the MD when the classes are balanced. The path moves between these
ends only over a rather narrow range and computationally useful bounds
are given that guarantee the tuning parameter is in either ending
state.

The SVM approach has been extended to image analysis data objects
lying in a curved manifold\index{manifold} (in the sense treated
in Chapter \ref{chapG:Manifold-Data-Analysis}) by \citet{Sen2008SVmanifold}. 

Another important issue is extension of the SVM to the multi-class
case $K>2$. A simple extension involves repeated use of the binary
SVM. There are two common approaches. First one trains a classifier
on each class versus all the rest and chooses the class with the best
result (in terms of projection onto the SVM direction vector). The
second one runs all of the pairwise classifications, and then chooses
a class by voting. As noted in \citet{Friedman1996PolyClass}, the
second option is generally better. \citet{LeeLinWahba2004MultiSVM}
proposed integrating the multi-class task into the SVM optimization
problem.

While the SVM has become a workhorse method for many machine learning
tasks, especially its hard margin form (as in Figure \ref{figJ.G:SVMIntro})
has a deficiency in high dimensions that results in some loss of classification
performance, as well as poor visualization properties. This is illustrated
in Figure \ref{figJ.H:DWDvSVM}, whose format and data set are very
similar to Figure \ref{figJ.C:HighDClass}. Again the same realization
is used in all panels and the data are essentially Gaussian with shifted
means, but this time the dimension is $d=50$. Each panel is a projection
of the data onto a direction in $\mathbb{R}^{50}$. 

The top left panel is the optimal direction, in the sense that this
is the direction in which the underlying population means have been
shifted (recall the population means have been shifted to $\pm2.2$.
Each panel displays an angle at the bottom, which is angle in degrees
to this direction. Projections onto the MDP direction are shown in
the top right panel. As expected from the discussion in Section \ref{secJ-A:ClassicalMethods},
the data in each class pile up completely at a single point.  To make
this happen, MDP must clearly exploit small scale noise artifacts
particular to this realization of the data, which should damage its
generalizability as a classification rule. This is reflected by the
large angle of $64^{\circ}$ to the optimal direction. 

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterJ/OODAbookChpJFigHDWDvSVM}

\caption{Illustration of data piling for MDP (top right panel), SVM (lower
left), DWD\index{Distance Weighted Discrimination|(} (lower right)
on 50 dimensional Gaussian data set with means shifted along the optimal
direction shown in the upper left. Angles to the optimal direction
in degrees are shown for each direction. Reveals that SVM shares some
of the data piling of MDP, which is avoided by DWD. \label{figJ.H:DWDvSVM} }
\end{figure}
The SVM direction shown in the lower left panel gives a better separation
between the projected data, because it optimizes the spacing between
classes (thus larger than MDP which must also make the data completely
pile up). This will clearly result in much better classification performance
because the angle to the optimal direction has now gone down to $37^{\circ}$.
However the SVM also reveals some data piling in the sense of both
classes having a fair number of points that project to a common value
at the SVM \emph{margin,} shown as the thin black dotted lines in
the right panel of Figure \ref{figJ.G:SVMIntro}. These are the support
vectors (basically ties for best in the SVM optimization) which tend
to be quite numerous in higher dimensions. As discussed for the MDP\index{Masximal Data Piling}
this also indicates a lesser extent to which the SVM feels small scale
noise artifacts which should also cause some loss in classification
performance. 

A related linear classification method, that has been motivated by
this problem with SVM is \emph{Distance Weighted Discrimination} (DWD)
developed in Section \ref{secJ-D:DWD}. Note that projection onto
the DWD direction (bottom right of Figure \ref{figJ.H:DWDvSVM}) has
a smaller angle with the optimal of only $29^{\circ}$ suggesting
generally better classification performance (which has been verified
using simulations in \citet{MarronToddAhn2007}). Another useful property
of DWD is the approximately Gaussian distribution of the projected
classes, in sharp contrast to the apparently opposing triangle distributions
for the SVM. As discussed in Section \ref{secJ-D:DWD} this has serious
implications for both data visualization and batch adjustment, as
discussed in \citet{Benito2004DWDadjustment} and \citet{LiuPerouMarron2009CrossPlatVisualization}.

\section{Distance Weighted Discrimination\label{secJ-D:DWD}}

The lower panels of Figure \ref{figJ.H:DWDvSVM} demonstrate the high
dimensional data piling problem that is endemic to SVM, and that DWD
can overcome it. The impact of this on data visualization is explored
in Figure \ref{figJ.I:BadSVMvis} using the Pan Cancer gene expression
data from Figure \ref{figB-R:PanCanDWD}. The left panel of Figure
\ref{figJ.I:BadSVMvis} is a zoomed in view of the bottom middle panel
of that figure, showing projections on the DWD directions trained
on BLCA vs. COAD and on BRCA vs. OV. The right panel shows projections
on the corresponding hard margin SVM directions. The big picture view
as to how the clusters relate to each other is the same. But quite
distracting to the human perception of these clusters in Figure \ref{figJ.I:BadSVMvis}
are the clear vertical lines of data points in the yellow (COAD) and
magenta (BLCA) clusters, as well as the nearly horizontal lineups
of points in the red (BRCA) and cyan (OV) clusters, which are all
consequences of SVM partial data piling. Note also that in each of
those clusters, the remaining data points all lie on the same side
of the line as expected from the ideas illustrated in the lower left
panel of Figure \ref{figJ.H:DWDvSVM}. This is why DWD is generally
recommended for visualization tasks where the focus is on potential
differences between pairs of given subgroups of data.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterJ/OODAbookChpJFigIBadSVMvisA}

\caption{Left panel is the scatterplot view (projections on DWD directions)
of the Pan Cancer data\index{Pan Cancer data} from the bottom middle
panel of Figure \ref{figB-R:PanCanDWD}. Corresponding SVM version
appears in the right panel. Shows how data piling strongly impacts
the use of SVM direction vectors in high dimensional data visualization\index{visualization|)}
by creating spurious lines of data points. \label{figJ.I:BadSVMvis}}

\end{figure}
As noted in \citet{MarronToddAhn2007}, DWD achieves this improvement
through a modification of the SVM optimization problem in (\ref{eqJ.E:SVMopt}).
In the separable case, the idea is to change the way that the residuals
$\check{r}_{i}$ impact the result, by letting all of them have some
influence\index{influence}, instead of just the smallest ones. However,
the smallest ones should play a more important role, which is achieved
by studying the inverse of the residual, $\frac{1}{\check{r}_{i}}$.
As any residual approaches 0, the inverse becomes large, so their
sum works like a set of poles pushing the separating hyperplane towards
the middle of the region between the data sets. Thus, the separable
version of DWD uses the optimization problem $\underset{\boldsymbol{w},\beta}{\min}\sum_{i=1}^{n}\frac{1}{\check{r}_{i}}.$
This is extended to the general case using the same slack variable
approach used for the SVM in Section \ref{secJ-C:SVM}. Employing
notation there, the DWD optimization problem is: given the tuning
parameter $\lambda$
\begin{equation}
\underset{\boldsymbol{w},\beta,\boldsymbol{\xi}}{\min}\left(\sum_{i=1}^{n}\frac{1}{\check{r}_{i}}+\lambda\boldsymbol{J}_{1,n}\boldsymbol{\xi}\right)\label{eqJ.F:DWDopt}
\end{equation}
subject to the constraints $\check{r}_{i}\geq0$, $\xi_{i}\geq0$,
for $i=1,\cdots,n$. Note that this is the linear version of DWD.
As for SVM, there are also appropriate kernelized versions. 

The numerical solution to (\ref{eqJ.F:DWDopt}) proposed by \citet{MarronToddAhn2007}
is based on the \emph{\index{second order cone programming}second
order cone} algorithm SDPT3, available at \citet{Toh2009SDPT3}. Like
the nominally simpler \emph{\index{quadratic programming}quadratic
programming} that underpins SVM, second order cone programming is
a fast greedy search algorithm. The computational speed of this approach
scales very well with respect to the dimension $d$, which was important
for earlier bioinformatics data sets. However, scaling with respect
to sample size $n$ is much worse. This has motivated the development
of \index{FastDWD}FastDWD by \citet{LamMarronToh2018FastDWD}, which
uses a different optimization approach based on a semiproximal \emph{alternating
direction method of multipliers}\index{alternating direction method of multipliers}.
The computational speed of FastDWD scales very well with respect to
both $d$ and $n$, and is sometimes even more efficient than the
highly optimized LIBLINEAR and LIBSVM implementations of SVM. 

Another issue is selection of the tuning parameter $\lambda$ in (\ref{eqJ.F:DWDopt}).
As noted in Section \ref{secJ-C:SVM}, the issue is quite important
for SVM, because depending on the particular data set at hand, best
performance can come from $\lambda$ anywhere between the hard margin
SVM and a trimmed version of the MD. The situation is much different
for DWD, mostly because the analog of the hard margin in that situation
gives much more broadly useful results. Thus \citet{MarronToddAhn2007}
recommend achieving this with the choice $\lambda=\frac{100}{d_{t}^{2}}$
(the number 100 is chosen as a ``large number''), where $d_{t}$
is a useful notion of \index{scale}scale of the data computed as
$d_{t}=median\left\{ \left\Vert \boldsymbol{x}_{i}-\boldsymbol{x}_{i'}\right\Vert :y_{i}=+1,y_{i'}=-1,\right\} $.
An interesting open problem is analysis of DWD of the type done in
\citet{CarmichaelMarron2017SVManal}.

As first shown in \citet{Benito2004DWDadjustment} DWD provides a
useful direction for the removal of batch effects from genetic data.
At first this good performance was just an empirical observation,
but the underlying reason was later revealed in \citet{LiuPerouMarron2009CrossPlatVisualization}.
The main idea is illustrated in Figure \ref{figJ.J:BatchAdjDWD}.
The left panel shows a toy data set that models a typical situation
in the area of batch adjustment. Symbols represent data objects with
$n_{+}=200$ plus signs data collected by one lab, while $n_{-}=200$
circles come from another. Colors highlight two subtypes whose difference
is the focus of the experiment. Note that for each lab separately
(both the plusses and the circles), there is a clear systematic difference
between the colors, but there is a $4:1$ imbalance between the subtypes.
To boost statistical power it is natural to pool the data over labs,
perhaps by subtracting the mean of the data from each lab. While this
works well when the subtypes (colors) are balanced within each lab
(symbols) it is challenging because of the imbalance here as shown
by the sample means (green x signs). Subtracting the mean from each
lab corresponds to bringing the labs together by shifting towards
each other along the MD\index{Mean Difference} direction (dashed
green line). The result of this manipulation is shown in the bottom
right panel. Note that this MD adjustment is profoundly unsatisfactory
because the difference between the colors (the goal of this study)
has been seriously reduced to the extent that the main body of the
subtypes overlap.

Better performance is available from use of DWD for batch adjustment.
The DWD direction is shown in the left panel of Figure \ref{figJ.J:BatchAdjDWD}
as the dashed pink line, which is much closer to the ideal vertical
direction than the green dashed MD direction. The result of sliding
the labs along the DWD direction is shown in the upper right panel
of Figure \ref{figJ.J:BatchAdjDWD}. This shows much better performance
of DWD for pooling the labs (symbols) as now there is a substantial
gap between the resulting colors. The DWD direction works better for
batch adjustment because it is driven by the DWD separating hyperplane
shown as the solid pink line. The $\frac{1}{\check{r}_{i}}$ poles
at each data point pushes this plane away, which results in a much
better direction for batch adjustment. The reason that DWD adjustment
was seen to work well for cancer genetic data in \citet{Benito2004DWDadjustment}
seems to be due to the strong heterogeneity of cancer, which results
in unknown unbalanced subtypes as illustrated in Figure \ref{figJ.J:BatchAdjDWD}.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterJ/OODAbookChpJFigJBatchAdjDWD}

\caption{Toy data set (left panel) showing value of DWD batch adjustment relative
to an MD approach. Colors represent important subtypes to be studied,
while symbols reflect irrelevant batch effects. MD shown in green,
and poor results of simply subtracting the MD from the batches appear
in the bottom right panel. DWD direction shown as a dashed pink line,
with much better result of sliding the batches along the DWD direction
in the upper right panel. \label{figJ.J:BatchAdjDWD}}
\end{figure}
Batch effects are becoming an increasingly important issue as there
is a growing awareness that understanding very diverse diseases such
as cancer need data sets larger than any single lab can collect. This
has led to collaborative efforts such as The Cancer Genome Atlas,
\citet{Weinstein2013TCGA}, to integrate data across a large number
of labs. While great care is taken about protocols and many other
aspects of data collection, it is still impossible to completely eliminate
unintentional batch differences. Handling these issues is the domain
of the statistical area of \emph{data heterogeneity}\index{data heterogeneity},
discussed in \citet{Marron2017BigData}. Current approaches include
\emph{explicit} approaches such as DWD and other approaches to batch
adjustment such as the empirical Bayes approach of \citet{Johnson2007BatchAdjust},
the Surrogate Variable Analysis method of \citet{Leek2012SVA}, as
well as \emph{implicit} approaches to statistical analysis that are
\index{robust methods}robust against data heterogeneity as proposed
by \citet{MeinshausenBuhlmann2015InhomoData,Buhlmann2015Inhomo}.

A sparse version of DWD can be found in \citet{WangZou2016SparseDWD}.
\citet{WangZou2018DWD} coined the term \emph{generalized DWD\index{generalized DWD}}
for a variation of DWD with important ramifications for batch adjustment
and robustness against heterogeneity, where the impact of the poles
$\frac{1}{\tilde{r}_{i}}$ created by the residuals is explicitly
controlled using a power $q$, i.e. $\frac{1}{\tilde{r}_{i}^{q}}$.
The case $q=1$ is conventional DWD as defined in (\ref{eqJ.F:DWDopt}),
and larger values of $q$ magnify the effect illustrated in Figure
\ref{figJ.J:BatchAdjDWD} in the sense of pushing the solid pink line
more towards horizontal. This results in a separating direction (dashed
pink line) closer to the optimal vertical direction, i.e. improved
robustness against heterogeneity. This also provides a connection
between DWD and SVM which seems to be essentially the limit as $q\rightarrow\infty$.

As discussed in Section \ref{secJ-C:SVM}, extension of DWD to the
multi-class case can be done by use of pairwise DWD in a one versus
all scheme, or in the preferred multiple pairwise vote scheme. An
intrinsically multi-class version of DWD (in the sense that the optimization
problem is formulated in a multi-class way) has been proposed by \citet{Huang2013MultiDWD}. 

As noted in Section \ref{secC2:Compositional-Data-Analysis}, \citet{XiongDittmerMarron2015}
addressed a challenge in virology using data objects lying on the
unit simplex. That challenge was a classification problem aimed at
using DNA to identify the presence of viruses in blood samples, by
aligning DNA in a given sample to the virus genome. The data objects
were vectors of counts of DNA-seq reads, that were divided by the
total counts to handle amplification effects. Such data objects are
compositional data\index{data objects!compositional} (i.e. lie on
the unit simplex, as discussed in Section \ref{secC2:Compositional-Data-Analysis})
as they sum to 1. The special challenge of that data was that while
positive cases tended to lie fairly close to the center of the unit
simplex, the negative cases tended to go off in many very different
directions. The setting is essentially a very high dimensional version
of the donut\index{Donut 2-d data} example of Figure \ref{figJ.E:DonutPolyEmbed}.
However it is not clear that any of the kernel methods of Section
\ref{secJ-B:KernelMethods} can handle the problem because test data
tends to go off in directions that are not present in the training
data. This concern is well justified in the left panel of Figure \ref{figJ.N:RadialDWD}
which is from Figure 6.2 of the PhD dissertation \citet{Xiong2015radialDWDdissertation},
reproduced as Figure 4 of \citet{XiongDittmerMarron2015}, showing
projections on the SVM separating direction in the radial basis function
kernel space of the data. That classifier is trained on some positive
cases shown as red plus signs and negative cases shown as blue circles.
As in many other graphics, the vertical coordinates of these points
are just jitter to give visual separation. The remaining symbols are
independent test cases. The pink asterisks are positive human samples
which all tend to be close to the center of the simplex and thus are
correctly classified. The green asterisks are approximately positive
samples representing viruses related to the positive class but from
non-human mammals. Note that these are also correctly classified (i.e.
to the right of the classification boundary shown as the vertical
dashed line). More problematic are the cases shown as gray x signs,
which are a large number of negative cases. While many of these are
correctly classified (to the left of the boundary) a large number
are not. These are cases that are far from the center of the simplex,
but in directions that are not present in the training data, and hence
are not reflected in the kernel space classification boundary. \citet{XiongDittmerMarron2015}
address this problem by replacing the separating hyperplane used at
many points above with a \emph{separating hypersphere}\index{separating hypersphere}.
This is an option to keep in mind other situations where one class
is fairly centrally located and the other is spread very widely around,
in a large number of other directions. \citet{XiongDittmerMarron2015}
formulate a DWD-like approach to find the center and radius of the
separating hypersphere (other linear classifier\index{linear classifier|)}s
could be adapted to spherical versions in a parallel way). This results
in \emph{Radial DWD\index{Radial DWD}} whose performance on the virus
data is shown in the right panel of Figure \ref{figJ.N:RadialDWD}.
The horizontal axis is essentially the radius of a polar coordinate
system, that has been shifted so that 0 is at the separating hypersphere.
Cases shown at the left end (labeled ``-inf'') are those with all
0 counts, which are not on the simplex (and thus are naturally classified
as negative). Colors and symbols are the same as in the left panels,
with red plus signs for positive training data and blue circles for
negative training cases. Once again, the positive test cases (pink
asterisks) and related positives (green asterisks) are all correctly
classified. But this time all of the negative test cases are correctly
classified, verifying the above idea about using a separating sphere
to correctly capture negative test cases that depart from the center
in directions not present in the training data. An even more dramatic
simulated example can be found in Figures 5 and 6 of \citet{XiongDittmerMarron2015}. 

\begin{figure}[H]
\includegraphics{ChapterJ/OODAbookChpJFigNRadialDWD}

\caption{Illustration of good performance of Radial DWD in virus detection
in the right panel, contrasted with poor performance of SVM using
a classical radial basis function\index{radial basis function|)}
kernel. Demonstrates ability of Radial DWD\index{Distance Weighted Discrimination|)}
to correctly classify test cases that are very different from anything
present in the training data. \label{figJ.N:RadialDWD}}
\end{figure}

\section{Other Classification Approaches }

There are other popular classification methods not treated in detail
here. \emph{Neural networks} were very popular in the 1980s and 1990s.
A great success of that era was voice recognition software which went
into widespread corporate use in the early 2000s. That success led
to gross over-advertisement of the methodology, which resulted in
many failed promises that in turn gave the whole field a bad name.
However neural networks have more recently come back in a very strong
way, under the new name of \emph{Deep Learning}\index{Deep Learning}.
While many of the basic ideas are the same, much of the current success
(which is becoming almost completely dominant in fields such as computer
vision) seems to be due to both size of modern data sets and also
the far greater computational resources that are currently available.
See \citet{Bengio2013RepresentationLearning} and \citet{Goodfellow2016Deep}
for good overview of this area. A current major challenge to deep
learning methods is \emph{interpretability}, i.e. the intuitive understanding
of the drivers of differences between classes in classification contexts.
In contrast, for simple linear methods, loadings\index{loadings!plots}
plots as shown using the Pan Cancer data\index{Pan Cancer data} in
Figure \ref{figB-X:PanCanDWDLoadings}, provide simple and direct
insights into the main drivers of class differences. 

Another approach to good interpretability is the tree based approaches,
such as the \index{CART}CART algorithm proposed in \citet{BFOS1984CART}
The basic tree algorithms tended to get stuck in local optima, motivating
many proposed solutions. The most widely accepted solution of this
problem is the \emph{random forest\index{random forest}} approach
of \citet{Breiman2001RandomForests}. This is still widely used and
is a frequent candidate in comparative classification studies, along
with kernel SVM\index{Support Vector Machine|)}.

A quite different approach to combining a set of tree classifiers
into a single more powerful classifier is one case of the idea of
``combining weak learners'' called \emph{\index{boosting}boosting}
by \citet{FreundShapire1995Boosting,FreundShapire1999ShortBoost}.
See \citet{Friedman2000Boosting} for an interesting statistical perspective
on boosting.

These many methods, all with their relative strengths and diverse
properties make classification an area where the ``every dog has
its day'' principle quite clear.

As noted above, all classification methods have potential for \emph{overfitting}\index{overfitting|)},
i.e. resulting in a classification rule\index{classification rule|)}
that is very effective on the training data but is not very generalizable,
i.e. useful for new data. This potential typically gets worse for
more complicated methods. This has led to a strong desire to measure
performance of a classification method. Basic measurement is typically
done using misclassification rates assessed on an independent \emph{test
data} set. This can be done either theoretically, or using a simulation
study. But misclassification rates assume a particular relative weighting
of classes. In situations where the desired weighting is not clear,
an approach that essentially considers all weightings is the \index{Area Under the Curve}Area
Under the Curve (AUC) summary of the Receiver Operating Characteristic
curve\index{Receiver Operating Characteristic curve} introduced in
Section \ref{secD2:Preprocess-Trans}.

For the analysis of a single real data set an objective approach is
\emph{\index{cross-validation}cross-validation}, see e.g. \citet{Stone1974CV},
\citet{Kohavi1995CV} and \citet{Arlot2010surveyCV}. The idea is
to divide the data into separate training and testing subsets, on
which the method is respectively trained and independently tested.
Random choice is generally a good idea for this, and Monte Carlo variation
should be damped by repeating the process and averaging. One approach
generates repeats through independent random train-test divisions
done at each step. Another common approach is $k$ fold cross-validation
where the data are partition\index{partition}ed into $k$ subsets
of approximate size $n/k.$ Repeats consist of each subset being used
in turn as the test set, with the others used for training. Choices
for $k$ can range from $2$ to $n$. The choice $k=2$ is generally
not recommended as the size of the training data set is only $\frac{n}{2}$
which can be a substantially different classification\index{classification|)}
problem than that being studied. The choice $k=n$ is essentially
\emph{leave one out} cross-validation. Common choices of $k$ are
4 or 5. Generally recommended is to explicitly choose random subsets
so that the classes are proportionally represented in the training
and testing sets, as opposed to using randomization to make this only
approximately true. Cross-validation is commonly used for two distinct
tasks: for evaluation of the performance of a classifier and for choice
of tuning parameters. When both are done together, it is important
to use none of the test set for the latter, but instead to do a separate
cross-validation completely within the training set\index{training data|)}.
While cross-validation is broadly useful, it is important to keep
in mind that it too can suffer strongly from sample variation, as
illustrated in the context of smoothing\index{smoothing} in Section
\ref{secN-C:SmooParSel}.
