
\chapter{Distance Based Methods\label{chapE:Distance-Based-Methods}}

This chapter is about OODA methods that are based only on distance\index{distance|(}s
between data objects. An advantage of such approaches is that they
are quite broadly useful, which can be important in exotic data spaces
such as manifold\index{manifold|(}s, or tree\index{tree}/graph\index{graph}
spaces, where simply developing data analytic methods can be challenging.
In particular, the only structure needed on the object space\index{object space}
for such methods is the presence of a \emph{metric\index{metric|(}}
(i.e. distance). In some situations, even less may be needed, e.g.
perhaps only a \emph{pseudo metric\index{metric!pseudo}} or even
a \emph{dissimilarity}\index{dissimilarity}, but these will not be
further pursued here. 

In this chapter, the symbol $\delta$ will be used to denote metrics,
i.e. distance functions. As in any elementary analysis textbook, as
apparently originally formulated by \citet{Frechet1906DefnMetrics},
a metric is a two argument function from some space to the real numbers
which is symmetric, nonnegative with the value 0 taken if and only
if the two arguments are the same, that satisfies the triangle inequality.
Given a set of data objects in an arbitrary space
\[
\left\{ \chi_{i}:i=1,\cdots,n\right\} ,
\]
and a distance $\delta$, the corresponding $n\times n$ symmetric
\emph{distance matrix\index{distance matrix|(}} is
\begin{equation}
D=\left[\begin{array}{cccc}
0 & \delta(\chi_{1},\chi_{2}) & \cdots & \delta(\chi_{1},\chi_{n})\\
\delta(\chi_{2},\chi_{1}) & 0 &  & \vdots\\
\vdots &  & \ddots & \delta(\chi_{n-1},\chi_{n})\\
\delta(\chi_{n},\chi_{1}) & \cdots & \delta(\chi_{n},\chi_{n-1}) & 0
\end{array}\right].\label{eq:ChpEeqC:DistanceMatrix}
\end{equation}
Working with data objects in this type of format is rather common
in the machine learning\index{machine learning} literature, see \citet{CristianiniShaweTaylor2000BookSVM},
\citet{ScholkopfSmola2002BookKernelSVM} and \citet{ShaweTaylorCristianini2004BookKernels}
for good overviews. Indeed a central idea in that area is called the
\emph{kernel trick\index{kernel trick}}, where one works with a matrix
of inner products (closely related to the usual distance in Euclidean
spaces), with the important goal of large computational benefits.

As noted in Sections \ref{secBA.1:Foundations_Terminology} and \ref{secB4:Further-Statistical-Tasks},
a quick and dirty default approach to OODA is to first find a metric,
and then simply carry out a distance based analysis. Such notions
of \emph{center} are discussed in Section \ref{secE.1:Fr=0000E9chet-mean}.
Methods for understanding variation about the center, in the spirit
of PCA, using only distances are explored in Section \ref{secE.2:Multi-dimensional-scaling}.
Clustering is another set of data analytic methods, that are frequently
based only on distances between data objects. Cluster analysis\index{cluster analysis}
is discussed in Chapter \ref{chapK:Clustering}. 

While the metric first approach gives useful insights of some types,
a drawback is that it tends to only yield \index{scores|(}scores
(displaying relationships between\index{data objects!relationships between|(}
data objects e.g. for the \index{Twin Arches data}Twin Arches data
in Figure \ref{fig:50d-ToyPCAScatMat}) and \emph{not} full modes
of variation\index{modes of variation!MDS} as defined in Section
\ref{subsecBA.1.1:FormalizingModes}. In particular, it does not easily
give the one dimensional indexed sets in the data space that convey
all of the insights available from full modes. A more careful use
of particular data space structure can result in full modes of variation
in many OODA contexts, as discussed in Chapters \ref{chapG:Manifold-Data-Analysis},
\ref{chapH:FDA-Curve-Registration} and \ref{chapI:Tree-Structured-Data}.

A critical aspect of metric based analyses is that choice of distance
has a major impact on the results. For a particularly dramatic example
consider the \emph{discrete metric\index{metric!discrete}}
\[
\delta_{D}\left(\chi,\varsigma\right)=\begin{cases}
\begin{array}[t]{c}
0\\
1
\end{array} & \begin{array}[t]{c}
\mathrm{for}\;\chi=\varsigma\\
\mathrm{for\;}\chi\neq\varsigma
\end{array}\end{cases}.
\]
This exists for any space, but is useless for OODA analyses, because
it contains no information about how the data objects relate to each
other. Section \ref{secG-B:CovMatDataObjects} contains some perhaps
surprising examples on the strong impact of metric choice in the context
of covariance matrices as data objects\index{data objects!covariance matrices}.

Generally a good choice of metric is very situation dependent. For
example, in the case of Euclidean data objects, say $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d}$,
the standard \emph{Euclidean} $L^{2}$ \emph{distance\index{distance!Euclidean}}
(a special case of the $L^{p}$ norms defined in (\ref{eqB.D:LpNorms}))
\begin{equation}
\delta_{2}\left(\boldsymbol{x},\boldsymbol{y}\right)=\left\Vert \boldsymbol{x}-\boldsymbol{y}\right\Vert _{2}=\left(\sum_{j=1}^{d}\left(x_{j}-y_{j}\right)^{2}\right)^{1/2}\label{eqCHp5EqA:EuclideanDistanceL2}
\end{equation}
is often very useful. However, when it makes sense to think in terms
of polar coordinates, and the important variation happens in the angular
direction with mostly distracting noise in the radial direction, a
more useful metric can be the \emph{cosine distance\index{distance!cosine}}
\begin{equation}
\delta_{C}\left(\boldsymbol{x},\boldsymbol{y}\right)=1-\frac{2}{\pi}\cos^{-1}\left(\frac{\boldsymbol{x}^{t}\boldsymbol{y}}{\left\Vert \boldsymbol{x}\right\Vert _{2}\left\Vert \boldsymbol{y}\right\Vert _{2}}\right),\label{eqChpE.G:CosDist}
\end{equation}
which is driven only by the angle between $\boldsymbol{x}$ and $\boldsymbol{y}$
projected onto the unit sphere. The cosine distance is actually a
pseudo metric\index{metric!pseudo} because it does not distinguish
between vectors that are multiples of each other. When outliers are
a major concern, the $L^{1}$ \emph{\index{distance L1}distance}
\begin{equation}
\delta_{1}\left(\boldsymbol{x},\boldsymbol{y}\right)=\left\Vert \boldsymbol{x}-\boldsymbol{y}\right\Vert _{1}=\sum_{j=1}^{d}\left|x_{j}-y_{j}\right|\label{eqCHp5EqA:L1DistanceL}
\end{equation}
is another useful approach because of its natural tendency to down-weight
their influence\index{influence|(}. The $L^{1}$ distance is also
very useful when sparsity\index{sparsity} is desired because of a
tendency to give optimization results with many $0$ entries. In some
situations, a drawback of the $L^{1}$ distance is a lack of rotation
invariance\index{rotation invariance|(}, i.e. a strong dependence
on the coordinate system, as shown in Figure \ref{figE.D:Toy-MedianRotInv}.
In situations where the variables have a relationship determined by
a non-singular covariance matrix $\Sigma$, a unit free metric is
the \emph{Mahalonobis \index{distance!Mahalanobis}distance}
\[
\delta_{M}\left(\boldsymbol{x},\boldsymbol{y}\right)=\left(\left(\boldsymbol{x}-\boldsymbol{y}\right)^{t}\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{x}-\boldsymbol{y}\right)\right)^{1/2}.
\]
The Mahalonobis distance to the mean is the multivariate generalization
of the one dimensional notion of number of standard deviations from
the mean.

\section{Fréchet Centers In Metric Spaces\label{secE.1:Fr=0000E9chet-mean}}

An important notion of \emph{center} of a set of data objects $\left\{ \chi_{1},\cdots,\chi_{n}\right\} $
in an arbitrary metric space $S$ (with distance $\delta$) is the
\emph{Fréchet mean}\index{Fréchet mean|(},
\begin{equation}
\arg\min_{\chi\in S}n^{-1}\sum_{i=1}^{n}\delta\left(\chi,\chi_{i}\right)^{2},\label{eq:ChpEEqBFrechetMean}
\end{equation}
from \citet{Frechet1948Elements}. The factor of $n^{-1}$ does not
affect the minimizer, but provides convenient notation below. This
is a direct generalization of the standard sample mean $\overline{\boldsymbol{x}}=n^{-1}\sum_{i=1}^{n}\boldsymbol{x}_{i}$
in Euclidean space $\mathbb{\mathbb{R}}^{d}$, because it is straightforward
to show that $\overline{\boldsymbol{x}}$ is the solution of (\ref{eq:ChpEEqBFrechetMean})
in the case of Euclidean distance\index{distance!Euclidean} (\ref{eqCHp5EqA:EuclideanDistanceL2}).
Insight as to how the Fréchet mean works is given in Figure \ref{figE-A:ToyEgFrechetMean}.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterE/OODAbookChpEFigAToyL2FrechetMeanCMYK}

\caption{Toy two dimensional example illustrating the Fréchet mean. Left panel
shows a (far from optimal) candidate point as a blue x, with line
segments representing the distance to each data point. Right panel
shows the optimal Fréchet mean as a green x, with again line segments
showing this point is much closer to the data points than the blue
candidate on the left. These issues are quantitated in each panel
using the Fréchet average of squared $L^{2}$ distances. \label{figE-A:ToyEgFrechetMean}}
 
\end{figure}
A toy data set of $n=5$ data points in $\mathbb{\mathbb{R}}^{2}$
are shown as black circles in both panels. The blue x in the left
panel is a candidate choice of mean. The Fréchet criterion in (\ref{eq:ChpEEqBFrechetMean})
(multiplied by $n=5$) is the sum of the squared lengths of the blue
line segments, whose numerical value for this candidate is seen to
be 2.64. The operation of solving the Fréchet optimization problem
can be thought of as moving the candidate point to minimize the Fréchet
criterion. The solution is shown as the green x in the right panel
(actually just the familiar sample mean vector $\overline{\boldsymbol{x}}$
in this special case), based on an overall shorter collection of line
segments and the smaller criterion value of 0.67.

Depending on the choice of metric, many familiar notions of center
can be viewed as Fréchet means. For example on the real line $\mathbb{R}^{1}$,
the median can be written in the form (\ref{eq:ChpEEqBFrechetMean}),
by taking $\delta=\delta_{2}^{1/2}$ (note that it is straightforward
to show that the square root of any metric is again a metric). Furthermore
on $\mathbb{R}_{+}=\left\{ x\in\mathbb{R}:x>0\right\} $ the geometric\index{geometric mean}\index{mean!geometric}
and harmonic\index{harmonic mean}\index{mean!harmonic} means are
Fréchet means, based on $\delta\left(x,y\right)=\left|\log x-\log y\right|$
and $\delta\left(x,y\right)=\left|x^{-1}-y^{-1}\right|$ respectively. 

As noted in \citet{Koenker2006FrechetMedian}, for any distance $\delta$
an important variation is the \emph{Fréchet median\index{Fréchet median|(}},
\begin{equation}
\arg\min_{\chi}n^{-1}\sum_{i=1}^{n}\delta\left(\chi,\chi_{i}\right),\label{eqE.H:FrechetMedian}
\end{equation}
which differs from the Fréchet mean only by replacing the power of
$2$ by $1$. This and other variations of the Fréchet mean have been
extensively studied in the field of \emph{robust\index{robust methods|(}
statistics}, see e.g. \citet{HampelEtAl2011Robust}, \citet{HuberRonchetti2009RobustStat},
\citet{StaudteSheather1990Robust} and \citet{Clarke2018RobustnessBook},
because both the median, and also distances $\delta$ which down-weight
points farther away, have good properties in terms of reduced sensitivity
to outliers. See \citet{Fletcher2009GeometricMedian} for interesting
applications of this approach in the context of Riemannian manifolds.
Deeper discussion of robustness in OODA contexts can be found in Chapter
\ref{chapO:Robust-Methods}.

Direct comparison between the Fréchet mean and median is provided
in Figure \ref{figE.C:Toy-L2FrechetMedian}. This example features
the same toy data set as in Figure \ref{figE-A:ToyEgFrechetMean},
but this time the center-point is calculated as the Fréchet median
based on $\delta_{2}$. Again a poor candidate point is shown as the
blue x, with the optimal solution shown in green. The key difference
is that the sum of lengths of the line segments (not their squares)
make up the Fréchet sums that are shown.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterE/OODAbookChpEFigCToyL2FrechetMedianCMYK}

\caption{Toy example illustrating the Fréchet median, based on the distance
$\delta_{2}$, for the same data set as in Figure \ref{figE-A:ToyEgFrechetMean}.
Again a sub-optimal candidate appears on the left, with the Fréchet
median on the right. This estimate is more robust in the sense of
reduced sensitivity to the outlying point in the data set. \label{figE.C:Toy-L2FrechetMedian}}
\end{figure}
First note the substantial difference between the Fréchet mean and
median in this case. In particular, the improved robust\index{robust methods|)}ness
property of the Fréchet median can be seen in terms of the green x
being in the middle of the convex hull of the 4 nearest points in
Figure \ref{figE.C:Toy-L2FrechetMedian}, while the outlier is far
enough away to pull it outside in Figure \ref{figE-A:ToyEgFrechetMean},
because it has much stronger influence when the distances are squared.

The impact of metric choice is again illustrated in Figure \ref{figE.B:Toy-L1FrechetMedian}.
This is also the same toy data as in Figures \ref{figE-A:ToyEgFrechetMean}
and \ref{figE.C:Toy-L2FrechetMedian}, but now the distance used in
the Fréchet median is the $L^{1}$ distance $\delta_{1}$. Because
this distance is the sum of absolute distances in each coordinate
direction, each is now represented as a horizontal and a vertical
line.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterE/OODAbookChpEFigBToyL1FrechetMedianCMYK}

\caption{Study of the impact of metric choice using the same toy data from
Figures \ref{figE-A:ToyEgFrechetMean} and \ref{figE.C:Toy-L2FrechetMedian}.
Again the left panel shows the Fréchet sum, but now with $\delta_{1}$,
for a poor candidate, and the right shows the optimum, i.e. the $L^{1}$
Fréchet median. This version of center also demonstrates reduced influence\index{influence|)}
of outliers. \label{figE.B:Toy-L1FrechetMedian}}
\end{figure}
This variation of the Fréchet median reveals some perhaps surprising
differences with the $\delta_{2}$ version. In particular, the solution
here is the point-wise median (i.e. this can be computed by simply
taking the median of each coordinate of the data vectors), which is
not true for the $\delta_{2}$ version of the Fréchet median. Thus
the $\delta_{1}$ variation is easier to compute, while an iterative
(but rather fast even in high dimensions) computation is needed in
the case of $\delta_{2}$. 

Another basis for the comparison of these notions of center is \emph{rotation
invariance}, as illustrated in Figure \ref{figE.D:Toy-MedianRotInv}.
The panels show different rotations of the same toy data set in $\mathbb{R}^{2}$,
together with both the $\delta_{1}$ and $\delta_{2}$ Fréchet medians.
Because the distance $\delta_{2}$ is rotation invariant, the greenish
blue x showing that Fréchet median has the same relative position
when the data are rotated. But because the distance $\delta_{1}$
is not rotation invariant, this is not true for that version, shown
as the red + sign. This example is deliberately constructed to show
that the $\delta_{1}$ Fréchet median (i.e. the coordinate-wise median,
which is here seen to be a naive choice of ``multivariate median'')
can be viewed as a poor notion of center, in the sense that it actually
lies on the boundary of the convex hull of the data.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterE/OODAbookChpEFigDToyMedianRotInvCMYK}

\caption{Two dimensional toy example studying rotation invariance\index{rotation invariance|)}.
Data in the right panel are a 45 degree rotation of the data on the
left. In each case the $\delta_{2}$ based Fréchet median is shown
as the blue-green x, with the $\delta_{1}$ Fréchet median appearing
as the red plus sign. This shows the latter is not rotation invariant.
\label{figE.D:Toy-MedianRotInv}}
\end{figure}
Despite these clear differences, note that both of these medians are
direct generalizations of the familiar one dimensional median in $\mathbb{R}^{1}$
in the sense that they both reduce to that quantity in the case $d=1$.
As discussed in Section \ref{secO.C:OtherRobust}, actually there
are a number of other quite different notions of center in $\mathbb{R}^{d}$
which also reduce to the median in the case $d=1$.

Neither the Fréchet mean nor median are guaranteed to be unique. This
is typically handled by working with \emph{mean set\index{mean!set}s}
(\emph{median sets}, resp.). For example, on $\mathbb{R}^{1}$ when
$n$ is even, the interval between the central two data points is
the median set. In that situation, and in others, sometimes a representative
is chosen for the set, e.g. the midpoint of the interval in the $\mathbb{R}^{1}$
median case. The extreme case of non-uniqueness is the discrete metric
$\delta_{D}$, with respect to which both the Fréchet mean and median
sets are the full data space, and thus are of no practical use.

From the mathematical statistical perspective an important issue is
just what is being estimated by all of these notions of center. Useful
models for this are based on the notion of randomness and probability
theory. The convention here and in the rest of this book is that a
tilde over notation indicates a random quantity as shown in Table
\ref{tab7-1:TildeNotation}.

\begin{table}[H]
\begin{tabular}{|c|c|}
\hline 
Type of Quantity & Notation for Random Version\tabularnewline
\hline 
\hline 
random number (scalar) & $\widetilde{x}$\tabularnewline
\hline 
random vector & $\widetilde{\boldsymbol{x}}$\tabularnewline
\hline 
matrix of random variables & $\widetilde{\boldsymbol{X}}$\tabularnewline
\hline 
element of arbitrary space & $\widetilde{\chi}$\tabularnewline
\hline 
\end{tabular}

\caption{\index{notation for random quantities}Notation for random quantities.
\label{tab7-1:TildeNotation}}

\end{table}
When the data objects are thought of as being a \emph{random sample\index{random sample}}
(i.e. drawn independently) from a probability distribution on the
object space\index{object space}, then appropriate notions of population
center are defined using Fréchet criteria based on expected values.
In particular, sample Fréchet means as defined in (\ref{eq:ChpEEqBFrechetMean})
can be usefully treated as estimates of the corresponding theoretical
(i.e. population) Fréchet means 
\begin{equation}
\arg\min_{\chi\in S}E\delta\left(\chi,\widetilde{\chi}\right)^{2},\label{eqChpE-F:FrechetMin}
\end{equation}
where $\widetilde{\chi}$ is a random variable with the population
probability distribution. Similarly for Fréchet medians. A good overview
of many results studying various types of asymptotic convergence of
estimates to their population versions, in the challenging contexts
of data objects on manifolds and stratified spaces\index{stratified space},
can be found in \citet{PatrangenaruEllingson2016}. That book also
provides a good overview of affine and projective shape manifolds.
For some early related discussion of statistics on stratified spaces
see \cite{Bhatetal13}. Some interesting more recent results can be
found in \citet{HuckemannEltzner2017BackwardsAsys} who derive asymptotics
for backwards methods\index{backwards methods} as developed in Chapter
\ref{chapG:Manifold-Data-Analysis}, and \citet{HuckemannHotz2016NonparametricAsys}
who explore nonparametric\index{nonparametric} methods on manifolds.
See \citet{HotzEtal2013StickyMeans} for a particularly unusual limiting
distribution theory for the Fréchet mean in a stratified space, involving
the nonstandard notion of \emph{stickiness\index{stickiness}}.

There are several synonyms of Fréchet mean. These include \emph{\index{Riemannian barycenter}Riemannian
barycenter} (really the center of mass in physics, when using the
Euclidean distance $\delta_{2}$), and \emph{geodesic mean} (when
$\delta$ is a geodesic \index{distance!geodesic}distance, say on
a curved manifold, as discussed in Chapter \ref{chapG:Manifold-Data-Analysis}).
The term Karcher mean\index{Karcher mean} is sometimes also used,
with a common distinction being that it refers to only a local minimizer,
instead of a global minimizer as in (\ref{eqChpE-F:FrechetMin}).
However this terminology is not universally accepted, see \citet{Karcher2014NameObjection}.
Similarly there are synonyms for the Fréchet median\index{Fréchet median|)}
including \emph{\index{geometric median}geometric median}, \emph{\index{spatial median}spatial
median} and $L^{1}$ \emph{M estimate.}\index{$L^{1}$ M-estimate}

Also useful is to define a \emph{weighted \index{Fréchet mean!weighted}Fréchet
mean} in an arbitrary metric space $S$ with fixed weights $w_{i}$,
$i=1,\cdots,n$, by solving 
\[
\arg\underset{\chi\in S}{\min}\sum_{i=1}^{n}w_{i}\delta(\chi,\chi_{i})^{2}
\]
where $\sum_{i=1}^{n}w_{i}$=1. The Fréchet mean given in (\ref{eq:ChpEEqBFrechetMean})
has equal weights $w_{1}=1$, $i=1,\cdots,n$. An example where the
weighted Fréchet mean is useful is the construction of an \emph{\index{interpolation path|(}interpolation
path} between two objects $\chi_{1}$ and $\chi_{2}$, where the weights
depend on a parameter $\alpha\in[0,1]$:
\begin{equation}
\arg\underset{\chi\in S}{\min}(1-\alpha)\delta(\chi,\chi_{1})^{2}+\alpha\delta(\chi,\chi_{2})^{2}.\label{eqE.I:FrechetInterp}
\end{equation}
Values $\alpha=0$ and $\alpha=1$ correspond to the two objects at
either end of the path, and the Fréchet mean\index{Fréchet mean|)}
of the two objects is in the middle of the path given by $\alpha=\frac{1}{2}$.
Insightful examples of interpolation paths on the manifold of covariance
matrices are given in Figure \ref{interp}, where quite different
behavior is observed for each choice of distance. For certain types
of distances (intrinsic\index{intrinsic} \index{distance!intrinsic}distances)
the interpolation path also provides a \emph{\index{minimal geodesic path}minimal
geodesic path} between points on a manifold\index{manifold|)}.

