
\chapter{Manifold Data Analysis\label{chapG:Manifold-Data-Analysis}}

As discussed in Section \ref{sec:Shapes-as-Data}, at several points
in Chapter \ref{chap1B:Breadth-of-OODA} and in Section \ref{subsec:Covariance-matrix-distances},
data objects that naturally lie on a curved manifold\index{manifold|(}
are an important and challenging part of OODA. The value of analyzing
data in this way is first illustrated using the perhaps most straightforward
example of data objects on the unit circle in Section \ref{secG.E:DirectionalData}.
A brief introduction to aspects of manifold geometry needed for OODA
appears in Section \ref{secG.F:IntroManifolds}. Section \ref{secG-A:ShapeDataObjects}
highlights a succession of improvements in the area of PCA-like analyses
of data objects lying on manifolds and Section \ref{sec:Landmark-Shapes}
considers the shape space of landmark\index{landmark|(}s in more
detail. Section \ref{subsecG-A-A:CLTon Manifolds} briefly explores
the Central Limit Theorem for probability distributions on manifolds.
A fundamental concept stemming from these developments is Backwards
PCA, as discussed in Section \ref{secG-C:BackwardsPCA}. Section \ref{secG-B:CovMatDataObjects}
contains a brief overview of covariance matrices as data objects,
another area where data analysis on manifolds is an important concept. 

\section{Directional Data \label{secG.E:DirectionalData}}

Perhaps the oldest and most deeply investigated area of statistics
involving non-Euclidean data objects is called \emph{directional data\index{directional data|(}}.
There the data objects are \emph{angles} such as wind or magnetic
field directions. The large body of methods developed for analysis
of such data are summarized in the monographs by \citet{Fisher1995Circular,Fisher1987Spherical},
\citet{MardiaJupp2000Directional} and \citet{Jammalamadaka2001CircularStats}.
Good recent overview of the area can be found in \citet{PewseyGarcia2021DirectionalOverview}.

As illustrated in Figure \ref{figG.A:ToyDirectionalData} it can be
very useful to carefully take the non-Euclidean nature of the space
of angles into account. Both panels of Figure \ref{figG.A:ToyDirectionalData}
show toy data sets with angle data objects as green dots on the unit
circle, $S^{1}$ (recall this notation from (\ref{eqB.C:Spheres})).
In the left panel of Figure \ref{figG.A:ToyDirectionalData}, the
shown angles are $8^{\circ}$, $14^{\circ}$, $342^{\circ}$ and $350^{\circ}$.
One approach to finding the center (a notion of \emph{mean}) of the
data is to simply average the four numbers, and the result of that
is shown as the red plus sign. Displaying this average on the unit
circle $S^{1}$ makes it clear this is a very poor notion of \emph{center}
of the set of green dots. A much better notion of center is shown
as the blue plus sign, which is the Fréchet mean\index{Fréchet mean|(}
with respect to the arc length metric\index{metric!arc length|(}
(also called geodesic distance\index{!geodesic}), as defined in (\ref{eq:ChpEEqBFrechetMean}).
This example shows how taking the unit circle (perhaps the simplest
non-trivial example of a curved manifold) structure into account is
essential to doing a reasonable statistical analysis of directional
data. In general, it is similarly useful to carefully consider curvature
of the space in the analysis of data objects lying in arbitrary manifolds.

\begin{figure}[H]
\includegraphics[scale=0.7]{ChapterG/OODAbookChpGFigAToyDirectionalDataCMYK}

\caption{Toy examples of directional data, shown as green dots in $S^{1}$
both panels. Left panel shows the conventional mean of the angles
as numbers (red) is a poor notion of center of data, while the Fréchet
mean (blue) is more sensible. Right panel compares the extrinsic resultant
mean\index{resultant mean|(} (red) with the intrinsic\index{intrinsic|(}
Fréchet mean (blue).\label{figG.A:ToyDirectionalData}}
\end{figure}
Another example of data \index{centering}centering issues can be
seen in the right panel of Figure \ref{figG.A:ToyDirectionalData}.
That data set (shown as green dots) consists of 7 data angles slightly
less than $90^{\circ}$ (some shown outside the circle to indicate
replication) together with 6 data angles slightly more than $270^{\circ}$
on $S^{1}$. A notion of mean called an \emph{extrinsic} mean is shown
as the red x sign. This is computed by first taking the mean of the
green dots (actually the angles on the unit circle) in $\mathbb{R}^{2}$,
shown as the red circle. Next the red circle is projected out to the
unit circle to get the estimated extrinsic mean angle (red x sign).
In directional data, this notion of mean is often called the \emph{resultant
mean} because it is the result of treating each angle as a vector
in $\mathbb{R}^{2}$, finding the vector sum, and then projecting
that back to the circle $S^{1}$. Advantages of this extrinsic mean,
relative to the Fréchet mean (called \emph{intrinsic} because all
of its operations happen within $S^{1}$ rather than in the embedding
space $\mathbb{R}^{2}$) include straightforward fast computation
and better uniqueness properties of the projection\index{projection|(}.
Recall from Section \ref{secE.1:Fr=0000E9chet-mean} that the Fréchet
mean is not unique (see Figure \ref{figG.C:BadTangentPlane} for an
example). An extreme example is an equally spaced data set on $S^{1}$
where the Fréchet mean set\index{mean!set|(} consists of all the
midpoints between consecutive pairs of data points. \index{intrinsic}
\index{extrinsic}

A perhaps unattractive feature of the resultant mean is that it does
not correspond to one's guess as to where a notion of mean should
be, and that it can be rather unstable with respect to small changes
in the data. The Fréchet mean with respect to the arc length metric\index{metric!arc length|)}
(thus an intrinsic notion of center) is shown using the blue x sign.
Many would view this as a more intuitive notion of center in this
case. It is also more stable with respect to changes in the data which
remain in the 1st and 4th quadrants. This comparison is not unlike
comparing the standard Euclidean $\mathbb{R}^{1}$ mean and median
for slightly unbalanced, strongly bimodal data. One more point worth
noticing (that does not happen in $\mathbb{R}^{1}$) is that if the
data are shifted slightly to the 2nd and 3rd quadrants, the Fréchet
mean is even less stable in the sense of jumping to the complete opposite
side of the circle (while the resultant mean\index{resultant mean|)}
moves relatively less far). 

This directional data space of angle data objects can also be viewed
as another example having the equivalence class\index{equivalence class|(}
(i.e. quotient) structure as the triangle shape space indicated by
Figure \ref{fig:ToyTriangles}. In particular, using polar coordinates
every point (except the origin) in $\mathbb{R}^{2}$ can be represented
by an angle and a distance \index{distance|(} from the origin. Identification
of points with the same angle is an equivalence relation, whose resulting
equivalence classes are rays from the origin, which are characterized
simply by the angles. This is the perhaps simplest example of how
spaces of equivalence classes tend to naturally be represented by
curved spaces.

\section{Introduction to Shape Manifolds \label{secG.F:IntroManifolds}}

As discussed above, a convenient mathematical representation of shapes
is based on equivalence classes as data objects. In the case of angles
as data objects in Section \ref{secG.E:DirectionalData}, the resulting
object space\index{object space} yields the best analyses when viewed
as a \emph{curved manifold}. 

The general version of this process starts with a group of transformations.
For the case of directional data, this group is the set of constant
multiples (i.e. scalings) in $\mathbb{R}^{2}$. For the triangles
(or more generally for shape configurations represented by landmarks)
depending on the particular application it can be a set with any of
translations, rotations, scalings or perhaps permutation of the vertex
labels. Identification of elements which can be reached by a transformation
in the group yields an equivalence relation. For directional data,
points lying on each ray in $\mathbb{R}^{2}$ are thus all equivalent
to each other. In the case of triangle shapes, all translations, rotations
and scalings of a given triangle are identified with each other. The
space of resulting equivalence classes (also called \emph{orbit\index{orbit}s},
which are typically the data objects of interest) is called the \emph{quotient
space\index{quotient space}}. For example angles represent equivalence
classes in directional data, while sets of triangles are the shape
equivalence class data objects in the triangle case. Generally selection
of the group of transformations is a critical issue in OODA.

Quotient spaces are usually \emph{curved manifolds} (e.g. the unit
circle in the case of \index{directional data|)}directional data)
which are fundamental to differential geometry. These curved surfaces
lie in a higher dimensional Euclidean space and are \emph{smooth}
in the sense that at each point the surface is approximated by a tangent
(hyper)plane in the (limit operation) sense of shrinking neighborhoods.
In particular, a \emph{manifold} $M$ is a topological space which
locally resembles a Euclidean space near each point. This is illustrated
in the left panel of Figure \ref{figG.F:FletcherManifoldEg} which
appears as Figure 2.2 in the PhD dissertation \citet{Fletcher2004Dissertation}.
It shows how at the point $p$ on the manifold $M$, that smooth surface
is approximated by the tangent plane $T_{p}M$ on neighborhoods of
$p$.

Further a \emph{Riemannian \index{manifold!Riemannian}manifold} is
a smooth manifold which has a positive definite inner product defined
on the tangent space $T_{p}M$ at $p$. The distance between two points
on a Riemannian manifold is measured by the length of the shortest
path (minimal \emph{geodesic}) in the manifold connecting the points.
The Riemannian structure allows defining the useful concept of the
\emph{exponential map\index{exponential map|(}} at $p$ in Riemannian
geometry, also shown in the left panel of Figure \ref{secG-A:ShapeDataObjects}.
That is the inverse of the \emph{\index{log map}log map}, which (on
neighborhoods of $p$) maps points on the curved surface of $M$ to
points on the tangent plane in such a way that the geodesic  distance
(to $p$) along the surface of $M$ is the same as the corresponding
distance in the plane $T_{p}M$. The exponential-log terminology comes
from complex analysis as illustrated in the right panel of Figure
\ref{figG.F:FletcherManifoldEg}. In particular the first quadrant
of the complex unit circle is shown there, which contains the point
$e^{i\theta}=\cos\theta+i\sin\theta$, for a given angle $\theta$.
Recall that the length of the arc from the point $1+0i$ to $e^{i\theta}$
(shown as the thick dashed curve) is $\theta$ radians. That is also
the length of the tangent line segment from $1+0i$ to $1+\theta i$
(shown as the thick dashed line segment). In this sense the complex
exponential maps the tangent line segment to the manifold which is
the unit circle so that distances (to $1$) on the tangent line are
mapped to distances along the manifold (unit circle). Hence the manifold
exponential map generalizes that of complex analysis.

\begin{figure}[H]
\includegraphics[scale=0.4]{ChapterG/OODAbookChpGFigAFletcherManifoldEg}\qquad{}\includegraphics[scale=0.33]{ChapterG/OODAbookChpGFigFLogExp}

\caption{Left panel illustrates the concept of a curved manifold, and its tangent
plane approximation at each point, with the exponential map. Right
panel shows how the exponential-log terminology comes from complex
analysis. \label{figG.F:FletcherManifoldEg}}
\end{figure}
The PGA (recall Principal Geodesic Analysis\index{Principal Geodesic Analysis|(}
from Section \ref{sec:Shapes-as-Data}) modes of variation\index{modes of variation!manifold|(}\index{modes of variation!PGA|(}
illustrated in Figure \ref{fig:BladProsRect-PGA} can be viewed as
starting with a tangent plane centered at the Fréchet mean. The data
are mapped into that plane using the log map where they are analyzed
with conventional Euclidean PCA. Mapping the resulting PC direction
vectors back to the manifold using the corresponding exponential map\index{exponential map|)}
results in the corresponding principal geodesic modes of variation.
Reconstruction of corresponding s-rep figures gives the one dimensional
set of object space\index{object space} members that make up each
mode of variation as defined in Section \ref{subsecBA.1.1:FormalizingModes}.

For a deeper introduction to manifolds and Riemannian geometry, see
\citet{Chavel2006Geometry} and \citet{Lee2018Manifolds}. A far more
in depth discussion of statistical analysis on Riemannian manifolds
can be found in \citet{PennecSommerFletcher2019Riemannian}. 

An interesting variant of the data objects lying on a manifold theme
of this chapter is \emph{manifolds as \index{data objects!manifolds as}data
objects} as studied in \citet{KangReimherr2017manifold}. This idea
has been extended to data objects which are a union of a manifold
together with a function on that manifold (motivated by medical imaging
applications) in \citet{Lila2016SmoothPCA,Lila2018Statistics,Lila2019JASA}.

Again the focus of this chapter is the statistical analysis of populations
of data objects that naturally lie on curved manifolds. However, there
is another area where the term manifold intersects with statistics,
which is now called \emph{manifold learning}\index{manifold learning}.
The main idea is that in the analysis of high dimensional data, there
are situations where the data may lie near a relatively low dimensional
curved manifold, and the goal is to estimate that manifold from the
data. An informative example, from \citet{Tenenbaum2000ManifoldLearning},
is a set of digital photographs of the face of a statue from different
angles, with different light sources. The vectorize\index{vectorize}d
version of such a set of images\index{data objects!images} (as done
in Section \ref{sec:Images-as-Data}) is in some sense low dimensional,
since the 3-d structure of the face is the same in each, and the variation
is easily summarized by a few parameters. But this variation is nonlinear
in the space of vectorized images, so methods like PCA will not properly
find the low dimensional structure. The field of manifold learning
grew out of the \emph{\index{principal curves}principal curves} idea
of \citet{HastieStuetzle1989PrincipalCurves}, in the case of fitting
a one dimensional manifold, i.e. a curve, to data giving a notion
of first mode of variation. An interesting issue was how to find a
second mode of variation, which can be hard to imagine when trying
to generalize PCA, as there is no useful notion of orthogonality in
the space of curves. A solution to this problem was provided by \citet{LeblancTibshirani1994PrincSurf}
who proposed finding an appropriate surface, i.e. a higher dimensional
manifold, as an appropriate higher dimensional generalization of PCA.
A large amount of research in manifold learning, involving many variations,
was started by the implementations proposed by \citet{Tenenbaum2000ManifoldLearning}
and \citet{RoweisSaul2000ManifoldLearn}. See \citet{WangMarron2008EffDim}
for an early approach to estimation of \emph{effective dimensionality}\index{effective dimensionality},
i.e. the dimension of an approximating manifold in the presence of
noise using scale space\index{scale space} ideas. 

\section{Statistical Analysis of Shapes \label{secG-A:ShapeDataObjects}}

A\index{data objects!shapes|(} fundamental aspect of the statistical
analysis of shapes as data objects is that the quotient operation
mathematically removes irrelevant aspects of the variation from the
analysis. As in Section \ref{sec:Amplitude-and-Phase} such irrelevant
variation (e.g. translation, rotation and scaling in the case of triangle
shapes illustrated in Figure \ref{fig:ToyTriangles}) is usefully
called nuisance variation\index{nuisance variation}. In that section
on curve registration\index{curve registration!nuisance variation}
ideas, it was noted that depending on the context either phase\index{phase!data objects}
or amplitude\index{amplitude!data objects} data objects could be
of primary interest or else could merely be thought of as nuisance
variation. Less well known is that there is a parallel situation for
landmark representations. In particular, the study of plate tectonics
and continental drift is also based on landmark data, as studied in
\citet{Chang1988JASA} and \citet{RoyerChang1991}. However, in that
context an opposite choice of data objects was most sensible. In particular,
shape variation was the nuisance, which was removed by essentially
the inverse quotient operation. Thus the transformations were the
data objects in that analysis. 

As can be seen in other statistical areas, there are several rather
different ideas as to the ``proper'' way to statistically analyze
data objects lying on a manifold. An important dichotomy is between
extrinsic and intrinsic methods, using terminology introduced in Section
\ref{secG.E:DirectionalData}. The main idea of extrinsic analysis
is to treat the data objects as points in the embedding Euclidean
space ($\mathbb{R}^{3}$ in the case of the 2 dimensional manifold
shown in the left panel of Figure \ref{figG.F:FletcherManifoldEg}),
do the analysis there, and then project the result back to the manifold.
In contrast, intrinsic methods base the statistical analysis on operations
within the manifold. One way to achieve intrinsic methods is to work
with distance methods, such as the Fréchet mean and MDS as discussed
in Sections \ref{secE.1:Fr=0000E9chet-mean} and \ref{secE.2:Multi-dimensional-scaling},
with respect to a metric\index{metric|(} using distances along the
curved manifold surface (hence intrinsic).

This topic was introduced and motivated, using bladder-prostate-rectum
data\index{Bladder-Prostate-Rectum data} objects, in Section \ref{sec:Shapes-as-Data}.
As discussed there, several rather different representations of shape
can be found in the literature, including landmark and boundary representations.
Another approach to shape representation\index{shape representation}
takes the data objects to be transformations, typically of some template,
such as the \index{Large Deformation Diffeomorphic Metric Mapping}Large
Deformation Diffeomorphic Metric Mapping of \citet{DupuisGrenanderMiller1998LDDMM}
and \citet{BegMillerTrouveYounes2005LDDMM}. A particularly attractive
shape representation, based on quotienting over potential reparameterizations
to again give \index{equivalence class|)}equivalence classes (i.e.
orbit\index{orbit}s) as data objects can be found in \citet{SrivastavaKlassen2016book}.
That idea also has provided a compelling approach to curve registration
as discussed in Chapter \ref{chapH:FDA-Curve-Registration}. Statistical
analysis of any of these shape representations entails working on
a curved manifold at some level. 

In the rest of this section the focus is on skeletal representation\index{skeletal representation}s
(s-reps), as illustrated in Figure \ref{fig:BlaProsRect-Srep}, because
of their strong intuitive appeal and their ability to effectively
summarize important aspects of shape as discussed in \citet{PizerMarron2017ObjectStats}
and \citet{Pizer2020ShapeSkeletal}. These data objects naturally
lie in a curved manifold because the representation includes a number
of angles on $S^{2}$. While it is possible to attempt analysis of
s-rep data in the ambient space, e.g. to treat points on $S^{2}$
as lying in $\mathbb{R}^{3}$, a major problem is that analytic methods
such as PCA tend to leave the space where the data lie. As noted above,
one approach to this is extrinsic analysis. The idea is to do the
statistical analysis (e.g. the mean or PCA) in the ambient space and
then project back to the curved manifold. This approach works well
when the data lie in a small region of the manifold where there is
not much curvature, so there is little distortion caused by the extrinsic
approach.

When the data are distributed more broadly across the manifold, curvature
matters more, which provides strong motivation for intrinsic analysis
methods. One of these is the Principal Geodesic Analysis\index{Principal Geodesic Analysis|)}
(PGA) of \citet{Fletcher2004PGA}. As noted in Section \ref{sec:Shapes-as-Data},
the main idea of PGA is to think of the standard Euclidean PCA basis
as a set of orthogonal lines that (sequentially) best fit the data.
In PGA these best fitting lines are replaced by best fitting geodesics
(e.g. great circles on $S^{2}$) which are a natural analog of lines
along the surface of the manifold, that naturally determine modes
of variation. A deliberate choice that was made in PGA was to consider
only geodesics passing through the Fréchet mean. That restriction
allowed straightforward computation of PGA as discussed in Section
\ref{secG.F:IntroManifolds}, using the corresponding tangent plane.
Modes of variation generated by a PGA were shown in Figure \ref{fig:BladProsRect-PGA}. 

More recent research has led to major improvements in statistical
methodology that have been realized through improved incorporation
of the underlying geometry. \citet{Huckemann2010GPCA} made the important
observation that the effectiveness of the PGA tangent plane analysis
of \citet{Fletcher2004PGA} was strongly tied to the quality of the
Fréchet mean as a notion of center-point. Figure \ref{figG.C:BadTangentPlane}
shows how that can be a serious limitation. The data objects are shown
as blue dots distributed along the equator of the ordinary sphere
$S^{2}$. Note that since they follow a geodesic this data set is
essentially one dimensional. The Fréchet means of this data set are
shown as red crosses. The Fréchet mean is a mean set\index{mean!set|)}
(consisting of both north and south poles) in this case, because there
is not a unique minimum of the Fréchet criterion (\ref{eq:ChpEEqBFrechetMean}).
This is a consequence of least squares penalizing large distances
most strongly. In particular, a candidate point along the equator
will have some residuals extending half way around the sphere, while
the poles are only a quarter of the way from each data point. If the
data only approximately lie along the equator, the Fréchet mean will
typically be unique, but it will still be very far from all data points
(essentially choosing just one of the north and south poles), and
it will similarly be very poorly representative of the data.

Not only are the red Fréchet means in Figure \ref{figG.C:BadTangentPlane}
an unintuitive notion of data center, they are also exceptionally
poor candidates as a point of tangency for a PGA analysis. This is
illustrated by the yellow tangent plane in Figure \ref{figG.C:BadTangentPlane}.
The log map\index{log map|)} projection of the data objects onto
the equator, shown as green dots, follow the green circle which is
centered at the geodesic mean. This has the unattractive property
that it requires two PGA modes of variation to fully represent essential
aspects of the data. Such a two dimensional representation seems inefficient
because the data distributed on the equator of $S^{2}$ are one dimensional
(in the sense of following a one dimensional curve in the high dimensional
space). In fact these data objects are actually describable using
just the single (nonlinear) mode of variation (as defined in Section
\ref{Definition:Mode_of_Variation}) shown as the green circle. That
will generally result in much more efficient statistical analysis
(e.g. Bayes model fitting for segmentation\index{segmentation} as
seen in \citet{JeongStough2008,Jeong2009thesis}; \citet{PizerJeong2005,PizerMarron2017ObjectStats}).

\begin{figure}[H]
\includegraphics[scale=0.25]{ChapterG/OODAbookChpGFigCBadTangentPlane}

\caption{Toy example on sphere $S^{2}$ demonstrating poor performance by tangent
plane PCA. Data distributed along the equator shown as blue dots,
Fréchet means shown as red crosses, with log data shown as green.
Note log data lies along green circle in the tangent plane, thus requiring
two PGA modes of variation to fully capture the data relationships.
\label{figG.C:BadTangentPlane} }
\end{figure}
The solution proposed by \citet{Huckemann2010GPCA} was to consider
all geodesics in modeling the data, thus moving beyond the restriction
of PGA to only geodesics going through the Fréchet mean. An interesting
point here is that the first Euclidean PCA component can be viewed
as the line that best fits the data, which will necessarily contain
the sample mean, as easily seen using an ANalysis Of VAriance (ANOVA\index{ANOVA})
decomposition of sums of squares (a generalization of the \index{Pythagorean Theorem}Pythagorean
Theorem). However, in non-Euclidean situations (e.g. data objects
lying on curved manifolds) this is no longer true, so a conscious
decision needs to be made. In particular the use of PGA entails the
restriction to geodesics which go through the Fréchet mean. The \emph{\index{Geodesic PCA}Geodesic
PCA} proposal of \citet{Huckemann2010GPCA} considers all geodesics,
resulting in a first mode of variation along the equator for this
toy example, which gives a much more appropriate one dimensional representation
(i.e. mode of variation) of this data set. 

The toy example in Figure \ref{figG.C:BadTangentPlane} of data objects
lying on the equator of $S^{2}$ may appear to be artificial, but
related modes of variation are actually frequently important to medial
and skeletal shape representations,\index{shape representation} as
demonstrated in Figure \ref{figG.B:SimulatedSreps}. The left panel
shows the distribution of a single spoke over a number of realizations
from the bladder-prostate-rectum simulator model\index{Bladder-Prostate-Rectum data}
of \citet{Jeong2009thesis}, which captured realistic modes of variation.
Note that the data are quite broadly spread over the sphere, as in
the above discussed toy example.

Figure \ref{figG.B:SimulatedSreps} also demonstrates a limitation
of Geodesic PCA\index{modes of variation!PGA|)}. This is that the
spoke variation does not actually follow a great circle (i.e. a geodesic
such as the equator in Figure \ref{figG.C:BadTangentPlane}), but
instead follows a small circle (think Tropic of Capricorn on the earth).
Hence the Geodesic PCA fit still requires two modes of variation as
shown in the center panel. This motivated the \emph{Principal Arc
Analysis\index{Principal Arc Analysis}} proposed by \citet{JungFoskey2011},
which generalizes Geodesic PCA to allowing both small and great circle
fits to the data as modes of variation\index{modes of variation!Principal Arc Analysis}.
The benefit of this is shown in the right panel of Figure \ref{figG.B:SimulatedSreps},
where only one mode is needed to fully model this data set. This is
an example of where careful exploitation of the curved structure of
the manifold gives much more efficient analyses than are available
from distance methods as discussed in Chapter \ref{chapE:Distance-Based-Methods}.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterG/OODAbookChpGFigBSimdSreps}

\caption{Variation of a single spoke, from a bladder-prostate-rectum simulator
model. Left panel shows the blue data objects distributed on the sphere.
Center panel shows the results of a Geodesic PCA summary, requiring
two modes of variation. Right panel shows the corresponding Principal
Arc Analysis, which is much more statistically efficient since only
one mode of variation is needed.\label{figG.B:SimulatedSreps}}
\end{figure}
Motivated by landmark based shape analysis, the idea of using small
circles to give an analog of PCA was extended to higher dimensional
spheres, $S^{d}$ for $d\geq2$, by \citet{JungPNS2012}, to give
a method called \emph{Principal Nested\index{nested analyses} Spheres}
(PNS).\index{Principal Nested Spheres|(} In the special case of $S^{2}$,
PNS gives the result shown in the right panel of Figure \ref{figG.B:SimulatedSreps}.
One step of the iterative calculation of PNS is illustrated in Figure
\ref{figG.D:PNSintro}, where data objects distributed on $S^{d}$
are shown as blue dots. The gold circle represents the best fitting
sub-sphere of dimension $d-1$. The green dots represent the projections
of the data objects onto the gold sub-sphere. Note that the projection
(defined as the closest point in the sub-sphere) of each blue dot
is calculated with respect to the metric of arc length along the sphere.
These lengths are graphically depicted using cyan arcs. The gold sub-sphere
is taken to be ``best-fitting'' in the sense of minimizing the sum
of the squares of these cyan lengths. The green dot projections are
then the best co-dimension 1 (i.e. one lower dimension) approximations
of the blue data objects. Also fundamental to PNS is keeping the signed
(depending on which side of the gold sub-sphere) lengths of the cyan
arcs as the highest level PNS scores\index{scores!PNS}. This process
is then repeated iteratively down through dimensions, to generate
a full set of PNS scores. The final projections to $S^{1}$ play the
role of PNS 1 scores. One more aspect of this process worth noting
is that the Fréchet mean of the PNS 1 scores is a compelling notion
of center called the \emph{\index{PNS mean}PNS mean} or the \emph{backwards\index{backwards methods|)}
mean}\index{backwards mean|(}. In particular, for the Figure \ref{figG.C:BadTangentPlane}
toy example of data distributed around the equator of $S^{2}$, the
backwards mean is the quite reasonable one dimensional Fréchet mean
computed along the equator. 

\begin{figure}[H]
\includegraphics[scale=0.3]{ChapterG/OODAbookChpGFigDPNSintro}

\caption{Illustration of one step of Principal Nested Spheres. Raw data shown
as blue dots on $S^{d}$. Green dots are projections (along the surface)
onto the best fit version of $S^{d-1}$ shown as a gold circle, representing
the $d-1$ dimensional approximations. Level $d$ scores are signed
lengths of the cyan lines. \label{figG.D:PNSintro}}
\end{figure}
As noted in \citet{pPzerJung2013,PizerMarron2017ObjectStats,Pizer2020ShapeSkeletal},
PNS has provided major improvements in the Bayes segmentation\index{segmentation}
methods discussed in Section \ref{sec:Shapes-as-Data}. For example,
representations that required 20 PGA components, could be done using
just 13 PNS components, meaning very substantial noise reduction and
corresponding improved segmentation. This dramatic improvement appears
to be mostly driven by the natural prevalence of small circle \index{modes of variation!PNS}modes
of variation in s-rep data, as illustrated in Figure \ref{figG.B:SimulatedSreps}.
Strong performance of PNS in landmark shape contexts was also demonstrated
by \citet{JungPNS2012} and is further discussed in Section \ref{subsec:Principal-Nested-Shape}
for landmark shapes. Also PNS performed well in the study of rotational
deformations by \citet{Schulz2015RotDef} and in medical imaging by
\citet{KimJung2017PNSMedImag}. In addition PNS has been used in the
area of Fisher-Rao curve registration\index{curve registration!Fisher-Rao}
as discussed in Section \ref{secH.B:PNSonFisherRao}, which maps curves
to a high dimensional sphere. The benefits of this for blood glucose
monitoring data was demonstrated by \citet{YuLuMarron2017}. Another
important contribution of PNS is that it motivated the development
of the \emph{Backwards PCA\index{Backwards PCA|(}} viewpoint\emph{,}
discussed in Section \ref{secG-C:BackwardsPCA}.

While PNS made major improvements in the statistical analysis of s-rep
shapes, there is still room for improvement because s-reps do not
lie on a single high dimensional sphere. Instead they lie on a manifold
that is a product of spheres called a polysphere. A direct approach
to polysphere data, called \emph{\index{Composite PNS}Composite PNS}
in \citet{pPzerJung2013}, is to use PNS to obtain Euclidean representations
for each component sphere, which are then concatenated into long vectors
that are analyzed using PCA. While this approach has given useful
results as noted above, a weakness is that it makes the perhaps dubious
assumption of linear dependence across spheres. A more compelling
approach is the \index{Polysphere PCA}Polysphere PCA methodology
proposed by \citet{EltznerJH2015PolysherePCA,EltznerHM2018Torus}.
The key idea is to better employ PNS by first distorting the polysphere
space into a higher dimensional sphere, and then doing PNS there.
This method has been seen to give some improvements, but finding a
truly intrinsic approach to PCA for data objects lying on a polysphere
space remains an interesting open problem. An important special case,
that arises in computational chemistry, is torus type spaces, $\left(S^{1}\right)^{d}$,
where the data objects are vectors of angles.

\section{Landmark Shapes\label{sec:Landmark-Shapes}}

\index{data objects!landmark configurations|(}

An example of a quotient space of the type discussed in Section \ref{secG-A:ShapeDataObjects} is the shape space of $k > m$ landmarks in $m$ real dimensions, where $m$ is usually 2 or 3. Landmark configurations  
were introduced in Figure \ref{fig:ToyTriangles} and in Section \ref{subsecI-C-C:Procrustes} and an important method for
comparing the shapes of the configurations as data objects is Procrustes analysis, see 
Sections \ref{subsecI-C-C:Procrustes} and \ref{subsec:Generalized-Procrustes-analysis}. Examples of landmark shape data were given in those sections, in particular the Digit 3 data in $m=2$ dimensions  and the DNA Molecule data in $m=3$ dimensions.

Let $\mathbf{X}$ be a $k \times m$ matrix representing a shape object of $k$ landmarks in $m$ dimensions. The shape of the  configuration $\mathbf{X}$  is what is left when translation, rotation and isotropic scale are removed.  The geometry of the space is complicated, see \cite{Kendetal99}. However, practical progress can be made by first transforming to the {\it pre-shape}, which involves the easier steps of removing  translation and scale by centering and dividing through by centroid size:  
\begin{equation}
\mathbf{Z} =  \frac{  \mathbf{C}\mathbf{X}  } { \| \mathbf{C}\mathbf{X} \|_F }  \in S^{(k-1)m-1}  \label{eq:preshape}
\end{equation}
where $\mathbf{C}$ is the centering matrix (\ref{centeringmatrix}) and $\| \mathbf{C}\mathbf{X} \|_F$ 
is the centroid size (\ref{eq:centroid-size}).  The pre-shape space is a sphere in $(k-1)m$ dimensions $S^{(k-1)m-1}$, as defined in (\ref{eqB.C:Spheres}). \index{pre-shape space} 

Rotations of $\mathbf{Z}$ on the pre-shape sphere constitute  an orbit or fiber of the pre-shape space $S^{(k-1)m-1}$, which is an equivalence class. Equivalence classes and orbits were discussed in the context of triangle shape objects in Section \ref{subsec1A.B.B:ShapeRepresentations} and also in Section \ref{subsecH.A.C:FisherRaoMath} for warping functional data objects. Fibers on the pre-shape sphere  correspond one to one with shapes in the shape space, and so we can think of  a fiber as an equivalence class representing the shape of a configuration.  The pre-shape sphere is partitioned into non-overlapping fibers by the rotation group $SO(m)$ and the fiber is the orbit of $\mathbf{Z}$ under the action of $SO(m)$.  
\index{equivalence class} \index{fiber}

The idea behind constructing a distance in the quotient space  is to minimize the great circle distance on the pre-shape sphere between points, one in each fiber.
This shape distance is the {\it Riemannian shape distance}, which is the geometrically natural distance inherited from the    projection of the fibers on the pre-shape sphere to points in the  shape space. 
The {Riemannian shape distance} $\rho$ is given by  \begin{equation} \rho(\mathbf{X}_1,\mathbf{X}_2) = 2 \; {\rm arcsin} \left(\frac{\delta_{PP}}{2} \right) , \; \; \; 0 \le \rho \le \frac{\pi}{2} ,  \label{rho} \end{equation} 
where $\delta_{PP}$ is the partial Procrustes distance of (\ref{dP}). The interpolation path 
(\ref{eqE.I:FrechetInterp}) using $\rho$ as the shape distance is the shortest geodesic path between $\mathbf{X}_1$ and $\mathbf{X}_2$, which has length 
$\rho( \mathbf{X}_1,\mathbf{X}_2)$.  D.G.~Kendall introduced this
Riemannian distance in 
his seminal paper \cite{Kendall1984Shape}. \index{shape space} \index{distance!Riemannian shape}
\index{Riemannian shape} 


The projection from the pre-shape sphere to shape space is {\it isometric}, because distances are preserved. The Riemannian distance is an {intrinsic distance} in the shape 
space, and recall that discussion of intrinsic versus extrinsic distance was discussed in 
Section \ref{secG.E:DirectionalData}.   The shape  space is a type of { quotient space}, where   the rotation has been quotiented out from the pre-shape sphere  using optimization. The quotient space optimization imposes nonlinear constraints on the configuration, and
hence a non-Euclidean distance is appropriate. \index{quotient space} \index{extrinsic} \index{intrinsic} 
\index{isometric} 

Recall that an alternative shape distance is the partial Procrustes distance (\ref{dP}) which is a type of extrinsic distance between two shapes, where the extrinsic distance
 is measured in an embedding of the manifold, see \cite[Chapter 4]{DrydenMardia2016}. 
The partial Procrustes distance is the shortest {\it chordal distance} measured in the embedding space between two fibers on the pre-shape sphere. 

An additional choice of shape distance is the {\it full Procrustes distance} \begin{equation} \delta_{FP} = {\rm sin}(\rho)  \label{dF} \end{equation} which is another extrinsic distance and
a variation of the Procrustes distance defined at (\ref{dP}). 
The full Procrustes distance involves optimizing the Euclidean distance 
 in an embedding of the pre-shape sphere over both scale and rotation. The term ``full'' Procrustes is used
as both rotation and scale are used to match pre-shapes, whereas ``partial'' involves 
just rotation. 

When shape changes are small all three distances (\ref{rho}), (\ref{dP}), (\ref{dF}) are similar. In Table \ref{Tab:Shapedist} we see distances between the pairs of Digit 3s from Figure \ref{Digit3} and the pairs of DNA molecules in Figure \ref{DNAfig1}. Clearly the shape distances are small for the DNA data, and so all three shape distances are similar. However, for the Digit 3 Data example the shape differences are larger, and thus differ more.


\begin{table}
\begin{center}
\begin{tabular}{|c|c||c|c|c|}
\hline
        & Figure & $\rho$ &   $\delta_{PP}$ &   $\delta_{FP}$\\ 
\hline
Digit 3s &   \ref{Digit3} & 0.5155 & 0.5098 & 0.4930\\
\hline
DNA molecules & \ref{DNAfig1} & 0.1464 & 0.1459 & 0.1463\\
\hline
\end{tabular}
\end{center}
\caption{ Shape distances between the pair of Digit 3s and the pair of DNA molecules. This shows that these 
distances tends to be similar, especially when the shapes are closer to each other.}\label{Tab:Shapedist}
\end{table}


Further discussion of geometrical properties of the landmark shape space is given by  \cite{Kendall1984Shape,Kendall89,Lekend93,Kendetal99,Kentmard01,Small1996,DrydenMardia2016}.  The geometry is particularly straightforward in $m=2$ dimensions, where complex arithmetic can be used. In this planar case the shape space is a homogeneous space called a complex projective space $\mathbb{C}P^{k-2}$ \citep{Kendall1984Shape,Kent94}. \index{complex projective space} 

For triangle shapes a toy example was given in Figure \ref{fig:ToyTriangles}. In the triangle 
case the shape space is the complex projective space $\mathbb{C}P^{1}$, which is equivalent to a 
sphere in three dimensions \citep[p.54]{Kendall83,DrydenMardia2016}. A plot of the triangle shape space is given 
in Figure \ref{fig:triangle-sphere}, and the equilateral triangles are at the North 
Pole (anti-clockwise labels) and the South Pole (clockwise labels). Some longitudes with isosceles triangles are indicated, and  
the equator consists of all flat triangles (where the three points are collinear). If the shape sphere has radius $1/2$ then the great circle  distance
on the sphere is the same as the intrinsic Riemannian shape distance of (\ref{dR}). \index{shape space!triangles} 

\begin{figure}[H] \begin{center} \includegraphics[width=8cm]{./ChapterG/Triangle-shape-space.png} \end{center}  \caption{Kendall's spherical shape space  for triangles in $m = 2$ dimensions. The shape coordinates  are the latitude $\theta$ (with zero at the North pole)   and the longitude $\phi$. 
The anti-clockwise labeled equilateral triangle is at the North 
Pole ($A$) with $\theta=0$, and the clockwise labeled equilateral triangle is at the South Pole ($B$) with $\theta=\pi$.
Longitudes with isosceles triangles are shown in black ($\phi = k\pi/3, k=0,\ldots,5$), and 
flat triangles at the equator are shown in red at $\theta=\pi/2$.}\label{fig:triangle-sphere} \end{figure}

For $m=2$ the shape space is a homogeneous space, which is a nice space to deal with as the \index{homogeneous space} 
structure of the local geometry is the same at each point in the space. 
However, for $m>2$ the shape space is not homogeneous in the sense that the local geometry changes depending on the location in the space. \cite{Lekend93} studied the Riemannian structure 
of shape spaces in detail, and in particular gave explicit formulae for the average section curvature in the space at each location. Also, for $m > 2$ there are singularities in the shape space  when 
the $k$ points lie in an $m-2$ dimensional subspace, see  \cite{Kendall89}, so we
assume that we are away from such singularities. Furthermore 
the tangent space to the shape space itself is also complicated \citep{Le91b}. 
As the shape space is complicated in general, it is convenient to carry out 
analysis on the pre-shape sphere $S^{(k-1)m-1}$ and then adapt methods so that the analysis does not
depend on the rotations of the configurations. 


\subsection{Shape Tangent Space}\label{sec:shapetangentspace}
The tangent space to the pre-shape sphere can be decomposed into two complementary sub-spaces: the {\it horizontal}  tangent space of dimension $(k-1)m-m(m-1)/2-1$ which does not depend on rotation, and the {\it vertical} tangent space  of dimension $m(m-1)/2$ which contains the rotation information. For shape analysis we work with coordinates in the horizontal  tangent space to the pre-shape sphere. \index{shape tangent space}

Constructing the Procrustes tangent space involves first choosing a pole $\mathbf{P}$, a $k \times m$ matrix, which  is assumed to be non-degenerate \citep{Kentmard01} which means that the eigenvalues of  $\mathbf{P}^t\mathbf{P}$ satisfy $\lambda_{m-1} > 0$ (i.e. the rank of $\mathbf{P}$ is at least $m-1$).  We carry out  Procrustes rotation of $\mathbf{X}$ onto $\mathbf{P}$ (where $\mathbf{P}$ and $\mathbf{X}$ have been centered)  to give $\mathbf{X}^P = \mathbf{X}\hat{\mathbf{\Gamma}}$,  where $\hat{\mathbf{\Gamma}}$ is the optimal rotation from Ordinary Procrustes Analysis in Section \ref{subsecI-C-C:Procrustes}.  

The partial Procrustes tangent coordinate matrix is then \citep{Kentmard01} 
\begin{equation}
\mathbf{V} = \mathbf{X} \hat {\mathbf{\Gamma}} - \alpha \mathbf{P} , \label{tangentshape}
\end{equation}
where $\alpha = \cos \rho(\mathbf{X},\mathbf{P}) > 0$, with $\rho(\mathbf{X},\mathbf{P})$ the 
Riemannian distance between $\mathbf{X}$ and $\mathbf{P}$, and $\mathbf{X}$ is also non-degenerate.  We can write  
\begin{equation} 
\mathbf{X} =  (\alpha \mathbf{P} +  \mathbf{V}) \hat{\mathbf \Gamma}^t , \label{backproject} 
\end{equation}
where $0 \le \rho < \pi/2$, which is helpful for projecting back to the configuration space.




Multivariate statistical methods can then be carried out in the tangent space on the vectorized coordinates of $\mathbf{V}$, for example PCA, and $\mathbf{P}$ is taken as the Fr\'echet mean, obtained from Generalized Procrustes Analysis (see Section \ref{subsec:Generalized-Procrustes-analysis}).  We compute the vectorized Procrustes  tangent space coordinates, evaluate the sample covariance matrix $\mathbf{S}_V$ and carry out an eigendecomposition of $\mathbf{S}_V$. 
For detailed discussion of PCA and variants  see Chapter \ref{chapP:PCA-Details}.
The eigenvectors of $\mathbf{S}_V$ are the PC loadings, and the eigenvalues are the variances explained by each of the PCs. In order to visualize the modes of variation revealed by the PCs we back-project using (\ref{backproject}) to configurations that are located at multiples along each PC, for example 
at $\pm 3$ standard deviations along each PC. We consider an example of shape PCA in 
Section \ref{subsec:PCA-Digit3} for the Digit 3 Data. \index{visualization} \index{modes of variation!shape PCA}






\subsection{Case Study: Digit 3 Data}\label{subsec:PCA-Digit3}
In Section \ref{subsecI-C-C:Procrustes}  we introduced the Digit 3 data  \index{Digit 3 data}
and carried out Generalized Procrustes Analysis in Section \ref{subsec:Generalized-Procrustes-analysis}.  In Figure \ref{PCA-Digit3} we see a 
plot of the Fr\'echet mean obtained from Procrustes analysis in the middle column. 
We also carry out PCA of the partial 
Procrustes tangent coordinates of (\ref{tangentshape}). 
In the left and right columns modes of variation are displayed using configurations drawn at $\pm 3$ standard deviations 
along the 1st (top row), 2nd (middle row) and 3rd (bottom row)  principal component directions.  This matrix of scatterplots is essentially the transpose of that used to display modes of variation in Figure \ref{fig:BladProsRect-PGA}.  In this example the percentage of variability explained by these modes of variation are 50.4\%, 15.4\%, 12.8\%, and so 
more than three quarters of the  shape variability is explained by the first three PCs. There are often complicated combinations of features present in each mode of variation which may not be simple to interpret. In this case the perhaps dominant aspects are:  the first mode reflects the length of the middle prong in the digit, the second mode reflects tall thin digits versus short fat digits, and the third mode reflects an asymmetry in slanting up at one side versus the other. 

\begin{figure}[H] \begin{center} 
\includegraphics[width=10cm]{./ChapterG/digit3-pca.png} 
\end{center} \caption{The first three PCs for the Digit 3 data. Each row represents a mode of variation  with the Fr\'echet mean in the middle and configurations at $\pm 3$ standard deviations along each PC on the left and right. We see that the first mode reflects the length of the middle prong in the digit, the second mode reflects tall thin digits versus short fat digits, and the third mode reflects an asymmetry in slanting up at one side versus the other.} \label{PCA-Digit3} \end{figure}

An alternative view of these same modes of variation \index{modes of variation!shape PCA} is shown in Figure \ref{PCA-Digit3B}. The piecewise colored line shows the mean (center column of Figure \ref{PCA-Digit3}), with colors used to indicate the ordering of the landmarks from $1$ to $k$. The variation in each mode is displayed by including the mean plus 3 standard deviations (right column of Figure \ref{PCA-Digit3}) as the black piecewise line, and the mean minus 3 standard deviations (left column of Figure \ref{PCA-Digit3}) as the gray piecewise. The line segments connecting the landmarks show the paths of each landmark through the mode of variation (recall from Section \ref{subsecBA.1.1:FormalizingModes}, that is a set of data objects that is one dimensional in some sense). 

\begin{figure}[H] \begin{center} 
\includegraphics[height=3.8cm]{./ChapterG/digit3-pca-fig1.png} 
\includegraphics[height=3.8cm]{./ChapterG/digit3-pca-fig2.png}
\includegraphics[height=3.8cm]{./ChapterG/digit3-pca-fig3.png}
\end{center} \caption{PCA for the Digit 3 data with the mean (shown as the colored piecewise line) and the first (left), second (center) and third (right) modes of variation.  Mean $\pm 3$ standard deviations are shown as the black and gray piecewise curves.} \label{PCA-Digit3B} \end{figure}


\subsection{Case Study: DNA Molecule Data}\label{subsec:PCA-DNA}
In Section \ref{subsecI-C-C:Procrustes} we introduced the DNA Molecule data and carried out Generalized Procrustes Analysis in Section \ref{subsec:Generalized-Procrustes-analysis}. We also now carry out PCA using the shape  tangent coordinates. \index{DNA Molecule data}
In Figure \ref{PCA-DNA} the colored piecewise lines are plots of the mean shape configuration obtained from Procrustes analysis using translation, rotation and scale, with colors indicating landmarks as in Figure \ref{PCA-Digit3B}. 
The larger spheres are located at the landmarks, with the perspective view making the closer spheres larger, and those further away smaller. In addition variation is reflected in a fashion similar to  Figure \ref{PCA-Digit3B}, with vectors drawn to configurations which are 3 standard deviations along the 1st (left), 2nd (middle) and 3rd (right)  modes of variation.  In this example the percentage of variability explained by the first three modes of variation  are 26.8\%, 19.7\%, 14.7\%, and so just over 60\% of the  shape variability is explained by the first three PCs. It is difficult to interpret the modes of variation here, but features such as bending and twisting of the molecule may be present, as well as more localized features. The first mode reflects change at the opposite ends of each strand (at red and cyan atoms), 
and the second mode has orthogonal changes at the same points. 

In our examples up to this point, the focus has been on shape, where all of translation, rotation and scale (size) have been removed by the quotient operation. Yet in some situations size matters. For example the part of biology called {\it allometry} studies relationships between size and shape. An approach to this is a scatterplot matrix comparing data object size, quantified as the centroid size defined in (\ref{eq:centroid-size}), with the first 3 shape PC scores shown in Figure \ref{PCA-DNA-scores}. It can be seen that each PC score has some (strongest for PC1) correlation with centroid size.  
Also the remaining (shape) scores are unstructured which is consistent with lack of correlation between cases caused by the data thinning.  


\begin{figure}[H] \begin{center} 
\includegraphics[height=4cm]{./ChapterG/dna-pca-fig1.png} 
\includegraphics[height=4cm]{./ChapterG/dna-pca-fig2.png}
\includegraphics[height=4cm]{./ChapterG/dna-pca-fig3.png}
\end{center} \caption{PCA for the DNA Molecule data with the mean and first three modes of variation. Most modes seem to be various contrasts of variation at the red end versus the cyan fold point.} \label{PCA-DNA} \end{figure}

\begin{figure}[H] \begin{center} 
\includegraphics[width=10cm]{./ChapterG/dna-withscaling.pdf} 
\end{center} \caption{A scatterplot matrix plot of the centroid size and first three PC scores for the DNA Molecule data. We see that these PC scores are all correlated with centroid size.} \label{PCA-DNA-scores} \end{figure}

A more extensive analysis of the mean and covariance structure for the DNA Molecule data for six different types of DNA which are undamaged or damaged has been carried out by  \cite{Drydenetal17}. Other related work for investigating variability in time-series of DNA molecules includes {\it time-orthogonal principal components}  
\citep{Dryden2010CovMat}. 


\subsection{Principal Nested Shape Spaces\label{subsec:Principal-Nested-Shape}}

PCA in shape tangent space is appropriate when the data object variability
is not too large and the variability has linear structure. An alternative
is an extension of the Principal Nested Spheres (PNS) methodology
of Section \ref{secG-A:ShapeDataObjects} (\citet{JungPNS2012}),
adapted to Principal Nested Shape Spaces (PNSS) by \citet{dryden2019}.
The method of PNSS involves fitting a nested sequence of pre-shapes
to the Procrustes registered data. In order to speed up computation
an initial PCA can be carried out, e.g. by first reducing to 20 PCs
as in the following examples. Thus a small amount of noise is removed
by linear PCA, and then the main computational effort is spent in
fitting a sequence of pre-shape spheres of successively lower dimension.
As for PNS one has the option to fit either a great sphere (with radius
1) or a small sphere (with a smaller fitted radius) at each level.
The lowest levels of PNSS are the primary interest, and it is desirable
that these explain high percentages of shape variability. The very
final zero dimensional level (PNSS0) corresponds to the backwards
mean\index{backwards mean|)} shape (an alternative to the Procrustes
mean) and the first PNSS mode of variation (PNSS1) is a circular arc
in shape space with PNSS1 \index{scores!PNSS}scores indicating the
position along this arc. In order to visualize the effects of PNSS
we can back-project the PNSS0 to a configuration, and then draw arcs
to configurations along each PNSS mode of variation. 

In Figure \ref{PNS-DNA} we see the results from PNSS modes of variation
carried out on the DNA Molecule data of Section \ref{subsec:Generalized-Procrustes-analysis}
using a format similar to that of Figure \ref{PCA-DNA}. The back-projections
to configurations show the effects of each PNSS mode of variation.
Again the rainbow color is the PNSS mean. The black configuration
is +1 standard deviations along the first component, the gray configuration
is -1 standard deviations along each component for DNA, and the path
along each component is plotted as an arc. Note that the DNA data
has some very small radii spheres/circles fitted in the PNSS sequence,
indicating that they are rather far back. For PNSS the small circular
arcs for the DNA reflect the fact that the distance between the atoms
on either strand remains nearly fixed in the dataset. These plots
should be compared to the linear shape PCA plots that were given in
Figure \ref{PCA-DNA}. The nonlinearity of these nodes of variation
is reflected by the small arcs for the atoms being curved here, in
comparison to being straight lines due to the linear analysis there. 

The percentage variation in the first 5 PNSS scores and PCs can be
seen in the top rows of Table \ref{Tab:PCvar}. It is clear that the
nonlinear circular variation captured by PNSS1 explains a much higher
percentage of shape variation (69.8\%) than the linear PCA method
(26.8\%), because the mode of variation is allowed to follow a curved
path. See Section \ref{secH.B:PNSonFisherRao} for a parallel phenomenon
in the context of curve registration. 

  \begin{figure}[H] \begin{center} 
\includegraphics[height=4cm]{./ChapterG/dna-pns-fig1.png} 
\includegraphics[height=4cm]{./ChapterG/dna-pns-fig2.png}
\includegraphics[height=4cm]{./ChapterG/dna-pns-fig3.png}
\end{center} \caption{DNA Molecule data. PNSS with the PNSS mean and first three modes of variation at 1 standard deviations along the arcs.} \label{PNS-DNA} \end{figure}


\begin{table}
\begin{center}
\begin{tabular}{|c|c||c|c|c|c|c|}
\hline
  Data      & Figure & 1st & 2nd   &  3rd & 4th & 5th\\ 

\hline
DNA molecule & PCA & 26.8 & 19.7 & 14.2 & 7.2 & 4.6\\
              & PNSS & 69.8 & 7.3  & 4.1 & 3.0 & 2.0\\
\hline
Digit 3 &  PCA & 50.4 & 15.4 & 12.8 & 7.5 & 4.3\\
              & PNSS & 59.3 & 13.2  & 9.1 & 6.7 & 3.9\\
\hline
Digit 3 &   PCA & 43.6 & 18.4 & 13.2 & 9.0 & 4.9\\
w/o outlier & PNSS & 64.1 & 12.1  & 8.9 & 4.9 & 3.3\\
\hline

\end{tabular}
\end{center}
\caption{The percentages of variation explained by the first five modes of variation
using PCA and PNSS for the DNA Molecule and Digit 3 Data. This shows that the PNSS component explains higher percentages of shape variability than 
PCA in the first few components.}\label{Tab:PCvar}
\end{table}

We now apply the PNSS method \index{Principal Nested Shape Space} to the Digit 3 data of Section \ref{subsec:Generalized-Procrustes-analysis}. \index{modes of variation!PNSS} 
The percentages of variability compared to PCA are given in the middle row of Table \ref{Tab:PCvar}. 
In this example there is not such a major difference between the percentages of variability explained. 
The back projected PNSS mean and first three modes of variation are given in Figure \ref{PNS-Digit3}. 
These plots should be compared to the linear shape PCA plots that were given in Figure \ref{PCA-Digit3B}, 
and the structures of PNSS and PCA are quite similar here. The small arcs reflecting variation of each landmark are perhaps surprisingly roughly linear here, suggesting this nonlinear analysis is not far from that of Figure \ref{PCA-Digit3}.  \index{Digit 3 data} 
In this dataset there is a large outlier (the first observation) and the first PNSS and PC 
modes of variation are both in the direction
of the outlier. This outlier (and the paucity of data between it and the bulk of the data) appears to be the cause of the nearly linear modes of variation in Figure \ref{PNS-Digit3}.  \index{modes of variation!PNSS}

\begin{figure}[H] \begin{center} 
\includegraphics[height=3.8cm]{./ChapterG/digit3-pns-fig1.png} 
\includegraphics[height=3.8cm]{./ChapterG/digit3-pns-fig2.png}
\includegraphics[height=3.8cm]{./ChapterG/digit3-pns-fig3.png}
\end{center} \caption{PNSS for the Digit 3 Data with the mean (rainbow colored) and first three modes of variation 
(gray and black configurations) at $\pm 2$ standard deviation along the arcs.  } \label{PNS-Digit3} \end{figure}

To explore this, an alternative analysis without the outlier is given in Figure \ref{PNS-Digit3out} with the percentages of variability given in the last rows of Table \ref{Tab:PCvar}. The first PNSS mode of variation now represents more variability (64.1\%) 
than the first PC (43.6\%) and the nonlinear aspects of variability are now more clearly seen in the
back-projection for the first PNSS 
in Figure \ref{PNS-Digit3out}. In particular the first mode now once again uses a curved path to explain much more variation.  

\begin{figure}[H] \begin{center} 
\includegraphics[height=3.8cm]{./ChapterG/digit3a-pns-fig1.png} 
\includegraphics[height=3.8cm]{./ChapterG/digit3a-pns-fig2.png}
\includegraphics[height=3.8cm]{./ChapterG/digit3a-pns-fig3.png}
\end{center} \caption{PNSS for the Digit 3 data with the outlier removed with the mean (rainbow colored) and first three modes of variation (gray and black configurations) at $\pm 2$ standard deviation along the arc. PNSS1 is on the left. Here the PNSS1 variability is clearly nonlinear.} \label{PNS-Digit3out} \end{figure}




\subsection{Size-and-shape space}

In applications we have a choice of data objects, and in shape analysis one particular choice is 
whether to retain size in the analysis or whether to remove it. For example with the DNA Molecule data in Section \ref{subsec:Generalized-Procrustes-analysis} we removed location, rotation and scale, but an alternative is to choose to 
retain the original scale of each molecule and just remove translation and rotation. 
In this case the space of interest is the {\it size-and-shape space} (see \cite[Chapter 5]{DrydenMardia2016}). \index{size-and-shape space}

The size-and-shape space is again a quotient space, \index{quotient space}
but the normalization by size to obtain the pre-shape is not carried out. The 
Riemannian size-and-shape metric is based on the Frobenius norm of the residual matrix after carrying out 
Ordinary Procrustes Analysis using translation and rotation only, and so the distance
 between $\mathbf{X}$ and $\mathbf{P}$ is 
\begin{equation} 
\delta_{SS} = \| \mathbf{X}\mathbf{\hat \Gamma} - \mathbf P \|_F = \sqrt{ S(\mathbf{X})^2
+ S(\mathbf{P})^2  - 2 S(\mathbf{X})S(\mathbf{P}) \cos \rho(\mathbf{X},\mathbf{P}) } \; ,
\end{equation}
where $S$ is the centroid size of (\ref{eq:centroid-size}) and $\rho$ is the Riemannian distance of ({\ref{rho}). \index{distance!Riemannian} \index{distance!size-and-shape} 
The Fr\'echet mean size-and-shape is given by the sample mean of the configurations after 
carrying out Generalized Procrustes Analysis using translation and rotation only.
The size-and-shape tangent coordinates are given by the residual matrix
after Procrustes rotation, and specifically this is given by (\ref{tangentshape}) but with $\alpha = 1$. 
Again multivariate statistical methods can then be carried out 
on these tangent coordinates, and  
for the DNA Molecule Data the analysis is similar to the shape space analysis.


The data object choice of whether to scale or not is up to the data analyst. Figure \ref{PCA-DNA-scores} took an allometry approach by plotting pairwise views  of the first three shape PC scores for the DNA Molecule data (where size was removed),  together with the centroid size as a separate variable of interest. An alternative approach is to retain the original sizes of the objects throughout,  and in that case we work directly with size-and-shapes.  

\subsection{Further methodology}
Statistical shape analysis has undergone extensive methodological developments, with a large array of applications. Some recent examples include Riemannian smoothing splines applied to peptides \citep{Kimetal21}; regression models for human movement data \citep{Drydenetal20}; and measurement error 
modeling and Bayesian regression for face landmarks \citep{Duetal15,Dryden2019bayesian}.   Good overviews of contemporary developments in statistical shape analysis include  \cite{GrenanderMiller07,Younes10,Bhatbhat12,Hamelrycketal12,Bookstein14,DrydKent15,PatrangenaruEllingson2016,DrydenMardia2016,SrivastavaKlassen2016book,Turaga2016riemannian,Zheng2017statistical,PennecSommerFletcher2019Riemannian}.



\index{data objects!shapes}\index{data objects!landmark configurations|)}\index{Fréchet mean|)}\index{landmark|)}

\section{Central Limit Theory on Manifolds \label{subsecG-A-A:CLTon Manifolds}}

An important property of the mean $\overline{\boldsymbol{x}}$ of
a random sample\index{random sample} $\widetilde{\boldsymbol{x}}_{1},\cdots,\widetilde{\boldsymbol{x}}_{n}$
in Euclidean space $\mathbb{R}^{d}$ is the Central Limit Theorem\index{Central Limit Theorem}
(CLT). That Gaussian limiting distribution property is shared by many
other notions of center such as the median. As defined in Section
\ref{secE.1:Fr=0000E9chet-mean} a random sample means independently
drawn from a common probability distribution. The CLT states that
under rather broad conditions on the underlying probability distribution
(e.g. merely the existence of some moments) the limiting distribution
of $\overline{\boldsymbol{x}}$ as $n\rightarrow\infty$ is Gaussian
in the sense that 
\begin{equation}
n^{1/2}\left(\overline{\boldsymbol{x}}-\boldsymbol{\mu}\right)\stackrel{d}{\rightarrow}N\left(\boldsymbol{0}_{d,1},\boldsymbol{\Sigma}\right),\label{eqG-A:CLT}
\end{equation}
where $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ are the population
mean vector and covariance matrix, using the notation $\boldsymbol{0}_{d,1}$
for a vector of zeros from (\ref{eqB.G:0sMatrix}). This CLT is why
the Gaussian distribution plays such a central role in statistical
inference as taught in elementary courses, and also is why the Gaussian
distribution is often called ``the Normal distribution'' (hence
the symbol $N$ in (\ref{eqG-A:CLT})). For extension of the CLT to
manifolds, a fundamental aspect is that the CLT convergence in distribution
happens on neighborhoods that shrink at the rate $n^{-1/2}$. Connecting
this idea with the tangent plane definition of manifolds illustrated
in Figure \ref{figG.F:FletcherManifoldEg} has resulted in a number
of analogs as surveyed in Chapter 5 of \citet{PatrangenaruEllingson2016}.
A good classical reference for CLTs on manifolds is \citet{Terras1985HarmAnSymVol1},
which has been updated to \citet{Terras2013HarmAnSym,Terras2016HarmAnSym}.
Early work on large sample theory for intrinsic and extrinsic means
on manifolds was given by \citet{Bhatpatr03,Bhatpatr05}. An important
CLT for PNS (the Principal Nested Spheres method developed in Section
\ref{secG-A:ShapeDataObjects}) and related methods has been established
by \citet{HuckemannEltzner2018CLT4PNS}. An intriguing variation that
appears to happen only on positively curved manifolds is the case
of \emph{smeariness\index{smeariness}} as defined by \citet{EltznerH2019Smeariness},
where limiting Gaussians with a rate slower than $n^{-1/2}$ have
been discovered.

While the issue is not critical to analogs of the Central Limit Theorem,
because that happens on shrinking neighborhoods, it is natural to
contemplate analogs of the Gaussian distributions on manifolds, for
example this could be a reasonable approach to the development of
a meaningful notion of covariance matrix on manifolds. The Gaussian
analog issue can be approached in several ways. An approach that is
explained well in \citet{Terras1985HarmAnSymVol1} is via the heat
diffusion equation, whose Green's function is the Gaussian density
in Euclidean space. While this is natural for isotropic diffusion
(essentially spherically symmetric covariance) it is not clear how
to construct it in a non-isotropic way (i.e. to develop aspherical
covariances). Another approach is via normalized \emph{random walks}.
In Euclidean space $\mathbb{R}^{d}$, an $n$ step random walk that
is normalized by $n^{-1/2}$ converges (as $n\rightarrow\infty$)
to a Gaussian distribution. Extending this to a random walk on a manifold
is another approach to a Gaussian distribution. Handling anisotropy
after the first step is a challenge that has been elegantly addressed
using the moving frame idea of \citet{Sommer2019ProbGeomStat}.

When either of these approaches to an analog of the Gaussian distribution
is followed on the unit circle $S^{1}$, the result is the wrapped
normal distribution, 
\[
\sum_{j=-\infty}^{\infty}\frac{1}{\sqrt{2\pi}\sigma}e^{\frac{-\left(x-\mu-2\pi j\right)^{2}}{2\sigma^{2}}},
\]
see page 50 of \citet{MardiaJupp2000Directional}. While this is similar
to the normal distribution when the scale parameter $\sigma^{2}$
is relatively small (so all but one term of the sum are negligible),
this approach entails a substantial price in that the distribution
is not an exponential family because of the summation. Hence, likelihood
methods are challenging for numerical implementation, for statistical
inference, and for mathematical analyses. For these reasons, most
of the classical approaches to directional data\index{directional data}
are based on variations of the von Mises distribution\index{von Mises distribution},
which has the simple exponential family closed form
\[
\frac{1}{2\pi I_{0}\left(\kappa\right)}e^{\kappa\cos\left(x-\mu\right)},
\]
for centerpoint $\mu$ and scale parameter $\kappa$ (where $I_{0}$
is the modified Bessel function of order $0$), and thus straightforward
likelihood properties. A contrast with the Gaussian is that the von
Mises distribution does not seem to be the limiting distribution of
any analog of the CLT. However, it is the conditional (on norm $1$)
distribution of a bivariate isotropic Gaussian distribution. Furthermore,
for smaller values of its spread parameter (larger values of $\kappa$)
it is very close to the wrapped normal (and to the corresponding tangent
plane Gaussian). Many more related examples, asymptotics and discussion
can be found in a series of papers by \citet{Hotz2013,HotzHuckemann2015,HuckemannHotz2014}.

Much more discussion of extrinsic versus intrinsic\index{intrinsic|)}
issues, together with many other aspects of manifold data including
a wide array of theoretical results, are in \citet{PatrangenaruEllingson2016}.
That book goes beyond manifold data object space\index{object space}s
to \emph{stratified spaces}\index{stratified space|(}. Typical stratified
spaces are richer data spaces than manifolds that essentially consist
of a union of manifolds of differing dimension that are appropriately
attached together. An important example is the case of covariance
matrices as data objects, discussed in more detail in Section \ref{secG-B:CovMatDataObjects}.
As noted there the set of covariance matrices of the same rank is
usefully treated as a manifold. But a given set of data covariance
matrices may contain members of different ranks, which thus lie in
manifolds of differing dimensions. These manifolds are naturally connected
to each other through limiting operations where some eigenvalues tend
to 0. Stratified spaces also play an important role in the analysis
of tree\index{tree} structured data objects\index{data objects!tree structured}
considered in Chapter \ref{chapI:Tree-Structured-Data}. In parallel
to the above CLT discussion, even stranger versions of the CLT occur
on stratified spaces, such as the \emph{\index{stickiness}stickiness}
first observed by \citet{HotzEtal2013StickyMeans} (also mentioned
near the end of Section \ref{secE.1:Fr=0000E9chet-mean}), where the
rate of convergence can be much faster than the conventional Euclidean
CLT.

\section{Backwards PCA \label{secG-C:BackwardsPCA}}

PCA\index{Principal Component Analysis!backwards|(} is typically
introduced as a \emph{forwards\index{forwards methods|(}} decomposition,
where one starts with the mean (best 0-d least squares approximation),
then finds the 1-d best fitting subspace (best fit line) and continues
sequentially building up higher dimensional approximating subspaces.
But PCA can also be developed in a \emph{backwards} way, starting
with a hyperplane whose dimension is the rank of the data, and iteratively
finding a decreasing sequence of best fitting subspaces. The terminology
of forwards and backwards here is consistent with that used in classical
variable selection for linear models where forwards (start with small
models and sequentially add variables that explain the most variation
at each step) and backwards methods (start large and sequentially
delete variables explaining minimal variation) were popular before
the advent of sparse methods\index{sparse methods}. A detailed study
of the backwards PCA\index{backwards PCA|)} idea can be found in
\citet{DamonMarron2014}.

Little attention has previously been paid to this forwards versus
backwards issue in the context of PCA, perhaps because in the case
of Euclidean data they give the same result. As noted in Section \ref{subsecP-A-B:SVD},
both decompositions are simply calculated from either a Singular Value
Decomposition\index{Singular Value Decomposition|(} (SVD) of the
object mean centered data matrix, or equivalently from an eigen analysis\index{eigen analysis}
of the sample covariance matrix. The forwards (backwards) decomposition
of the data matrix consists of projection\index{projection|)}s onto
the increasing (decreasing respectively) sequence of affine spaces
generated by the eigenvectors. As noted in Section \ref{secG-A:ShapeDataObjects}
this equivalence is usefully thought of as a consequence of the \index{ANOVA}ANOVA
decomposition of sums of squares (Pythagorean Theorem\index{Pythagorean Theorem})
that underlies the eigen spectral analysis, which no longer holds
in non-Euclidean contexts. 

In such non-Euclidean cases, the Pythagorean Theorem no longer applies,
so backwards methods are no longer the same as forwards. \citet{DamonMarron2014}
review these ideas noting that the backwards approach seems to frequently
give more useful methodologies. An intuitive basis for this observation
was developed through viewing PCA in terms of a \emph{nested series
of constraints}. That idea can be understood through considering the
SVD of a data matrix, described in detail in Section \ref{subsecP-A-B:SVD}.
SVD is studied here instead of the usual mean centered PCA for notional
simplicity (otherwise subtraction of means obscures the main point).
Such an SVD results in a sequence of subspaces of $\mathbb{R}^{d}$
having rank $k=1,\cdots,r$ (where $r$ is the rank of the data matrix)
\[
S_{k}=\left\{ \boldsymbol{x}:\boldsymbol{x}=\sum_{j=1}^{k}c_{j}\boldsymbol{v}_{j},c_{j}\in\mathbb{R}\right\} ,
\]
where the $\boldsymbol{v}_{j}$ are the eigenvectors and the $c_{j}$
play the role of scores. Now to find the next smaller subspace $S_{k-1}$
one could simply change the index above the summation to $k-1$, but
much more insight comes from writing 
\[
S_{k-1}=\left\{ \boldsymbol{x}\in S_{k}:\boldsymbol{x}^{t}\boldsymbol{v}_{k}=0\right\} ,
\]
where $\boldsymbol{x}^{t}\boldsymbol{v}_{k}$ denotes the Euclidean
inner product. This highlights the fact that $S_{k-1}$ is easily
calculated from $S_{k}$ by incorporating a single linear constraint.
The general success of backwards methods is interpretable in terms
of it being relatively easy to sequentially find constraints, as done
by backwards methods. On the other hand, to compute a forwards analog
of PCA one needs to know the full sequence of constraints in advance,
and then sequentially relax them, which takes much more effort. 

Many variations of PCA have been proposed for a wide variety of purposes.
Some of those that focus on non-Euclidean contexts are discussed here,
while others are in Section \ref{secP-A:PCA-Views}. A common theme
in many non-Euclidean extensions is an unfortunate tendency to move
away from nested approximations. In particular, while one can usually
have a rank $k$ approximation for desired $k$, there is typically
no relationship between say the rank $k$ and $k-1$ approximations,
making multi-scale (in the sense of rank) interpretation quite challenging.
Nested versions of PCA analogs are also essential for developing modes
of variation that yield insightful visualization\index{visualization}s
of population structure, as detailed in Section \ref{secB2:Data-Visualization},
and also used above in many places. A major strength of backwards
approaches to non-Euclidean PCA is that they easily address both of
these concerns, using the fundamental idea of defining the method
through a nested series of constraints.

In very high dimensions, because backwards approaches are stepwise,
they can be computationally slow relative to conventional linear calculations
of PCA, which are done with fast Singular Value Decomposition\index{Singular Value Decomposition|)}s
(SVD) as discussed in Section \ref{subsecP-A-B:SVD}. This can be
mitigated by doing an initial SVD with a high enough rank to only
eliminate noise and retain almost all the signal followed by using
the backwards method to give good compression of the remaining signal
part of the data.

In addition to PNS,\index{Principal Nested Spheres|)} and PNSS, contexts
where backwards ideas have had (or apparently could have), a useful
impact include:
\begin{itemize}
\item \emph{Non-Negative Matrix Factorization}\index{Non-Negative Matrix Factorization}.
NMF was introduced in Section \ref{sec F.D:AlternativeDirns}. Recall
it is an approximation of the (uncentered) data matrix by a product
of loadings and scores\index{scores!NMF} matrices subject to the
constraint that all matrix entries are nonnegative. There is a substantial
literature on various ways to implement NMF (which reference the original
\citet{LeeSeung1999NNMF} paper). However most of them are non-nested
in the sense that one can request any rank of approximation, but the
results are completely different when a different rank is requested.
An exception is \emph{Nonnegative Nested Cone Analysis} (NNCA), proposed
by \citet{ZhangLuMarron2015NNCA}. NNCA specifically uses a backwards
implementation, starting with the nonnegative cone generated by the
data, and successively reducing the dimension through a sequence of
linear cone constraints, keeping that level of scores\index{scores!NNCA}
at each step.
\item \emph{Manifold Learning.} As discussed above, this approach is very
useful for high dimensional data sets that tend to lie quite close
to, or perhaps even on, a lower dimensional curved manifold. A drawback
to all current approaches is that they are not nested, so there is
no relationship between a rank $k$ and rank $k-1$ representation
of the data, and thus no decomposition into modes of variation. An
interesting open problem is to find a nested\index{nested analyses}
version of manifold learning\index{manifold learning}, where one
iteratively finds a rank $k-1$ approximation as a submanifold of
the one of rank $k$. This would make the most sense by minimizing
the sum of squared geodesic distances at each step. This would result
in a more natural analog of PCA scores\index{scores!manifold learning}
distribution plots\index{scores!distribution plot} than the current
PCA typically performed on unrolled versions of the $k$ dimensional
manifold.
\item \emph{Robust PCA.} Robust methods\index{robust methods} that are
generally useful in OODA are discussed in Chapter \ref{chapO:Robust-Methods}.
One overall approach is to use methods based on the $L^{1}$ norm.
Such methods typically draw robustness from behavior similar to the
median as discussed in Section \ref{secE.1:Fr=0000E9chet-mean}. Motivated
mostly by computational considerations, \citet{BrooksDulaBoone2013LIPCA}
have proposed a backwards approach to $L^{1}$ based PCA.
\item \emph{Principal Flows.\index{Principal Flows}} Another backwards
approach to PCA for data objects lying on a manifold is the Principal
Flow approach of \citet{PanaretodPhamYao2014PrincFlows}. This method
results in nongeodesic components with an interesting ability to locally
track maximal variation.
\end{itemize}
An important contribution to the discussion of backwards versus forwards
PCA is the \emph{\index{Barycentric Subspace Analysis}Barycentric
Subspace Analysis} idea of \citet{Pennec2016barycentric}. In addition
to providing an interesting general framework for manifold PCA\index{Principal Component Analysis!backwards|)},
that paper points out that both backwards and forwards\index{forwards methods|)}
methods are essentially greedy searches, which one expects can be
improved by solving an overall optimization, although that comes at
the price of increased complication. An interesting open problem is
how much gain is available from a globally optimized overall method,
relative to a simple backwards\index{backwards methods|)} approach.
\index{modes of variation!manifold|)} As noted in Section \ref{subsecG-A-A:CLTon Manifolds},
an approach to asymptotic analyses in this direction can be found
in \citet{HuckemannEltzner2018CLT4PNS}.

\section{Covariance Matrices as Data Objects \label{secG-B:CovMatDataObjects}}

Covariance matrices as data objects\index{data objects!covariance matrices|(}
appear in several applications. The one that has received the most
research effort has been \emph{diffusion tensor imaging}\index{diffusion tensor imaging|(},
started by \citet{Basser1994DiffTens}. In that area connectivity
within the brain is studied using 3-d images\index{data objects!images|(},
where at each voxel\index{voxel} the distribution of diffusion of
water molecules is summarized in terms of a covariance matrix, as
shown in Figure \ref{figG.E:DTIeg}. These 3-d covariance matrices
are represented by shaded ellipsoids representing a contour of the
corresponding Gaussian distribution in Figure \ref{figG.E:DTIeg},
after \citet{Dryden2009CovMat}. In fluid regions of the brain this
covariance is essentially a multiple of the identity matrix\index{identity matrix},
reflecting the ability of the molecules to diffuse freely in various
directions. In regions with fibers, such as neuronal brain connections,
one eigenvalue of the covariance matrix is much larger than the other
two, indicating only one degree of freedom of diffusion along the
direction of the fiber. Covariance matrices as data objects also play
an important part of the speech sounds as data objects discussed in
Section \ref{sec:Sounds-as-Data}. Additional examples appear in the
survey of \emph{Object Oriented Spatial Statistics\index{Object Oriented Spatial Statistics}}
by \citet{MenafoglioSecchi2017OOSpatialStats}.

\begin{figure}[H]
\includegraphics[scale=0.18]{ChapterG/real-dti-data}

\caption{Example of a coronal slice of diffusion tensors estimated from a human
brain. Tensors are represented by ellipsoids. Those stretched in a
single direction indicate fibers running through that location.\label{figG.E:DTIeg} }
\end{figure}
A number of approaches to data analysis of covariance matrices have
been proposed. Since matrices are just arrays of numbers they could
be vectorize\index{vectorize}d into long vectors, and then analyzed
with standard linear methods. But that approach generally leads to
severe complications because linear analyses such as PCA tend to leave
the space of positive semi-definite matrices. As noted in Chapter
\ref{chapE:Distance-Based-Methods} PCA distance methods provide a
useful approach to this challenge. An interesting aspect of covariance
matrices is that there are quite a few metrics of interest, with quite
different properties as seen in Section \ref{subsec:Covariance-matrix-distances}.
The log-Euclidean metric\index{metric!log-Euclidean} (\ref{dL})
was popularized by \citet{ArsignyPennec2006LogEuclideanDTI,Arsignyetal07}.
\citet{FletcherJoshi2004RiemannDTI} point out benefits of a Riemannian
metric\index{metric!Riemannian} (\ref{dR}) approach and \citet{Dryden2009CovMat}
advocate the Procrustes metric\index{metric!Procrustes} (\ref{dP}).
Properties of the metrics have been revealed by studying widely varying
interpolation path\index{interpolation path}s generated by different
metrics, for example in \citet{Dryden2009CovMat} and \citet{PigoliAstonDryden2014a},
and some examples were given in Figure \ref{interp}. 

A measure calculated from a covariance matrix that is very commonly used in diffusion tensor imaging is {\it Fractional Anistropy} \index{fractional anisotropy} 
$$FA = \left\{ \frac{k}{k-1} \sum_{i=1}^k (\lambda_i - \bar \lambda)^2 /  \sum_{i=1}^k \lambda_i^2 \right\}^{1/2} $$ where $ 0 \le FA \le 1$ and $\lambda_i$ are the eigenvalues of the diffusion tensor matrix (in which case $k=3$), or alternatively a covariance matrix.  Note that $FA \approx 1$ if $\lambda_1 >> \lambda_i \approx 0, i > 1$ (very strong principal axis of the confidence ellipse) and $FA = 0$ for isotropy. \cite{Basserpierpaoli96} developed the use of $FA$ in the context of diffusion tensor imaging. High values of $FA$ correspond to  white matter fibers which connect different parts of the brain. The study of the arrangement of these fibers is a key objective  of diffusion weighted MR imaging. 

An alternative to $FA$ for finding voxels indicating fibers is to use the full Procrustes shape distance $\delta_F$ (\ref{dF}) to isotropy.  This gives the {\it Procrustes Anisotropy} of  \cite{Dryden2009CovMat} \begin{equation} PA =  \left\{ \frac{k}{k-1} \sum_{i=1}^k (\sqrt \lambda_i - \overline { \sqrt \lambda } )^2 /  \sum_{i=1}^k \lambda_i \right\}^{1/2} , \end{equation} where $ \overline {\sqrt \lambda } = \frac{1}{k} \sum \sqrt{\lambda_i} $.  So $ 0 \le PA \le 1$, with $PA = 0$ indicating isotropy,  and $PA \approx 1$ indicating  a very strong principal axis.  \index{distance!Procrustes}

A final isotropy measure based on metrics $\delta_L$ or $\delta_R$ is the {\it Geodesic Anisotropy}  $$ GA  =  \left\{  \sum_{i=1}^k (\log \lambda_i - \overline {\log \lambda} )^2 \right\}^{1/2} ,  $$ where $0 \le GA < \infty$   
\citep{Arsignyetal07,Fillardetal07,Fletjosh07}, which has been used in  diffusion tensor analysis in medical imaging with $k = 3$. 

To compare the measures of anisotropy consider the diffusion tensors obtained from diffusion  weighted images in the brain that were plotted in Figure \ref{figG.E:DTIeg}. In Figure \ref{DTI} we see a slice from the brain with the $3 \times 3$ tensors displayed. This  image is a coronal view of the brain, and the corpus callosum and cingulum  can be seen.  
\begin{figure}[htbp] \begin{center} \includegraphics[width=12cm]{./ChapterG/DTI-3images.png} \end{center} \caption{The anisotropy measures $FA$ (left), $PA$ (middle), $GA$ (right) for the diffusion tensors, from Figure   \ref{figG.E:DTIeg}. This shows the differences in displays using the three measures, with arguably $PA$ having better contrast within the (relatively bright) corpus callosum.} \label{DTI} \end{figure}
At first sight all three measures appear broadly similar, but the $GA$ image has fewer brighter areas than $PA$ or $FA$.  \cite{Dryden2009CovMat} argue that $PA$ is slightly preferable here, offering  more contrast than the $FA$ image in the highly anisotropic region - the corpus callosum. The data were provided by Paul Morgan and for related discussion see \cite{Tench2002white}. 

Important\index{diffusion tensor imaging}\index{data objects!images|)}
work on the statistical analysis of covariance matrices and covariance
operators as data objects can be found in \citet{Dryden2009CovMat},
\citet{PigoliAstonDryden2014a,PigoliAstonDryden2014b} and \citet{AstonPigoliTavakoli2015}.
Analysis of diffusion tensor data, using local polynomial \index{smoothing}smoothing
methods can be found in \citet{YuanZhu2012DTI}, and a varying coefficient
model approach is given in \citet{YuanZhu2013DTI}.\index{distance|)}

The above works demonstrate that it has been very useful to understand
covariance matrices as data objects lying on a curved manifold. However,
an even more appropriate mathematical context is a stratified space\index{stratified space|)}.
As noted in Section \ref{subsecG-A-A:CLTon Manifolds} that is a connected
set of manifolds of different dimensions. Stratified spaces are appropriate
for covariance metric\index{metric!covariance}s of varying rank.
For each given rank $r$, the natural data space is a manifold whose
dimension is $r\left(r+1\right)/2$. These manifold\index{manifold|)}s
are naturally connected across rank through limiting operations where
eigenvalues tend to $0$.

While geodesics reveal interesting aspects of metric\index{metric|)}s,
they are also useful for interpolation between covariance matrices.
Also of interest for interpolation purposes are the non-geodesic paths
of \citet{JungSchwartzGrossier2015SPD} and \citet{GrossierJungSchwarz2017ScalRotSPDs}
which better correspond to intuitively expected interpolation for
diffusion tensors. 

Additional useful analytic methods for covariance matrices as data
objects, including hypothesis testing and even kriging can be found
in \citet{Secchi2013HDLSSinference}, \citet{PigoliAstonDryden2014b,Pigoli2016Kriging}
and \citet{Cabassi2017PermutationTest}.\index{data objects!covariance matrices|)}