Values of $\alpha$ in (\ref{eqE.I:FrechetInterp}) which are outside
the interval $[0,1]$ lead to an \emph{\index{extrapolation path}extrapolation}
path. In particular if we choose $\alpha<0$ we extrapolate beyond
$\chi_{1}$ in an extension of the path from $\chi_{2}$ to $\chi_{1}$,
and if $\alpha>1$ we extrapolate beyond $\chi_{2}$ in an extension
of the path from $\chi_{1}$ to $\chi_{2}$. 

\section{Multi-Dimensional Scaling For Object Representation\label{secE.2:Multi-dimensional-scaling}}

As with notions of center, there are many ways of quantitating variation
about the center, based solely on metrics. A very simple one is the
\emph{Fréchet variance\index{Fréchet variance}} which is just the
minimum value attained in (\ref{eq:ChpEEqBFrechetMean}).

But generally more useful for data analytic tasks, such as those discussed
in Sections \ref{secB2:Data-Visualization} - \ref{secB4:Further-Statistical-Tasks},
are distance based analogs of PCA. One approach to this is \emph{Multi-Dimensional
Scaling\index{Multi-Dimensional Scaling|(}} (MDS). MDS has been very
popular in the psychometrics literature, and at least the nomenclature
is usually attributed to \citet{Torgerson1952MDS,Torgerson1958MDSbook}
and \citet{Gower1966DS}. However, the underlying mathematics is substantially
older, see \citet{Eckart1936Approximation} and \citet{Young1938MDS}.
MDS is also called \index{Principal Coordinates Analysis}Principal
Coordinates Analysis (abbreviated PCoA) in the biological literature.

In its simplest form, MDS starts with a set of $n$ data objects,
with a known set of pairwise distances between them summarized as
a distance matrix $D$ as in (\ref{eq:ChpEeqC:DistanceMatrix}), and
seeks to represent the objects as a set of points $\boldsymbol{x}_{1},\cdots,\boldsymbol{x}_{n}\in\mathbb{R}^{d}$
for some $d$, in such a way that the Euclidean distances $\delta_{2}\left(\boldsymbol{x}_{i},\boldsymbol{x}_{j}\right)$
approximate the elements of $D$ as well as possible, in various senses.
When the input distance matrix $D$ is itself composed of pairwise
Euclidean distances $\delta_{2}$, typical basic algorithms return
$\boldsymbol{x}_{i}$ as the vector of the first $d$ PC scores. In
that sense MDS extends PCA to cases where only distances are known.
A toy example is shown in Figure \ref{figE.E:DistColsDataObjectsPart1}.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterE/OODAbookChpEFigEDistColsCMYK1}

\caption{Toy 2-d Gaussian example, illustrating MDS. Raw data is shown using
+ signs in the left panel. Right panel shows the corresponding MDS
scores as circles (with the same coloring), calculated from the Euclidean
distance matrix, which are essentially the same as the PC scores for
the data on the left.\label{figE.E:DistColsDataObjectsPart1}}
\end{figure}
The raw data in the left panel are an elongated Gaussian point cloud,
colored along the major axis of elongation. To construct the plot
in the right panel, the $\delta_{2}$ distance matrix $D$ was computed,
and then the corresponding classical (\citet{Torgerson1952MDS,Torgerson1958MDSbook})
MDS coordinates were plotted in the right panel. Note this looks much
like what would be expected from PCA scores for this data set, which
is consistent with the above discussion. As discussed in Section \ref{secF.C:Scatterplot-Matrix-Views}
for PCA scores, each of these coordinates is only determined up to
an arbitrary sign flip. In particular the colors suggest a reversal
along the long axis. 

Figure \ref{figE.H:EquatorialDataMDS} studies how the choice of metric
is critical to MDS. The data are called ``Equatorial'' because they
are distributed roughly along the equator of $S^{2}$ as shown in
the left panel. That is done using a spread Gaussian distribution
to determine the longitude together with latitudes from a tight Gaussian.
One approach to the distance between data objects is the Euclidean
distance in the embedding space $\mathbb{R}^{3}$. The MDS scores
with respect to that Euclidean distance matrix are shown in the center
panel. It is not surprising that the distribution looks much the same
as the projections of the data into the plane determined by the equator,
except they are centered at the mean. This follows from the fact that
MDS based on Euclidean distance is essentially PCA (of the data in
$\mathbb{R}^{3}$). A much different geodesic \index{distance!geodesic}distance
is used in the right panel. There distance is measured by arc lengths
along the surface. That distance naturally follows the spread of Gaussian
points, so that the distance between points in opposite tails is nearly
$2\pi$. That results in a far different set of scores, giving a more
intuitively useful view of the distribution of the data, that is much
more indicative of how they were generated.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterE/OODAbookChpEFigHSphereHorseshoe}

\caption{Equatorial data on $S^{2}$, for illustrating impact of metric on
MDS. Euclidean distance MDS is shown in the center panel, exhibiting
a strong horseshoe shape. The right panel displays the MDS using geodesic
distance (along the surface of the sphere), which eliminates the metric
induced distortion. \label{figE.H:EquatorialDataMDS}}
\end{figure}
As noted in Section 14.5.2 of \citet{MardiaKentBibby1979MultivarAnal},
point cloud shapes such as those in the center panel arise when we
can measure distance accurately when objects are close together, but
not when they are far apart. This situation occurs surprisingly often
in practice. That general pattern was given the name \emph{horseshoe}
by \citet{Kendall1970Horseshoe}. \citet{DiaconisHolmes2008Horseshoes}
pointed out that horseshoe effect\index{horseshoe effect|)}s tend
to occur in situations where distant points are less far away than
would be expected from linear extrapolation of the set of local distances.
Figure \ref{figE.H:EquatorialDataMDS} gives a canonical example of
this. Another context where far away points have distances less than
expected by extrapolation appears in the study of micro-biomes as
detailed in \citet{Morton2017Horseshoe}. They coined the term \emph{distance\index{distance!saturation}
saturation} to describe the situation where far distances are less
than expected by linear extrapolation. Horseshoe effects, together
with a higher order version, are discussed again in Section \ref{subsecH.A.A:ExampleShiftedBeta}.

A real data example that may be due to this horseshoe effect can be
seen in the Drug Discovery data\index{Drug Discovery data} shown
in Figure \ref{figD-K:DrugVDnonBinaryPCA}. This suggests that an
alternate representation of that data set could be useful. 

There are many generalizations of the basic MDS idea. An important
case in the psychometric literature extends the input of a distance
matrix to merely a \emph{dissimilarity\index{dissimilarity!matrix}
matrix}, whose entries are not required to satisfy the metric properties
such as the triangle inequality, see \citet{Kruskal1964MDS} for good
discussion. More overviews of this large area can be found in \citet{CoxCox2000MDSbook},
\citet{BorgGroenen2005MDSbook}, \citet{Buja2008VisMDS} and Chapter
3 of \citet{Zhai2016PCAphylo}.

For OODA, an important generalization of MDS is to replace the embedding
into Euclidean space $\mathbb{R}^{d}$, with embeddings into curved
spaces. This is useful for data spaces that are strongly curved, such
as phylogenetic tree\index{phylogenetic tree}\index{tree} space,
as discussed in Section \ref{subsecI-A-C:Phylogenetics}. 

A weakness of distance and dissimilarity matrix based methods, is
that they only give analyses of the data set at hand, and are challenging
to extend to contexts requiring generalization to additional data,
such as the classification / discrimination tasks considered in Chapter
\ref{chapJ:Classification}. An interesting approach to this is the
\emph{\index{out of sample MDS}out of sample MDS} ideas of \citet{TrossetPriebe2008OutSampMDS}.
See Section 3.4 of \citet{Zhai2016PCAphylo} for a more detailed overview
of out of sample MDS.

As noted above, from the perspective of modes of variation\index{modes of variation!MDS}
from Section \ref{subsecBA.1.1:FormalizingModes}, a perhaps more
serious weakness of MDS type analyses is that, while they can provide
useful and insightful \index{scores!MDS}scores, they fail to provide
an analog of loadings\index{loadings} that would be needed to give
full modes of variation. This means that understanding of the drivers
of modes, e.g. as done with loading plots in Figure \ref{figB-X:PanCanDWDLoadings}
for the Pan Cancer data\index{Pan Cancer data}, is not available
from this approach.

MDS has also been used as a way of representing landmark\index{landmark|(}
shapes, using the set of all distances between pairs of landmarks
to represent the shape (up to a reflection). This approach has been
used by \citet{Bandpatr05}, \citet{Dryden2008} and \citet{Bhattacharya08}.
Important choices to be made are: how to average the distance matrices
and how to project back to the space of configurations.

\citet{Duin2005Dissimilarity} survey another way of handling data
presented only in terms of a distance matrix $D$ from the machine
learning\index{machine learning} literature. That is to simply take
the columns of $D$ as data objects, i.e. to let the rows of $D$
define the features (using the object - feature terminology from Section
\ref{secBA.1:Foundations_Terminology}). In other words the trait\index{trait}s
become the distances to the other objects. The impact of this approach
on the relationship between this choice of data objects is studied
in Figure \ref{figE.E:DistColsDataObjectsPart2}, where direct understanding
comes from using the same data (and colors for keeping track of individual
data points) as in Figure \ref{figE.E:DistColsDataObjectsPart1}.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterE/OODAbookChpEFigEDistColsCMYK2}

\caption{PCA of data using columns of the Euclidean distance matrix as data
objects, for the same data and colors as in Figure \ref{figE.E:DistColsDataObjectsPart1}.
Shows quite strong distortion of the original relationships between
data objects. \label{figE.E:DistColsDataObjectsPart2}}
\end{figure}
This PCA scores scatterplot matrix\index{scatterplot matrix}\index{scores!scatterplot}
is in the same format as shown several times in Chapter \ref{chapBB:ExpConfAnal}.
The colors in the upper left panel, indicate that the first PC scores\index{scores|)}
are much as expected. However, the PC1 vs. PC2 scatterplot in the
middle panel of the top row shows structure that appears to be the
same horse-shoe effect studied in Figure \ref{figE.H:EquatorialDataMDS}.
In particular, points at both the magenta and red ends are closer
in this data object representation, than they are in either panel
of Figure \ref{figE.E:DistColsDataObjectsPart1}, which generates
a very strong distortion in the relationships between\index{data objects!relationships between|)}
the data objects. 

Insight into this distortion again comes from the ideas of \citet{DiaconisHolmes2008Horseshoes}
as illustrated in Figure \ref{figE.I:Horseshoe}. Think of a data
set of numbers equally spaced on $[0,1]$. The dashed (dot-dashed)
line segments represent distance matrix column data objects for the
points $0.4$ ($0.8$ respectively), which have distance $\delta=0.4$
between them. The $L^{1}$ distance between the piecewise linear data
object representations is the cyan shaded area. To first order, this
is $\delta\times1=\delta.$ However, a closer look shows the area
is more precisely $\delta-\delta^{2}/2$ (half the area shown as a
small square box needs to be subtracted). A very similar calculation
shows the corresponding $L^{2}$ distance is approximately $\sqrt{\delta^{2}-\delta^{3}/2}.$
Both of these have the property that the dependence on $\delta$ is
sub-linear, which fits exactly into the \citet{DiaconisHolmes2008Horseshoes}
explanation of horseshoe effect\index{horseshoe effect|)}s. For most
purposes these issues seem like a serious weakness of the practice
of using columns of a distance matrix\index{distance matrix|)} as
data objects. However in the spirit of kernel methods described in
Section \ref{secJ-B:KernelMethods} there may conceivably be situations
where this type of distortion of the relative positions of the data
objects could be useful. 

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterE/OODAbookChpEFigIHorseshoeCMYK1}

\caption{Reveals reason for horseshoe distortion, caused by using columns of
a distance matrix as data objects, seen in Figure \ref{figE.E:DistColsDataObjectsPart2}.
Half the area of the box shows how typical distances are sublinear.
\label{figE.I:Horseshoe} }
\end{figure}

\section{Important Distance Examples\label{secI-C:ImportantDists} }

\subsection{Conventional Norms\label{subsecI-C-A:ConventNorms}}

As noted at the beginning of this chapter, a broadly useful family
of distances on $\mathbb{R}^{d}$ are based on the $L^{p}$ norms,
defined in (\ref{eqB.D:LpNorms}). The $L^{p}$ norms are also directly
extendable to metrics on function spaces as 
\[
\delta_{p}\left(f,g\right)=\left\Vert f-g\right\Vert _{p}=\left(\int_{-\infty}^{\infty}\left|f(x)-g(x)\right|^{p}\right)^{1/p}.
\]
There are a number of other useful norms on the space of matrices,
including the \emph{Frobenius norm\index{Frobenius norm|(}} which
is the square root of the sum of squared matrix entries:

\begin{equation}
\| \mathbf{X} \|_F = \sqrt{ {\rm trace} (\mathbf{X}^t\mathbf{X}) } = \left( \sum_i \sum_j x_{ij}^2 \right)^{1/2}  \label{eq:Frobeniusnorm}
\end{equation}which is also sometimes known as the Euclidean norm or Hilbert-Schmidt
norm\index{Hilbert-Schmidt norm}. The Frobenius norm is used in Procrustes
analysis in Sections \ref{subsecI-C-C:Procrustes} and \ref{subsec:Generalized-Procrustes-analysis},
in the analysis of covariance matrix data objects in Section \ref{subsec:Covariance-matrix-distances},
and in the analysis of networks as data objects in Section \ref{secI-B:NetworkData}. 

Other matrix norms include the \emph{operator norm\index{operator norm}}
(the largest singular value), the \emph{\index{one norm}one norm}
(maximum absolute column sum of the matrix), the \emph{\index{infinity norm}infinity
norm} (maximum absolute column row of the matrix) and the \emph{\index{maximum norm}maximum
norm} (maximum modulus of all elements). Finally, the \emph{\index{entrywise norm}entrywise
norm} of the $m \times n$ matrix $\mathbf{X}$ is the $L^p$ norm
(\ref{eqB.D:LpNorms}) of the vectorize\index{vectorize}d version
of $\mathbf{X}$, i.e. 

\begin{equation}
\| \mathbf{X} \|_{EW} = \left( \sum_{i=1}^m \sum_{j=1}^n |\mathbf{X}_{ij}|^p \right)^{1/p}. \label{eq:EWnorm}
\end{equation} For $p=2$ the entrywise norm is the Frobenius norm. Entrywise norms
are particularly appropriate for computation with high-dimensional,
sparse matrices. An example using several of these matrix norms is
given in Section \ref{subsec:corpuslingustics} when comparing high-dimensional,
sparse networks in the context of natural language processing. 

\subsection{Wasserstein Distances\label{subsecI-C-B:Wasserstein}}

A particularly useful metric on the space of probability measures
is the \emph{Wasserstein} (i.e. Kantorovich--Rubinstein) \emph{distance\index{distance!Wasserstein}}.
In the case of discrete measures this is also called the \emph{Earth
Mover's \index{distance!Earth Mover's}distance} (usually in Computer
Science). That name nicely expresses the intuition that this distance
essentially quantifies how much probability mass needs to be moved
to get from one measure to the other. For this reason it is not surprising
that computation of Wasserstein distances in general involves optimal
transport formulations. In theoretical studies of the bootstrap a
popular criterion is \emph{Mallows \index{distance!Mallows}distance},
proposed in \citet{Mallows1972Distance} and shown to be equivalent
to the Wasserstein Distance by \citet{LevinaBickel2001EarthMovers}.
A comprehensive resource for many important aspects of Wasserstein
distances is \citet{PanaretosZemel2020invitation}. Theorem 1.5.1
of that book shows that in the case of probability measures on the
real line, the Mallows distance is just the $L^{2}$ norm of the \index{quantile}quantile
functions (defined in Section \ref{subsecBA.2.1:ProbDistns}). Hence
the good properties of the Wasserstein metric\index{metric!Wasserstein|(}
are very consistent with the conclusion from Figures \ref{figB-V:ProbDensObjects}
and \ref{figB-W:QuantFuncObjects} that quantiles are frequently useful
representatives of probability distributions as \index{data objects!probability distributions}
data objects. The ideas of this section suggest that another useful
approach to tasks like PCA \index{Principal Component Analysis} in
that case is MDS with respect to the Wasserstein distance. 

As discussed in Section \ref{secBA.2:OverviewObjFeat}, the optimal
transport aspect of the Wasserstein distance accomplishes in a more
rigorous way the goals of the \index{Visual Error Criterion}Visual
Error Criterion of \citet{MarronTsybakov1995VisErr}. This is seen
in Figure \ref{figE.G:Wasserstein}, which shows the main example
in that paper (which was motivated by one in \citet{KooperbergStone1991Logspline}).
The example consists of a target curve (shown as the solid line type)
and two candidate estimates (dashed and dotted). The target curve
has two distinct bumps, one quite thin, and the other broader. Especially
in situations where bumps are critical, Estimate 2 seems preferable
as it has similar bumps, although the thin bump is not quite in the
right location. On the other hand Estimate 1 has smoothed out that
perhaps important feature. Now if effectiveness of the curve estimation
is measured using the $L^{2}$ norm (the same holds for the $L^{1}$
norm), Estimate 1 is closer to the target. This is because the differences
between the curves are measured vertically, so Estimate 2 is penalized
twice for having its bump in the wrong location, while Estimate 1
only pays the peak penalty once. However, measuring discrepancy using
the squared error version of the Wasserstein metric gives the opposite
result. In that sense Estimate 2 is closer to the target than Estimate
1.

\begin{figure}[H]
\includegraphics[scale=0.7]{ChapterE/OODAbookChpEFigGWasserstein}

\caption{Kooperberg-Stone example showing curve estimation context that is
conceptually challenging to the $L^{2}$ norm, but nicely quantified
using the $W^{2}$ version of the Wasserstein metric. \label{figE.G:Wasserstein}}
\end{figure}
A useful application of MDS\index{Multi-Dimensional Scaling} based
on the Wasserstein metric\index{metric|)}, in the context of FDA\index{functional data analysis}
curve registration\index{curve registration!Wasserstein metric} of
the Shifted Betas data\index{Shifted Betas data} is shown in Figure
\ref{figH.N:ToyBadWasser}.

\subsection{Procrustes Distances\label{subsecI-C-C:Procrustes}}

\index{data objects!landmark configurations}\index{distance!Procrustes|(}

An introduction to the shape of landmark configurations as data objects was given in Section \ref{subsec1A.B.B:ShapeRepresentations}. Much more discussion appears in Section \ref{sec:Landmark-Shapes}. We now consider suitable distances that measure differences 
between the shapes of two landmark configurations. 

Consider a configuration of $k$ landmarks in $m$ dimensions represented by a $k \times m$ matrix $\mathbf{X}$, where $k > m$. Usually we have $m=2$ or $m=3$.  Procrustes analysis involves matching a set of such configurations  using translation, rotation and possibly scale using least squares techniques.  Procrustes analysis was developed initially for applications of factor analysis in psychology. The name  was given by \cite{Hurlcatt62} following the loose analogy  with Greek mythology where Procrustes would stretch his victims to fit exactly to a bed if too short, or chop off their limbs if they were too tall.   The technique can  be traced back to \cite{Boas05} and \cite{Mosier39}. Other important references  include  \cite{Gower75},  \cite{Kendall1984Shape} and  \cite{Goodall91}, and a summary is given by  \cite[Chapter 7]{DrydenMardia2016}. 

{\it Ordinary Procrustes Analysis} \index{Procrustes analysis!ordinary} is where one configuration is fitted to another by minimizing the sum of squared distances of landmarks between  the configurations, using translation, rotation and possibly scale. 

\subsubsection*{Example: Digit 3 Data} 
The {\it Digit 3 data} are \index{Digit 3 Data}
available in the {\tt shapes} package in R \citep{Dryden2019RShapes} and further discussed in Sections 
\ref{subsec:Generalized-Procrustes-analysis}, \ref{sec:Landmark-Shapes} and \ref{subsec:PCA-Digit3}. 
The landmarks on each digit were located by hand on images of envelopes containing British postcodes by  \cite{Anderson97},  
and there are $n=30$ handwritten digits each with $k=13$ landmarks in $m=2$ dimensions. 

In Figure \ref{Digit3} we see two of the handwritten digit 3s in red and blue, 
where each digit has $k=13$ landmarks in $m=2$ dimensions. 
The left-hand plot shows the unaligned data. The 
right-hand plot shows the two digits after they have been centered and the blue digit has been rotated and rescaled to match the red digit as closely as possible using Ordinary 
Procrustes Analysis. 


\begin{figure}[H] \begin{center} \includegraphics[trim=0cm 5cm 0cm 5cm, clip=true, width=10cm]{./ChapterE/digit3.png} \end{center} \caption{An example of two handwritten digits from the Digit 3 data. 
The left-hand plot shows the original data. The 
right-hand plot shows the digits after they have been centered and in addition the blue digit has been rotated and rescaled to match the red digit as closely as possible.}\label{Digit3} \end{figure}

\subsubsection*{Example: DNA Molecule Data} 
The {\it DNA Molecule data} form \index{DNA Molecule data} 
part of the study by  \cite{Drydenetal17} and are further discussed in Sections 
\ref{subsec:Generalized-Procrustes-analysis}, \ref{sec:Landmark-Shapes} and \ref{subsec:PCA-DNA}. 
The DNA Molecule dataset consists of $n=50$ DNA molecules, each with $k=22$ landmarks in $m=3$ dimensions. 
More specifically the dataset is the TFC molecule in \citet{Drydenetal17} (a type of damaged DNA molecule) provided 
by Charlie Laughton and the observations have been temporally thinned to each 50th observation. The landmarks are located at the phosphorous atoms in the molecule.

In Figure \ref{DNAfig1} we see two of the DNA molecules.
Each landmark is marked by 
small red and blue spheres in the left-hand plot for each molecule. 
Lines are joined along each of the two strands, and also between the 
strands to indicate the base pairs. The familiar helix shape of DNA 
can be seen. In this perspective view of the DNA molecule the spheres are smaller for atoms further away and larger for atoms that are closer to the viewer.  

In the 
right-hand plot the blue molecule is translated, rotated and rescaled to minimize the sum of squared Euclidean distances between pairs of atoms. We say that the blue DNA molecule has been Procrustes rotated onto the  red DNA molecule.  

\begin{figure}[H] \begin{center} \includegraphics[height=7cm]{./ChapterE/DNA-TFC-pair5.png}  \includegraphics[trim=0cm 1cm 0cm 0cm, clip=true,height=7cm]{./ChapterE/dna-pair-withscale.png} \end{center} \caption{An example of two molecules in red and blue from the DNA Molecule data, each consisting of $k=22$ landmarks in $m=3$ dimensions. The left-hand plot shows a view of the two DNA molecules 
in arbitrary positions. In the right-hand plot   the molecules have been matched using Ordinary Procrustes analysis with the blue molecule translated, rotated and scaled to match the red molecule.}\label{DNAfig1} \end{figure}


We can measure the distance between two shapes using Procrustes analysis. 
We first transform the configurations to pre-shapes by centering and rescaling to 
unit norm (see (\ref{eq:preshape})), and then optimize over rotation.  
The {\it rotation group} $SO(m)$ is the set of all rotation matrix operators on $\mathbb{R}^{m}$, which is the set of all $m\times m$  matrices with orthonormal columns  (simply summarized as $\boldsymbol{U}^{t}\boldsymbol{U}=\boldsymbol{I}_{d}$) and determinant $1$.
The {\it partial Procrustes distance} is a shape distance given by  \begin{equation} \delta_{PP}( \mathbf{X}_1,\mathbf{X}_2) = \min_{\mathbf{\Gamma} \in SO(m)} \| \mathbf{Z}_2 - \mathbf{Z}_1 \mathbf{\Gamma} \|_F ,  \label{dP} \end{equation} where  $\mathbf{Z}_j = \mathbf{C}  \mathbf{X}_j / \| \mathbf{C}  \mathbf{X}_j\|_F, \; j=1,2$ are called pre-shapes;  
\begin{equation} 
\mathbf{C} =  \mathbf{I}_k -  \mathbf{J}_{k,1}  \mathbf{J}_{k,1}^t/k \label{centeringmatrix}
\end{equation}
is the centering matrix; \index{matrix!centering}  and $\| \mathbf{X} \|_F$ is the 
Frobenius norm (\ref{eq:Frobeniusnorm}). Other alternative 
distances include the full Procrustes distance (\ref{dF}) and Riemannian shape distance (\ref{rho}) which are described in Section \ref{sec:Landmark-Shapes}.  
For small shape changes there is little difference between these choices of shape distance, and the distances between the 
pairs of Digit 3s (a large distance) and the pairs of DNA molecules (a small distance) are given in Table \ref{Tab:Shapedist}. 



\subsection{Generalized Procrustes Analysis\label{subsec:Generalized-Procrustes-analysis}}

In order to compute the mean from a sample of shapes the most commonly used method in practice  is Generalized Procrustes Analysis, e.g. see  \cite{Gower75},  \cite{Goodall91}  and  \cite{DrydenMardia2016}. \index{Procrustes analysis!generalized}
The method of Generalized Procrustes Analysis involves translating, rotating  and rescaling the configurations  relative to each other so as to minimize a total sum of squares \begin{equation}   \sum_{i=1}^n  { \| (\beta_i \mathbf{X}_i \mathbf{\Gamma}_i + 1_k \boldsymbol{\gamma}_i^t)  - \boldsymbol{\mu} \|_F^2 }   \end{equation}  with respect to scale $\beta_i$, rotation $\mathbf{\Gamma}_i$, translation $\boldsymbol{\gamma}_i, i=1,\ldots,n$, and an overall mean $\boldsymbol{\mu}$, subject to  an overall size constraint (which avoids everything shrinking to the origin), such as \begin{equation}   \sum_{i=1}^n S^2(\beta_i  \mathbf{X}_i \mathbf{\Gamma}_i + \mathbf{1}_k \boldsymbol{\gamma}_i^{\rm T}) = \sum_{i=1}^n S^2( \mathbf{X}_i ) . \label{constraint2} \end{equation} where 
\begin{equation}
S(\mathbf{X}) = \| \mathbf{C} \mathbf{X} \|_F  \label{eq:centroid-size}
\end{equation} 
is the centroid size of $\mathbf{X}$. \index{size!centroid} \index{size} 
 If we write $\hat{\boldsymbol{\gamma}}_i, \hat{\mathbf{\Gamma}}_i, \hat\beta_i$ for the optimal translation, rotation, scale then the Procrustes mean configuration   after matching is     \begin{equation}   \hat{\boldsymbol{\mu}} = \frac{1}{n} \sum_{i=1}^n  \hat\beta_i \mathbf{X}_i \hat{\mathbf{\Gamma}}_i + 1_k \hat{\boldsymbol{\gamma}}_i^t .   \end{equation} 
The translations $\hat{\boldsymbol{\gamma}}_i$ are simply obtained by object centering. However, the rotations  and scales must be obtained by an iterative algorithm. A summary of some algorithms is given by  \cite[Section 7.4]{DrydenMardia2016} and an implementation  in R is given by the function {\tt procGPA} in the shapes package  \citep{Dryden2019RShapes}. 

The shape of $\hat{\boldsymbol{\mu}}$ is also equal to the Fr\'echet mean (\ref{eq:ChpEEqBFrechetMean}) with respect to the full Procrustes distance (where $\delta_{FP}$ is given in (\ref{dF})), which is obtained from     $$ \hat{\boldsymbol\mu} = {\rm arg} \min_{\boldsymbol{\mu}}  \sum_{i=1}^n \delta_{FP}( \mathbf{X}_i,  \boldsymbol{\mu})^2.$$  

Methods for studying variation in shape are studied in Section \ref{sec:shapetangentspace}. In particular, after alignment by Generalized Procrustes Analysis, Principal Components
Analysis and other techniques can be carried out to explore the modes of 
shape variation.

\subsubsection*{Case Study: Digit 3 Data}
An example of Generalized Procrustes Analysis for the full Digit 3 data (introduced in Section \ref{subsecI-C-C:Procrustes}) is given in Figure  \index{Digit 3 data} 
\ref{Digit3-full}. Here there are $n=30$  handwritten digit 3s with $k=13$ landmarks 
in $m=2$ dimensions. The original digits are viewed in the left-hand plot in red.  
In the right-hand plot we seen  the Procrustes registered digits in green, which are obtained using Generalized Procrustes Analysis. In both plots the shapes of the digits are the same, but in the right hand plot the unimportant rotations, translations and scales have been removed. 

\begin{figure}[H] \begin{center} \includegraphics[trim=0cm 5cm 0cm 5cm, clip=true,width=10cm]{./ChapterE/digit3-all.png} 
\end{center} \caption{The 30 handwritten digits from
the Digit 3 data in red (left) and the   digits 
matched using Generalized Procrustes Analysis in green, where translation, rotation and scale have been removed (right). This registration gives a more focussed view of the population variation.}\label{Digit3-full} \end{figure} 


\subsubsection*{Case Study: DNA Molecule Data} 
An example for the DNA Molecule data (introduced in Section \ref{subsecI-C-C:Procrustes}) is given in Figure \ref{DNAfig2}. Here there are $n=50$  DNA molecules  \index{DNA Molecule data} 
of $k=22$ atoms in $m=3$ dimensions. The 
unaligned molecules have been pictured in the left hand panel in red.  In the right hand panel one can see the Procrustes registered molecules in green, which are obtained using Generalized Procrustes Analysis. Here  
the unimportant rotations, translations and scalings have been removed in the right-hand plot. Alignment can be thought of as a form of noise reduction which results in distinctly visible landmarks in this case.

\begin{figure}[H] \begin{center} \includegraphics[height=6cm]{./ChapterE/DNA-TFC-all1.png} \includegraphics[height=6cm]{./ChapterE/dna-all50-withscale.png} \end{center} \caption{The $n=50$ unaligned molecules from the DNA Molecule data in red (left). The   molecules matched using Generalized Procrustes Analysis in green, where translation, rotation and scale have been removed (right).}\label{DNAfig2} \end{figure} 
\index{landmark|)}

\subsection{Covariance Matrix Distances\label{subsec:Covariance-matrix-distances}}

\index{data objects!covariance matrices|(}\index{distance!Procrustes|)}\index{distance!on covariance matrices|(}

Following \cite{Dryden2009CovMat} let us compare and discuss some possible choices of distances between covariance matrices $\boldsymbol{S}_{1},\;\boldsymbol{S}_{2}\in PSD_{k}$ (recall the set of symmetric and positive semi-definite  $k \times k$ matrices from Section \ref{secBA.2:OverviewObjFeat}).
The Euclidean (or Frobenius) distance  between two matrices is given by \begin{equation} \delta_E( \mathbf{S}_1, \mathbf{S}_2) = \| \mathbf{S}_1 -\mathbf{S}_2 \|_F = \sqrt{ {\rm trace}\{ (\mathbf{S}_1 -\mathbf{S}_2)^t(\mathbf{S}_1 - \mathbf{S}_2) \} } \; , \label{dE} \end{equation} where 
$\| \mathbf{X} \|_F$ is the Frobenius norm (\ref{eq:Frobeniusnorm}).  However, $PSD_{k}$ is more naturally considered as a curved manifold. 
This is because when linearly extrapolating beyond the data, for example when using (\ref{eqE.I:FrechetInterp}) with $\alpha<0$ or $\alpha>1$ one can leave  $PSD_{k}$. This undesirable feature of using a linear extrapolation in a Euclidean space is also seen in the warping example of Figure \ref{figH.L:ToyNoSphere}, where the 
extrapolated paths using Euclidean PCA are no longer strictly monotone, sometimes 
leaving the space of warps. 


Another drawback of the Euclidean metric noted by  \cite{ArsignyPennec2006LogEuclideanDTI,Arsignyetal07} is that this distance suffers from a swelling effect when averaging {\it diffusion tensors}, which are $3 \times 3$ covariance matrices that represent the local movement of water molecules in white matter fibers of the brain. The swelling effect occurs when the average of two tensors using the Fr\'echet mean (\ref{eq:ChpEEqBFrechetMean}) has a bigger volume than either of the two individual tensors. 
An example which displays interpolation paths between common diffusion tensors is given in Figure \ref{interp} 
using the  Fr\'echet interpolation (\ref{eqE.I:FrechetInterp}) for a selection of distances. The swelling effect can be clearly 
seen for the Euclidean distance in the first row. \index{interolation path} 

\begin{figure}[H] \begin{center} \includegraphics[width=12cm]{./ChapterE/interpolation-paths.png} \end{center} \caption{Four different interpolation paths between the two common
diffusion tensors ($3 \times 3$ covariance matrices) in the left-most and right-most columns. Each covariance matrix is represented by an ellipsoid at the 95th percentiles of a multivariate normal distribution with that covariance matrix.  The interpolation paths are obtained from (\ref{eqE.I:FrechetInterp}) using Euclidean $\delta_E$ (1st row),  log-Euclidean $\delta_L$ (2nd row), Cholesky $\delta_C$ (3rd row) and Procrustes $\delta_S$ (4th row) distances.  } \label{interp} \end{figure}

These drawbacks can be overcome by using any of several covariance matrix distances, which essentially results in the covariance matrix objects lying in a curved manifold. One such distance is based
on the matrix logarithm. We write the logarithm of a  positive definite (all positive eigenvalues) covariance matrix $\mathbf{S}$ as follows.  Let $\mathbf{S} = \mathbf{U} \mathbf{\Lambda} \mathbf{U}^t$ be the usual spectral (i.e. eigenvalue) decomposition with $\mathbf{U} \in O(d)$ (the {\it orthogonal group} which consists of the {\it orthogonal  matrices}, i.e. those $d\times d$ matrices which have orthonormal columns with determinant $\pm1$, which consists of both rotation and axis flip operators) and   $\mathbf{\Lambda}$ diagonal with strictly positive entries. See Section \ref{subsecP-A-B:SVD} for further discussion of eigenanalysis of 
matrices. Let $\log \mathbf{\Lambda}$ be a diagonal matrix with  logarithm of the diagonal elements of $\mathbf{\Lambda}$ on the diagonal. The {\it matrix logarithm} of $\mathbf{S}$ is  given by $\log \mathbf{S} = \mathbf{U} ( \log \mathbf{\Lambda} ) \mathbf{U}^t$ and likewise the {\it matrix exponential} of $\mathbf{S}$ is     $\exp \mathbf{S} = \mathbf{U} (\exp \mathbf{\Lambda}) \mathbf{U}^t$.  \cite{ArsignyPennec2006LogEuclideanDTI,Arsignyetal07} propose the use of the log-Euclidean  distance, where Euclidean distance between the logarithm of covariance matrices is used for  statistical analysis, i.e. \begin{equation} \delta_L(\mathbf{S}_1,\mathbf{S}_2) = 
\|  \log(\mathbf{S}_1) - \log(\mathbf{S}_2)  \|_F  \; . \label{dL} \end{equation}
Another logarithm-based distance uses a {\it Riemannian metric} in the space of  square symmetric positive definite matrices, also known as the affine-invariant metric: \begin{equation} \delta_R(\mathbf{S}_1,\mathbf{S}_2) = \| \log (\mathbf{S}_1^{-1/2}\mathbf{S}_2 \mathbf{S}_1^{-1/2}) \|_F . \label{dR} \end{equation} which has been explored by several authors including   \cite{Pennecetal06},  \cite{Moakher05},  \cite{Schwartzman06},  \cite{Lengletetal06} and  \cite{Fletjosh07}. A disadvantage of both the log-Euclidean and Riemannian metrics is that the covariance matrices must be strictly positive definite: non-full-rank matrices cannot be considered. 
Advantages are that the log-based distances do not suffer from the extrapolation problem - all extrapolations remain positive definite. Also, the interpolated covariance matrices are much less prone to swelling.  In the example in Figure \ref{interp} the second row shows the interpolation path using the log-Euclidean distance. Note that the volume of the tensor in the middle of the path is smaller
than when using the other metrics. The Riemannian metric interpolated path is visually 
very similar for this example. \index{distance!log-Euclidean} 

We have many other choices of distance for comparing covariance matrices. For the one dimensional case we can look at the distance between variances (Euclidean distance), log variances (log-Euclidean or Riemannian distances) or square root of variances. 
The square root scale is a reasonable choice that works well in many applications. But there is non-uniqueness in choosing a suitable root in higher dimensions. 
One choice is the symmetric square root decomposition where $\mathbf{S}_i = \mathbf{L}_i^2, i=1,2$ and $\mathbf{L}_i$ is symmetric, i.e. the symmetric  square root matrix of $\mathbf{S} = \mathbf{U} \mathbf{\Lambda} \mathbf{U}^t$ is $\mathbf{L} = \mathbf{U} \mathbf{\Lambda}^{1/2} \mathbf{U}^t$.  The square root distance is given by  \begin{equation} \delta_{1/2}(\mathbf{S}_1,\mathbf{S}_2) = \| \mathbf{L}_1 - \mathbf{L}_2 \|_F  \; . \label{dS} \end{equation}
The square root metric is less prone to swelling than the Euclidean metric, but unlike the logarithm-based metrics it can 
be applied to non-full-rank covariance matrices. \index{distance!square root Euclidean} 

Another distance (which requires full rank matrices) is based on the Cholesky decomposition   
\citep{Wangetal04},  where $\mathbf{S}_i = \mathbf{L}_i \mathbf{L}_i^t$ and $\mathbf{L}_i = {\rm chol}(\mathbf{S}_i)$ is lower triangular with positive diagonal entries. The Cholesky distance is given by  \begin{equation} \delta_C(\mathbf{S}_1,\mathbf{S}_2) = \| {\rm chol}(\mathbf{S}_1) - {\rm chol}(\mathbf{S}_2) \|_F  \; . \label{dC} \end{equation}
The Cholesky decomposition is commonly carried out for efficient numerical computation.

A final choice is the best root in terms of minimizing distance over rotation and reflection between two root matrices. Such
a transformation includes permutations of the rows/columns of the covariance matrix and sometimes it is sensible to not care about this ordering in applications. This motivates the choice of the Procrustes  size-and-shape metric \citep{Dryden2009CovMat} between two $k \times k$ covariance   matrices $\mathbf{S}_1$ and $\mathbf{S}_2$ which is defined as  \begin{equation} \delta_S(\mathbf{S}_1, \mathbf{S}_2) = \mathop{\min_{\mathbf{R} \in O(k)}}   \|  \mathbf{L}_1 - \mathbf{L}_2 \mathbf{R} \|_F , \label{dSS} \end{equation} where $\mathbf{L}_i$ is a decomposition of $\mathbf{S}_i$  such that $\mathbf{S}_i = \mathbf{L}_i \mathbf{L}_i^t$ or $\mathbf{S}_i = \mathbf{L}_i^t \mathbf{L}_i$ ($i=1,2$)  (the latter representation leads permutations of rows/columns). This metric is generally quite similar to the square root metric, and has similar properties. If a permutation of the rows and columns is appropriate then the 
Procrustes metric could be useful, as the permutation matrices are a subset of the orthogonal matrices. This property is particularly useful for comparing unlabeled networks as discussed in \ref{subsec:Labeled-versus-unlabeled}. 

The last two rows of Figure \ref{interp} show the interpolation paths for the Cholesky and Procrustes  size-and-shape metrics. The paths look rather different due to the extra rotation in the Procrustes  method. The metrics are both prone to some swelling in the interpolation, but much less than the Euclidean metric.

A broad family of power metrics was introduced by \cite{Dryden2009CovMat}. 
\cite{Drydenetal10power} selected from such metrics with 
Box-Cox transformations, and the square root transformation 
was appropriate in their diffusion tensor application.

A promising metric for interpolating between diffusion tensors has been given by  \cite{JungSchwartzGrossier2015SPD} and  \cite{GrossierJungSchwarz2017ScalRotSPDs}, and involves decomposition of the tensors  using rotation and scaling. An advantage with their approach is that there is no swelling effect, and the interpolation paths are visually appealing. 

A further example with covariance matrices as data objects is given in Section \ref{secG-B:CovMatDataObjects} in the context of diffusion weighted MR imaging.

\index{data objects!covariance matrices|)}\index{metric|)}\index{distance!on covariance matrices|)}\index{interpolation path|)}\index{Frobenius norm|)}
