
\chapter{FDA Curve Registration\label{chapH:FDA-Curve-Registration}}

As\index{data objects!curves|(}\index{functional data analysis|(}
discussed in Sections \ref{sec:Amplitude-and-Phase} and \ref{secD3:Preprocess-Register},
amplitude\index{amplitude!variation|(} and \index{phase!variation|(}phase
variation can play an important role in FDA. As noted there and demonstrated
in Figure \ref{figH.H:Comparisons}, many curve registration\index{curve registration|(}
methods work well on the data sets for which they were designed, yet
typically require intensive manual tuning for application to a different
data set. OODA ideas, including the use of equivalence class\index{equivalence class}es
as data object representations (which were quite useful for quantifying
``shape'' variation in Sections \ref{sec:Shapes-as-Data}, \ref{secG.F:IntroManifolds}
and \ref{secG-A:ShapeDataObjects}), are seen in Section \ref{secH.A:Fisher-Rao-Reg}
to provide a mathematically rigorous and hence broadly applicable
general solution to this challenge. In particular, that section discusses
the Fisher-Rao methodology of \citet{Srivastava2011FisherRao}, which
was never published in the statistical literature apparently due to
the mathematical overhead involved. Section \ref{secH.B:PNSonFisherRao}
shows how the Fisher-Rao approach can be enhanced using Principal
Nested Spheres as introduced in Section \ref{secG-A:ShapeDataObjects}.

Figure \ref{figH.H:Comparisons} is a subset of Figure 1.2 of \citet{Srivastava2011FisherRao},
that was constructed by Wei Wu. In each row, the left hand column
shows a set of curve data objects which present different challenges
to peak alignment. The data in the top row is similar to the Shifted
Betas in Figure 10.2 below. The middle two rows consider other types
of shift variation. That in the bottom row comes from \citet{WuSrivastava2014}.
The middle three columns show analyses of each data set by a method
from the literature with publicly available software. In particular,
the second column uses the method of \citet{LiuMueller2004TimeWarp},
the third column is based on the self modeling warp implementation
of \citet{GerviniGasser2004Warping}, while the fourth column shows
the moment based approach of \citet{James2007CurveAlign}. These give
various degrees of alignment of the peaks, which are quite example
dependent. It is important to note that with sufficient manual tuning,
most of these can give good results. The last column shows the corresponding
analysis using the Fisher-Rao curve registration method described
in Section \ref{secH.A:Fisher-Rao-Reg}. For the reasons given in
that section, unlike the other methods shown in Figure \ref{figH.H:Comparisons},
simple default settings give high quality peak alignment in a fully
automatic way.

\begin{figure}[H]
\includegraphics[scale=0.65]{ChapterH/OODAbookChpHFigHComparisons}

\caption{Comparison of several curve registration methods (columns) on four
data sets (rows). Shows the Fisher-Rao approach (last column) provides
the most stable fully automatic alignment of peaks and valleys in
curves. \label{figH.H:Comparisons}}
\end{figure}

\section{Fisher-Rao Curve Registration\label{secH.A:Fisher-Rao-Reg}}

Further motivation for Fisher-Rao curve registration comes from the
Shifted Betas example in Section \ref{subsecH.A.A:ExampleShiftedBeta}.
The impact of misalignment on the sample mean was shown in Figure
\ref{figD-M:ToyShiftsnFRReg}. Alignment issues even more strongly
affect PCA\index{Principal Component Analysis|(} as seen in Figure
\ref{figH.A:ToyBadPCA} below. The basics of warping\index{warping|(}
functions that generally underlie curve registration are discussed
in Section \ref{subsecH.A.B:WarpIntro}. Details of the Fisher-Rao
approach appear in Section \ref{subsecH.A.C:FisherRaoMath}.

\subsection{Example: Shifted Betas Data \label{subsecH.A.A:ExampleShiftedBeta}}

The Shifted Betas\index{Shifted Betas data|(} raw data object curves
are shown in the top left panel of Figure \ref{figH.A:ToyBadPCA}.
These curves are $n=29$ Beta($\alpha,\beta$) probability densities,
with $\alpha=20-\beta=4.4,4.8,\cdots,15.6$. There is a very strong
phase mode of variation\index{modes of variation!amplitude and phase|(},
which is approximately tracked by the means of these distributions,
$EX=\frac{\alpha}{\alpha+\beta}$, and coded with a rainbow color
scheme. The name Shifted Betas is not entirely appropriate as these
densities all have support $[0,1]$. Nevertheless it does convey the
clear left to right shifting of probably mass. The colors clearly
indicate a single nonlinear mode of variation. The same colors are
used at many other points in this chapter. Because this single mode
is nonlinear it can admit insightful decomposition into other (sometimes
linear) modes. For example there is a clear amplitude (linear) mode
of variation in terms of heights. There is also a nonlinear mode of
variation in the widths since as probability densities each curve
has area 1 underneath so change in height entails change in width.
The presence of these modes of variation (in addition to the obvious
phase mode suggests that an intuitive analysis of this data set should
show 3 modes of variation. The sample object mean is shown in the
top center panel. For the reasons given in the discussion around Figure
\ref{figD-M:ToyShiftsnFRReg} this mean is not representative of the
data objects. In particular the object mean is only a very broad peak,
that is much lower than any individual, as highlighted using the same
vertical axis as the top left panel. The mean centered residuals are
shown in the top right panel, which are then decomposed into modes
of variation by PCA in the remaining rows.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigAToyBadPCA1}

\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigAToyBadPCA2}

\caption{Shifted Betas example, illustrating the impact of phase variation
on straightforward Euclidean FDA using PCA. Format similar to Figure
\ref{fig:10d-ToyPCA}, with data, object mean and mean residuals in
the top row. Modes of variation show PCA does a very poor job of displaying
the visually obvious shift in peaks as the dominant mode. This motivates
alternative approaches to understanding horizontal variation. \label{figH.A:ToyBadPCA}}
\end{figure}
The left plot in the second row shows the first PCA mode of variation.
This is disappointingly unlike either the phase shift or the peak
height modes discussed above. Instead it only shows that the magenta
and blue curves are higher on the left and lower on the right, while
the yellow and red curves exhibit opposite behavior. The residuals
in the center are harder to interpret. However, notice that use of
common vertical axes, as for the Tilted Parabolas\index{Tilted Parabolas data}
in Figure \ref{fig:10d-ToyPCA}, visually suggests that these residuals
represent a substantial share of the variation, in fact 39.7\% of
the sum of squares about the mean. The color scheme of PC1 scores\index{scores!harmonic pattern|(}
in the right panel at least suggests that this mode is mostly driven
by the shift in peaks, sweeping from left to right. 

The third row shows the PC2 mode of variation (28.4\% of the sum of
squares). The scores on the right seem to be a doubling of the frequency
of the sweep through the colors, going right to left and back. The
colors of the mode of variation on the left show only the later plotted
green to red curves, because these overplot the (hence invisible)
magenta to green curves. The 3 bump general pattern can be interpreted
as the second order correction of the crude 2 bump approximation of
PC1. The next such correction in approximation seems to be the 4 bump
pattern shown in PC3 (8.4\% of the variation) in the bottom row left.
These increasing bumps approximations are very reminiscent of the
behavior of orthogonal polynomial basis systems over increasing degree.
The corresponding scores in the bottom right again represent an increase
in frequency of back and forth sweeps. 

The main point of this example is that just a single mode of phase
variation is very poorly captured by PCA, with the signal power of
that mode being spread widely across the PC spectrum. Using terminology
from \citet{MarronRamsay2015}, PCA is a \emph{vertical} (amplitude)
analysis method, which only weakly provides insights about \emph{horizontal}
(phase) variation.

It is interesting to study the impact of the systematic back and forth
swings of scores using the \index{scatterplot matrix}scatterplot
matrix view (this device was used e.g. to visualize the Twin Arches
data\index{Twin Arches data} in Figure \ref{fig:50d-ToyPCAScatMat}).
This is done in Figure \ref{figH.B:ToyBadScores}, where the first
3 diagonal plots are the same as in the right column of Figure \ref{figH.A:ToyBadPCA}.
The sequential sweeping effects of the scores over different PC modes
generates some interesting patterns. Perhaps initially hardest to
interpret may be the PC2 vs. PC4 panel because the green to red points
again exactly overplot (and thus obscure) the magenta to green points.
Particularly attractive is the knot tied in the PC3 vs. PC4 panel.
Together these scores plots indicate that the scores follow a one
dimensional curve that bends around through the scores space. This
is a sense in which the data can be thought of as a single mode of
variation. That is consistent with the visual impression of the bundle
of curves in the upper left panel of Figure \ref{figH.A:ToyBadPCA},
and with the fact that those are essentially a one parameter family
of curves. This is another way to understand that while PCA can provide
good decompositions of vertical variation as seen at many points above,
it has its limits in the presence of horizontal variation, which motivates
the coming ideas.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigBToyBadScores}

\caption{Corresponding PC scores scatterplot matrix\index{scores!scatterplot}
for the Shifted Betas data of Figure \ref{figH.A:ToyBadPCA}. Shows
very systematic patterns that are common in applications of PCA to
data with a strong component of pure phase variation. \label{figH.B:ToyBadScores}}
\end{figure}
The second panel of the top row of Figure \ref{figH.B:ToyBadScores}
is an example of the horseshoe effect\index{horseshoe effect} studied
in the context of MDS in Section \ref{secE.2:Multi-dimensional-scaling}.
The explanation (from \citet{DiaconisHolmes2008Horseshoes}) given
there applies here, because curves that are farther from each other
have much less overlap, and thus are less far apart than might be
expected from linear extrapolation of the shorter distances. This
same issue lies at the heart of the shifting support example given
in \citet{Morton2017Horseshoe} in a micro-biome context. The distinctive
patterns of the \index{scores!harmonic pattern|)}scores in Figure
\ref{figH.B:ToyBadScores} seems to be a type of higher order horseshoe
effect, that could be seen in other examples when higher order components
are viewed. Another situation where such patterns have emerged is
the time varying chemical spectra studied in \citet{MarronWendel2004TimeSeries}.
These appear often enough to suggest the existence of a general mathematical
theory (perhaps something akin to Sturm-Liouville theory) behind this
phenomenon. That might be a generalization of the concept of \emph{harmonics\index{harmonics}\index{Functional Data Analysis!harmonics}}
in Fourier\index{Fourier!analysis} analysis.

The upper left panel of Figure \ref{figH.A:ToyBadPCA} suggests that
these generalized horseshoe effects are caused by the $L^{2}$ norm
doing a poor job of quantifying the shifting probability mass that
is a very important aspect of the clear mode of variation. Because
these are probability densities, this is another manifestation of
the issue discussed using a family of Gaussian densities in Section
\ref{subsecBA.2.1:ProbDistns}, that densities often provide a poor
representation of probability distributions as \index{data objects!probability distributions}data
objects. The approach to this that was taken there was to represent
the data objects using quantile\index{quantile} functions. In the
case of one dimensional distributions as data objects, that was seen
in Section \ref{subsecI-C-B:Wasserstein} to be equivalent to using
an optimal transport based MDS using the Wasserstein metric\index{metric!Wasserstein|(}.
Both approaches gave a less distorted representation of this major
mode of variation. A Wasserstein MDS appears in the left panel of
Figure \ref{figH.N:ToyBadWasser}. Note that common axes are used
to reflect the fact that there is much more variation in the MDS 1
scores\index{scores!Wasserstein MDS|(} than in MDS 2. In particular,
this Wasserstein decomposition of the variation seems to be much more
informative than the PCA shown in Figure \ref{figH.B:ToyBadScores}.
The colors suggest that MDS 1 is driven by the clear shifting of mass
in the upper left panel of Figure \ref{figH.A:ToyBadPCA}. The existence
of a smaller, but still visible second mode of variation suggested
by MDS 2, is consistent with the fact that the peak heights in the
density start high, decrease towards the middle, then increase again.
Unfortunately it is not simple to visually verify these ideas using
loadings\index{loadings} (as for the Spanish Mortality data\index{Spanish Mortality data}
in the left panels of Figures \ref{fig:Mortality-PC1} and \ref{fig:Mortality-PC2},
or for the Pan Cancer data\index{Pan Cancer data} in Figure \ref{figB-X:PanCanDWDLoadings})
because the MDS type of analysis only provides scores, not full modes
of variation. One approach is to plot individuals at key points of
the MDS scatterplot. For example the first MDS mode can be understood
from the magenta and red curves (recall corresponding colors), shown
in the top panel of Figure \ref{figH.A:ToyBadPCA}, as a shift mode
of variation. The second mode is a contrast between the light green
curve and those other two indicating variation in peak height and
width.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigNToyBadWasser}

\caption{Wasserstein MDS analysis of the Shifted Betas example. Left panel
shows \index{scores!Wasserstein MDS|)}scores giving much more intuitive
decomposition of the variation than the PCA shown in Figure \ref{figH.B:ToyBadScores}.
Scree plots (on $\log_{10}$ scale) in the right panel also shows
the MDS gives a much more efficient representation from another viewpoint,
\label{figH.N:ToyBadWasser}}
\end{figure}
The right panel of Figure \ref{figH.N:ToyBadWasser} gives a further
contrast of PCA with the Wasserstein MDS analysis of the Shifted Betas
data, using scree plot\index{scree plot}s (defined in Section \ref{subsecBA.A.C:ScreePlots})
based on the first 4 components. Because the numbers involved vary
over several orders of magnitude, the $\log_{10}$ of the proportion
of variation explained by each of the first 4 components is shown.
The green curve clearly shows how conventional Euclidean PCA spreads
the energy of the signal broadly across the spectrum, as discussed
near Figure \ref{figH.A:ToyBadPCA} above. On the other hand, the
magenta Wasserstein scree plot shows the great bulk of the energy
appears in MDS 1, with almost all the rest in MDS 2. The remaining
very small energy in the other two components seems to be at the level
expected from errors in the MDS optimization process. Thus MDS with
respect to the Wasserstein metric\index{metric!Wasserstein|)} gives
a far more insightful decomposition of the Shifted Betas data. Its
optimal transport motivation is clearly useful to keep in mind when
facing serious registration issues. This seems to have been done in
the microbiome literature where \citet{Fukuyama2012UnifracWasserstein}
discuss connections between the Wasserstein distance\index{distance!Wasserstein}
and the Unifrac metric\index{metric!Unifrac}. However, it still has
the drawback of not providing full modes of variation, because of
a general lack of loadings. The latter consideration provides motivation
for the Fisher-Rao methodology in the rest of this section. 

\subsection{Introduction to Warping Functions \label{subsecH.A.B:WarpIntro}}

Warping functions are a very useful concept for understanding approaches
to curve alignment and to amplitude - phase decompositions of the
type illustrated using the Bimodal Phase Shift data\index{Bimodal Phase Shift data}
in Figure \ref{fig:Sim-Warps}. For numerical convenience these are
typically defined on compact intervals. Allowing for linear transformation\index{transformation|(},
these intervals can be taken without loss of generality to be $[0,1].$
In this chapter, the data objects are taken to be functions $f:\left[0,1\right]\rightarrow\mathbb{R}$.
The function $\gamma:[0,1]\rightarrow[0,1]$ is a \emph{warp function\index{warp function|(}}
when it
\begin{itemize}
\item Is \emph{onto} {[}0,1{]} (i.e. surjective), 
\item Is strictly increasing, hence \emph{one to one} (injective). These
last two properties result in $\gamma$ being invertible (bijective),
\item Is a \emph{diffeomorphism}, meaning both $\gamma(x)$ and its inverse
are differentiable.
\end{itemize}
Let $\Gamma$ denote the set of all such warps. As noted in previous
sections, useful intuition comes from thinking of warps in terms of
local stretching and compression of the horizontal axis, as illustrated
in Figure \ref{figH.C:Warps}. The horizontal axis in each panel shows
an equally spaced grid on $[0,1]$. The function $\gamma(x)$ is shown
as a curve from $(0,0)$ to $(1,1)$. The thin lines indicate how
$\gamma(x)$ induces a stretching and compressing of the equally spaced
grid on the vertical axis. The upper left panel shows the diagonal
line, i.e. the \emph{identity warp\index{identity warp|(}}, $\gamma_{I}(x)=x$,
which leaves the spacing between grid points unchanged. The version
of $\gamma(x)$ in the upper right panel is convex, which results
in a compression of smaller values and stretching of larger values.
The lower left panel shows a concave $\gamma(x)$ which results in
the opposite type of stretching and compression. The $\gamma(x)$
in the lower right has both convex and concave parts, resulting in
compression at both ends, with stretching in the middle.

\begin{figure}[H]
\includegraphics[scale=0.3]{ChapterH/OODAbookChpHFigCWarps}

\caption{Illustrations of warping functions $\gamma(x):[0,1]\rightarrow[0,1]$,
showing how they are usefully thought of as stretchings and compressions
of the horizontal axis. \label{figH.C:Warps} }
\end{figure}

\subsection{Fisher-Rao Mathematics\label{subsecH.A.C:FisherRaoMath}}

There\index{data objects!amplitude - phase|(} were several major
contributions to the curve registration literature made by \citet{Srivastava2011FisherRao}.
The presentation here is aimed at maximal understanding of this important
set of ideas, which is done with various types of approximation. In
particular several function space technicalities are deliberately
ignored here. See that paper for the fully rigorous version, which
uses concepts such as generalized functions and related approximations
for full mathematical precision.

The first contribution of \citet{Srivastava2011FisherRao} was to
distinguish amplitude and phase variation using a quotient structure
parallel to that employed in shape analysis in Section \ref{secG-A:ShapeDataObjects}.
In particular, the set of all warps is a group of transformations
that plays a similar role as the translation - rotation - scaling
group of transformations in landmark based shape analysis. Following
that idea, two curves are \emph{identified} when one can be reached
from the other by a warping function, i.e. $f_{1}(x)\sim f_{2}(x)$
when there is a warp $\gamma(x)\in\upsilon$ so that $f_{1}(x)=f_{2}(\gamma(x))=f_{2}\circ\gamma(x)$.
That equivalence relation induces a quotient structure, i.e. a set
of equivalence class\index{equivalence class}es which are usefully
thought of as an insightful type of data objects, playing a role parallel
to the shape equivalence classes in Section \ref{secG-A:ShapeDataObjects}.
In this chapter, given a function $f:\left[0,1\right]\rightarrow\mathbb{R}$,
let its orbit\index{orbit|(} be $\left[f\right]=\left\{ f\circ\gamma:\gamma\in\Gamma\right\} $,
i.e. the set of all warps of $f(x).$ Any member of $\left[f\right]$
is called a \emph{representer\index{representer|(}} of that equivalence
class. Intuitive data analytic understanding of orbits will come from
careful choice of representers.

Additional insights into warp equivalence classes come from the examples
shown in Figure \ref{figH.D:WarpOrbits}. The left panel shows a set
of curves that are all in the orbit of a single curve $f(x)$ with
three bumps. Not only does each each curve share all 3 bumps, each
bump also has the same height. Furthermore the bumps are in the same
order, so they can reach each other through smooth compressions and
extensions in the spirit of the lower right panel of Figure \ref{figH.C:Warps}.
Thus they are all elements of $\left[f\right]$. A much different
example appears in the right panel of Figure \ref{figH.D:WarpOrbits}.
The curves there all have peaks of differing heights, so none of them
can be warped into any other, i.e. none of these are warp-equivalent.
Thus $\left[f_{i}\right]\neq\left[f_{j}\right]$ for $i\neq j$ in
this example.

\begin{figure}[H]
\includegraphics[scale=0.6]{ChapterH/OODAbookChpHFigDWarpOrbits}

\caption{Toy example of some elements of a single warp equivalence class (left
panel), and of bimodal curves none of which are warp-equivalent (right
panel). \label{figH.D:WarpOrbits}}
\end{figure}
The notion of warp-equivalence classes allows a formal, mathematically
rigorous, definition of the amplitude of a functional data object:
\begin{description}
\item [{Definition:}] The \emph{amplitude\index{amplitude|(}} of a function
$f$ is its warp equivalence class\index{equivalence class} $\left[f\right]$.\label{Definition:Amplitude}
\end{description}
The left panel of Figure \ref{figH.D:WarpOrbits} indicates that any
amplitude data object representer of $f$ contains all of the vertical
aspects of $f,$ such as the heights of all peaks, and depths of all
valleys, in addition to the relative relationship between peaks and
valleys. A parallel definition of phase data object is a deeper concept,
and requires development of some intermediate notions.

The second fundamental contribution of \citet{Srivastava2011FisherRao}
to the curve registration literature was the proposal of the Fisher-Rao
approach to quantifying the difference between amplitudes (i.e. orbits).
A fundamental aspect of that is \emph{warp invariance}\index{warp invariance|(}.
Such \index{metric|(}metrics on functions induce a proper metric
on the quotient space\index{quotient space|(} of orbits. The key
idea is illustrated in Figure \ref{figH.E:WarpDirection}, which is
Figure 3.3 of the PhD dissertation \citet{Lu2013Dissertation} that
was reproduced as Figure 2 in \citet{YuLuMarron2017}. This demonstrates
the inconsistency\index{inconsistency} of a common approach to the
estimation of warping functions, which is to find the warp $\gamma(x)\in\Gamma$
to minimize 
\begin{equation}
\left\Vert f_{1}(x)-f_{2}\circ\gamma(x)\right\Vert _{2}.\label{eqH.A:ClassicL2Warp}
\end{equation}
In Figure \ref{figH.E:WarpDirection} the functions $f_{1}$ and $f_{2}$
(called the Step Function toy example) are essentially step functions
as shown in the top left panel (actually these are just parallel coordinate
plot\index{parallel coordinate plot}s on a dense grid, so each has
one very steep jump between consecutive grid points). The warp function
$\gamma(x)$ that minimizes (\ref{eqH.A:ClassicL2Warp}) is shown
in the top center panel, as the blue (dashed) function that essentially
shifts the jump point in $f_{2}$ until it coincides with the jump
in $f_{1}$. Because $f_{1}$ is unwarped in the computation of (\ref{eqH.A:ClassicL2Warp}),
a red diagonal line, representing the identity warp $\gamma_{I}$
is also shown in the top center panel. The bottom center panel gives
insight into the minimum value of (\ref{eqH.A:ClassicL2Warp}), in
terms of the area of the rectangle shown with a light blue green color
(actually the $L^{2}$ norm is the height times the square root of
the length). The right column shows the result of minimizing the version
of (\ref{eqH.A:ClassicL2Warp}) where $f_{1}$ and $f_{2}$ switch
roles, so that the jump in the red $f_{1}$ is shifted over to the
jump in $f_{2}$, using the red warping function in the top right
panel (with the blue dashed diagonal line indicating $\gamma_{I}(x)$,
i.e. no warping for $f_{2}$ in this case). This results in a minimizing
value of (\ref{eqH.A:ClassicL2Warp}) represented using a pink rectangle
in the lower right panel. Notice this $L^{2}$ minimum value is much
different than in the center column showing that criteria of the form
(\ref{eqH.A:ClassicL2Warp}) cannot provide a metric (i.e. distance\index{distance|(}
as discussed in Chapter \ref{chapE:Distance-Based-Methods}) on the
quotient space of amplitudes. In particular, Figure \ref{figH.E:WarpDirection}
reveals that an attempt to base a quotient space metric\index{metric!quotient space}
on (\ref{eqH.A:ClassicL2Warp}) would fail due to a lack of symmetry,
i.e. $\delta\left(\left[f_{1}\right],\left[f_{2}\right]\right)\neq\delta\left(\left[f_{2}\right],\left[f_{1}\right]\right)$.
The lower left panel is explained in detail below.

\begin{figure}[H]
\includegraphics[scale=0.5]{ChapterH/OODAbookChpHFigEWarpDirection}

\caption{Step Function example illustrating failure of classical approaches
to provide a metric on the space of warp equivalence class\index{equivalence class}es.
Data objects in the upper left, with aligning warps in the other upper
panels. Lower right panels shows impact of warps on the $L^{2}$ norms
(which are quite different). Lower left shows the impact of the intermediate
Fisher-Rao warp estimates $\gamma_{1,K}$, $\gamma_{2,K}$. \label{figH.E:WarpDirection}}
\end{figure}
A natural approach to developing a metric on the space of amplitudes
(i.e. the quotient space) is to start with a metric $\delta$ on the
original space of functions that is \emph{warp invariant} in the sense
that 
\begin{equation}
\delta\left(f_{1},f_{2}\right)=\delta\left(f_{1}\circ\gamma,f_{2}\circ\gamma\right),\label{eqH.B:WarpInvariance}
\end{equation}
for all warps $\gamma\in\Gamma$. That property of a metric is usefully
understood as implying that the equivalence classes are all ``parallel''
with respect to $\gamma$ as illustrated in Figure \ref{figH.F:WarpInvariance}.
In particular, moving through both equivalence classes by the same
warp, leaves the distance unchanged. The property (\ref{eqH.B:WarpInvariance})
ensures that $\delta$ induces a well defined metric on the warp-equivalence
quotient space as 
\begin{equation}
\delta\left(\left[f_{1}\right],\left[f_{2}\right]\right)=\underset{\gamma\in\Gamma}{\min}\;\delta\left(f_{1}\circ\gamma,f_{2}\right).\label{eqH.D:WrpInvMetric}
\end{equation}
In particular the warp invariance (\ref{eqH.B:WarpInvariance}) entails
that $\delta\left(\left[f_{1}\right],\left[f_{2}\right]\right)$ does
not depend on the representers $f_{1}$ and $f_{2}$. 

\begin{figure}[H]
\includegraphics[scale=0.35]{ChapterH/OODAbookChpHFigFWarpInvariance}

\caption{Illustration of how a warp invariant metric results in parallel quotient
space \index{orbit|)}orbits. \label{figH.F:WarpInvariance}}
\end{figure}
The Fisher-Rao analysis of amplitude and phase variation is motivated
by the area called \emph{information geometry\index{information geometry}},
that was started by \citet{Rao1945Information}, and more fully developed
by \citet{Amari1985DiffGeomStat}, which studies parametric likelihood
estimation from a geometric perspective. The Fisher-Rao approach achieves
warp invariance by approaching differentiable functions $f$ through
their \emph{Square Root Velocity Function\index{Square Root Velocity Function|(}}
(SRVF) 
\begin{equation}
q_{f}(t)=sign\left(f'(t)\right)\sqrt{\left|f'(t)\right|},\label{eqH.C:SRVF}
\end{equation}
i.e. a signed version of the square root of the derivative. The SRVF
transformation can be viewed as representing the variation in $f$
in the sense that 
\[
f(x)=f(0)+\int_{0}^{x}q_{f}(t)\left|q_{f}(t)\right|dt,
\]
where the integrand is similarly thought of as a signed version of
the square of $q_{f}$. For a given warp $\gamma\in\Gamma$, integration
by substitution shows that the SRVF representation is very useful
for curve registration because it gives the $L^{2}$ norm a type of
warp invariance\index{warp invariance|)} in the sense that 
\[
\left\Vert q_{f_{1}\circ\gamma}-q_{f_{2}\circ\gamma}\right\Vert _{2}=\left\Vert q_{f_{1}}-q_{f_{2}}\right\Vert _{2}.
\]
 The Fisher-Rao metric\index{metric!Fisher-Rao} on the quotient space
is nearly the $L^{2}$ norm on SRVFs, but an alignment of $q$ functions
is still needed. In particular, inducing a metric on the space of
amplitudes as in (\ref{eqH.D:WrpInvMetric}) gives the \emph{amplitude
\index{distance!amplitude}distance} between functions
\begin{equation}
\delta_{A}\left(f_{1},f_{2}\right)=\min_{\gamma\in\Gamma}\left\Vert q_{f_{1}\circ\gamma}-q_{f_{2}}\right\Vert _{2}.\label{eqH.E:AmplitudeDistance}
\end{equation}
It is straightforward to show that $\delta_{A}$ induces a proper
metric on the space of amplitudes, and also that for $f_{1},f_{2}\in\left[f\right]$,
$\delta_{A}\left(f_{1},f_{2}\right)=0$. 

Given a set of functions $f_{1},\cdots,f_{n}$, Fisher-Rao analysis
starts with the \index{Fréchet mean|(}Fréchet mean (called the \emph{Karcher
mean\index{Karcher mean|(}} in \citet{Srivastava2011FisherRao},
despite the objection of \citet{Karcher2014NameObjection} discussed
in Section \ref{secE.1:Fr=0000E9chet-mean}) of the amplitudes $\left[f_{1}\right],\cdots,\left[f_{n}\right]$
with respect to $\delta_{A}$. In principle, using (\ref{eq:ChpEEqBFrechetMean})
this gives 
\[
\overline{\left[f\right]}=\arg\min_{\left[f\right]}\sum_{i=1}^{n}\delta_{A}\left(\left[f\right],\left[f_{i}\right]\right)^{2}.
\]
But $\overline{\left[f\right]}$ is only defined up to warps, leaving
open the question of the choice of best representer. For now, let
$\overline{f}_{C}$ denote a \emph{candidate} representer of $\overline{\left[f\right]}$,
and let its SRVF be $q_{\overline{f},C}$. For $i=1,\cdots,n$ let
\[
\gamma_{i,C}=\arg\min_{\gamma\in\Gamma}\left\Vert q_{f_{i}\circ\gamma}-q_{\overline{f},C}\right\Vert _{2}
\]
 denote the corresponding warps of the $f_{i}$ into $\overline{f}_{C}$. 

The distance $\delta_{A}$ could also form the basis of an MDS (from
Section \ref{secE.2:Multi-dimensional-scaling}) version of PC scores
on \index{amplitude|)}amplitudes. However, direct computation of
$\delta_{A}$ is quite challenging, so the approach detailed below
is numerically preferable, and also results in fully interpretable
modes of variation.

The third contribution of \citet{Srivastava2011FisherRao} is a parallel
definition of phase distance\index{distance!phase}, which integrates
appropriately with $\delta_{A}$ to give the Fisher-Rao approach to
curve alignment, which is based on also applying the SRVF idea to
the warp functions. In particular, given a differentiable warp $\gamma\in\Gamma$,
define its SRVF 
\[
\psi_{\gamma}(t)=\sqrt{\gamma'(t)}.
\]
The definition of $\psi_{\gamma}$ does not need the $sign$ function
used in (\ref{eqH.C:SRVF}) because $\gamma$ is strictly increasing.
A very important property of the SRVFs of warps is 
\begin{equation}
\left\Vert \psi_{\gamma}\right\Vert _{2}^{2}=\int_{0}^{1}\psi_{\gamma}(t)^{2}dt=\gamma(1)-\gamma(0)=1,\label{eqH.F:WarpsOnSphere}
\end{equation}
i.e. the SRVFs of warps $\psi_{\gamma}$ all lie on the unit sphere
in function space. A direct analog of the arc length distance on $S^{d}$
(which is the arc cosine of the inner product of the corresponding
vectors in $\mathbb{R}^{d+1}$) is the distance between warps 
\[
\delta_{P}\left(\gamma_{1},\gamma_{2}\right)=\cos^{-1}\left(\int_{0}^{1}\psi_{\gamma_{1}}(t)\psi_{\gamma_{2}}(t)dt\right).
\]
Given $f_{1},\cdots,f_{n}$, a candidate center $\overline{f}_{C}$
and corresponding candidate warps $\gamma_{1,C},\cdots,\gamma_{n,C}$,
result in a candidate notion of the \emph{\index{distance!phase}phase
distance} between functions
\[
\delta_{P,C}\left(f_{1},f_{2}\right)=\cos^{-1}\left(\int_{0}^{1}\psi_{\gamma_{1,C}}(t)\psi_{\gamma_{2,C}}(t)dt\right).
\]
Again an analog of PC scores could be derived by MDS applied to the
distance $\delta_{P,C}$, but a derivation of full modes of variation,
as well as sensible choice of the candidate $\overline{f}_{C}$ and
efficient computation are described below. 

A natural sense of average representer of $\overline{\left[f\right]}$
is the Karcher mean function $\overline{f}_{K}$, which is the element
of $\overline{\left[f\right]}$ with the property that the corresponding
warps $\gamma_{i,K}$ of each $f_{i}$ into $\overline{f}_{K}$, are
centered at the identity warp, $\gamma_{I}(x)$. Recall $\gamma_{I}(x)$
was shown as the red diagonal line in the middle panel (and the blue
dashed line in the right panel) of Figure \ref{figH.C:Warps}. In
particular, $\overline{f}_{K}$ is the member of $\overline{\left[f\right]}$
which gives warps whose SRVFs have Fréchet mean 
\[
\psi_{\gamma_{I}}=\arg\min_{\psi}\sum_{i=1}^{n}\left\Vert \psi-\psi_{\gamma_{i,K}}\right\Vert ^{2}.
\]
Choosing the candidate representer $\overline{f}_{C}$ to be $\overline{f}_{K}$
in the above calculations results in the final version of the phase
distance between functions $\delta_{P}$. 

The fourth major contribution of \citet{Srivastava2011FisherRao}
is an iterative algorithm for numerical approximation of these quantities,
based on dynamic programming ideas. This is a special case of the
general approach to registration problems described in \citet{SrivastavaKlassen2016book}.
Efficient software can be found in \citet{TuckerWuSrivastava2013},
with a corresponding R version on CRAN in the package $\mathtt{fdasrvf}$. 

That computation results in a useful amplitude and phase decomposition
of a data set of curves $f_{1},\cdots,f_{n}$ into a couple of useful
variations on the theme of data object.
\begin{itemize}
\item A Karcher mean function, $\overline{f}_{K}$, which is a useful version
of overall center. In the Step Function Example of Figure \ref{figH.E:WarpDirection},
the Karcher mean is shown as the black dotted curve in the lower left
panel. Because that curve takes into account both $\delta_{A}$ and
$\delta_{P}$ it contains both phase and amplitude notions of center
of the input curves in the top left panel. For the Shifted Betas data
in Figure \ref{figH.A:ToyBadPCA} the Karcher mean is in the right
panel of Figure \ref{figH.I:ToyVertPCA}.
\item A set of warps, $\gamma_{1,K},\cdots,\gamma_{n,K}$ of each $f_{i}$
into $\overline{f}_{K}$. These represent the horizontal variation
in the data, and can be thought of as \emph{phase data objects}\index{phase!data objects}.
Such curves comprised the phase mode of variation of the Bimodal Phase
Shift data\index{Bimodal Phase Shift data} in the lower left panel
of Figure \ref{fig:Sim-Warps}. For the Shifted Betas data set in
Figure \ref{figH.A:ToyBadPCA}, these warps are shown in the top left
of Figure \ref{figH.J:ToyWarpPCA}, where these in turn are decomposed
into further insightful modes of variation.
\item A set of aligned curves, $f_{1}\circ\gamma_{1,K},\cdots,f_{n}\circ\gamma_{n,K}$,
that represent \index{amplitude!data objects}amplitude variation.
Similarly these are usefully considered to be \emph{amplitude data
objects}. For the Step Function example, these are the colored curves
in the bottom left panel of Figure \ref{figH.E:WarpDirection}. Such
curves for the Bimodal Phase Shift data appear in the upper right
panel of Figure \ref{fig:Sim-Warps}. The top left panel of Figure
\ref{figH.I:ToyVertPCA} shows these aligned curves in the case of
the Shifted Betas example.
\end{itemize}
As noted above, PCA of the \index{amplitude!variation}amplitude data
objects results in insightful vertical modes of variation, which is
illustrated in Figure \ref{figH.I:ToyVertPCA}, based on a Fisher-Rao
analysis of the Shifted Betas data in Figure \ref{figH.A:ToyBadPCA}.
The aligned curves, $f_{1}\circ\gamma_{1,K},\cdots,f_{n}\circ\gamma_{n,K}$,
are shown, using the same colors as in Figures \ref{figH.A:ToyBadPCA}
and \ref{figH.B:ToyBadScores}, in the top left panel. As in those
figures, most panels show only the green to red curves, as the others
have disappeared due to overplotting. Their existence does appear
in the scores\index{scores} plot in the bottom right panel, which
follows a similar pattern to the PC2 scores in Figure \ref{figH.A:ToyBadPCA}.
That pattern makes intuitive sense in view of the vertical down then
up motion of the heights of the peaks in the original data set in
the top left of Figure \ref{figH.A:ToyBadPCA}. The PC1 residuals
(bottom center panel) are essentially 0 (the inexactness reflects
the discrete approximation used in the the numerical algorithm which
is also visible in the upper right panel), showing that there is only
one major vertical mode of variation in this data set.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigIToyVertPCA1}

\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigIToyVertPCA2}

\caption{Vertical analysis of the Shifted Betas data from Figure \ref{figH.A:ToyBadPCA}.
PCA of the aligned curves reveals primarily just one vertical mode
of variation. \label{figH.I:ToyVertPCA}}
\end{figure}
The nonlinear horizontal objects, which are the set of warps $\gamma_{1,K},\cdots,\gamma_{n,K}$,
can be decomposed into further modes of variation using a PCA as done
in Figure \ref{figH.J:ToyWarpPCA}. The warps (for the Shifted Betas
example) are shown in the upper left panel, again using the same common
colors. Note that the shifting peak mode of variation is nicely reflected
in the family of warps which goes from something like the convex warp
in the upper right of Figure \ref{figH.C:Warps} to the concave warp
in the lower left. It is not surprising that the mean of the warps
is $\gamma_{I}(x)=x$, as that is the criterion that was used to pick
the Karcher mean as a representer\index{representer|)} of the Fréchet
mean\index{Fréchet mean|)} equivalence class\index{equivalence class}
$\overline{\left[f\right]}$. The first mode of variation (center
left panel) of these curves appears to be some circular arcs, that
are challenging to interpret, although the color scheme, and the scores
in the center right panel, show it is clearly about the strong mode
of shifting peaks apparent in the original data. The first PC residuals
in the center panel are similarly hard to interpret. The bottom row
shows that practically all of the remaining variation appears in the
second PC component. Note that the PC2 scores (bottom right) follow
the same pattern as, i.e. are strongly correlated with, the vertical
PC scores in the bottom left of Figure \ref{figH.I:ToyVertPCA}. This
makes sense because the original curves are probability densities
(thus having area one underneath), so lower curves are wider. However,
this notion of width is hard to understand in terms of warps. The
PC2 residuals in the bottom center plot reflect the fact that there
cannot be a pure shift mode of variation in densities supported on
$[0,1]$.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigJToyWarpPCA1}

\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigJToyWarpPCA2}

\caption{PCA of warps of the Shifted Betas toy data set from Figure \ref{figH.A:ToyBadPCA}.
Shows two modes of variation, whose scores are clear, but loadings
are hard to interpret. \label{figH.J:ToyWarpPCA}}
\end{figure}
The interpretation problems with the warp modes of variation in the
left column of Figure \ref{figH.J:ToyWarpPCA} are addressed in Figure
\ref{figH.K:ToyHorizPCA}. That gives a more interpretable version
of the families of warps, by taking each curve, turning it back into
an actual warp by adding $\gamma_{I}$ to it, and then applying that
warp to the Karcher mean $\overline{f}_{K}$. The results of this,
for each $\gamma_{i,K}$, once more using the common colors are shown
in Figure \ref{figH.K:ToyHorizPCA}. Note that all of these modes
are simply warps of the Karcher mean. The left panel very clearly
shows that the first mode of variation is a shift in the peak location.
The center panel shows that the second mode is the widening that is
due to the varying heights together with the curves being densities.
The very strong connection of this mode with the dominant phase mode
in Figure \ref{figH.I:ToyVertPCA} is shown by comparing those scores
with the bottom right panel of Figure \ref{figH.J:ToyWarpPCA}. The
third panel of Figure \ref{figH.K:ToyHorizPCA} shows only the last
red curve, because it overplots all the others, since there is essentially
no variation left for the third component, which is consistent with
the PC2 residuals in the bottom center panel of Figure \ref{figH.J:ToyWarpPCA}
being essentially 0.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigKToyHorizPCA}

\caption{Horizontal modes of variation for the Shifted Betas example, shown
by using the modes of variation in Figure \ref{figH.J:ToyWarpPCA}
(added back to $\delta_{I}$) as warps applied to the Karcher mean.
Accurately and intuitively reflects two modes of horizontal variation
as expected from the raw data in the upper left of Figure \ref{figH.A:ToyBadPCA}.
\label{figH.K:ToyHorizPCA}}
\end{figure}
Fisher-Rao analysis has been quite useful in a number applications.
These include the examples following \citet{MarronRamsay2014}, including
\citet{Tucker2014EJS}, \citet{LuEJS2014Proteomics,LuEJS2014juggling,LuEJS2014SpikeTrain},
\citet{WuEJS2014SpikeTrain}, \citet{KurtekEJS2014Juggling}, \citet{StaicuLu2014AneuRisk}
and \citet{XieEJS2014AneuRisk}, as well as the modeling of seasonal
variation in electricity prices in \citet{Chen2019Elec}. Thinking
about probability distributions as \index{data objects!probability distributions}
data objects as in Section \ref{subsecBA.2.1:ProbDistns}, in situations
where the focus must be on densities and not \index{quantile}quantile
functions (shown to be preferred in the analyses of Figures \ref{figB-V:ProbDensObjects}
and \ref{figB-W:QuantFuncObjects}), Fisher-Rao analysis can also
be very useful. 

However, the Fisher-Rao analysis is not the only warp invariant approach
to the analysis of warp-equivalence quotient spaces. See Section 3.2
of the PhD dissertation \citet{Yu2017Dissertation} for introduction
to a broader class of metrics based on \emph{general root derivative}
transformation\index{general root derivative transformation}\index{transformation|)}s.
Some improvements in the sense of consistency\index{consistency}
(as $n\rightarrow\infty$) are proposed there. However, the Fisher-Rao
approach remains a workhorse method because of its excellent properties
in terms of peak alignment in curve registration, as shown in Figure
\ref{figH.G:NoPinch}.

The value of the $L^{2}$ norm applied to SRVFs, as in (\ref{eqH.E:AmplitudeDistance}),
for peak alignment is demonstrated in Figure \ref{figH.G:NoPinch},
based on an example constructed by Xiaosun Lu. That is based on two
Gaussian curves (not densities as they are just vertical rescalings),
$f_{1}$ (dashed) and $f_{2}$ (solid), shown in blue in the bottom
center panel. The left side of the top panel compares $2\left\Vert f_{1}-f_{2}\circ\gamma_{\alpha}\right\Vert $
(solid black curve) with $\left\Vert q_{f_{1}}-q_{f_{2}\circ\gamma_{\alpha}}\right\Vert _{2}$
(dash dotted curve) for a family of warps indexed by $\alpha\in\left(-1,0\right]$.
The factor of 2 is used to put these on approximately the same vertical
scale. The family of warps (one member of which appears as the solid
red piece-wise line in the middle left panel) is piece-wise linear,
with a strong vertical jump of size $\alpha$ over a very small interval
near the center. Those warps squish the interval $\left(\frac{1-\alpha}{2},\frac{1+\alpha}{2}\right)$
into a very small interval, so that $f_{2}\circ\gamma_{\alpha}$ (shown
as the solid curve in the bottom left) suffers what has been called
a ``pinching effect''. The red line in the top panel shows the value
of $\alpha$ that minimizes $2\left\Vert f_{1}-f_{2}\circ\gamma_{\alpha}\right\Vert $,
and the corresponding warp $\gamma_{\alpha}$ is the solid red curve
in the middle left panel. It is clear from comparing the red curves
$f_{1}$ (dashed) and $f_{2}\circ\gamma_{\alpha}$ in the bottom left
panel that they are visually much closer together than are the blue
curves (corresponding to $\alpha=0$) in the bottom center panel,
which is quite concerning for peak alignment methods based on the
$L^{2}$ norm between functions. This comparison is quantified in
the top panel, where the solid curve is lower at the red line than
at the blue line). However, note that the norm of the corresponding
SRVFs $\left\Vert q_{f_{1}}-q_{f_{2}\circ\gamma_{\alpha}}\right\Vert _{2}$
(dot dashed curve) has the opposite behavior. This is because the
red solid curve has very steep derivatives near the center of the
plot, so the $L^{2}$ norm of the SRVFs strongly penalizes against
such pinching. 

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterH/OODAbookChpHFigGNoPinch}

\caption{Top row shows the $L^{2}$ distance\index{distance|)} between original
functions $f_{1}$ and $f_{2}$ (solid) versus SRVFs $q_{1}$ and
$q_{2}$ (dot-dashed) for a family of pinching warps (red) indexed
by $\alpha\in(-1,0]$ in the left hand column and a family of stretching
warps (green) indexed by $\alpha\in[0,1)$ in the right column. The
middle column shows the $\alpha=0$ identity warp\index{identity warp|)}
in blue and the original $f_{1}$, $f_{2}$ below. Shows how the $L^{2}$
norm applied to SRVFs gives a much more natural notion of peak alignment.
\label{figH.G:NoPinch}}
\end{figure}
The left side of the top panel does a slightly different comparison,
this time of $2*\left\Vert f_{1}\circ\gamma_{\alpha}-f_{2}\right\Vert _{2}$
versus $\left\Vert q_{f_{1}\circ\gamma_{\alpha}}-q_{f_{2}}\right\Vert _{2}$
for a different family of piece-wise warps $\gamma_{\alpha}$, for
$\alpha\in\left[0,1\right)$. That warping family has a flat spot
of length $\alpha$ in the center, which stretches (instead of compressing
as for $\alpha<0$) the width of the peak. The warp is applied to
$f_{1}$ instead of $f_{2}$ in this case because that gives smaller
norm values. The solid norm of the functions in the top panel shows
a local minimum indicated by the green line. That local minimum is
intuitively understood by noticing the green curves in the bottom
right panel are again closer than shown in the blue panel. Once again
applying the norm to the SRVFs results in a smaller value at $\alpha=0$
(blue curves) this time because the warped dashed curve has a too
small derivative near the center.

\citet{WagnerKneip2019} point out that a weakness of the Fisher-Rao
numerical approach is that when $\overline{f}_{K}\equiv0$ the method
will fail. Essentially the dynamic programming algorithm requires
a non-constant template mean for convergence. They propose a method
that does work in such cases, although it is not as efficient at peak
alignment in situations where that is a priority. Other approaches
based on the SRVF include Bayesian functional alignment, where prior
distributions are placed on the warping functions. \citet{Chengetal16Bayesian}
use a Dirichlet prior on the simplex and \citet{Luetal17gaussianprior}
set a Gaussian process prior on the sphere for the warp functions.

\section{Principal Nested Spheres Decomposition\label{secH.B:PNSonFisherRao}}

While the interpretation of the Fisher-Rao approach to curve registration\index{curve registration|)}
has already benefited in a major way from OODA concepts such as the
warp equivalence class\index{equivalence class} quotient space\index{quotient space|)},
in some situations there is also technical benefit from other OODA
ideas. In particular, the analysis of the warping functions done in
(\ref{eqH.F:WarpsOnSphere}) shows that the warps $\gamma_{1,K},\cdots,\gamma_{n,K}$,
which quantify horizontal variation in a Fisher-Rao analysis, all
have corresponding SRVFs that lie on the surface of a high dimensional
sphere in function space. The decomposition into modes of variation
in Figure \ref{figH.J:ToyWarpPCA} was essentially a Euclidean PCA
done in the tangent space centered at the Karcher mean\index{Karcher mean|)}
$\overline{f}_{K}$. As noted in Section \ref{secG-A:ShapeDataObjects},
this can result in the problem of a distorted analysis which can be
substantially improved using PNS (recall Principal Nested Spheres\index{Principal Nested Spheres}
from Section \ref{secG-A:ShapeDataObjects}) on the corresponding
SRVFs. This motivated the proposal of \citet{YuLuMarron2017} to couple
PNS with a Fisher-Rao analysis to understand phase variation.

In many situations, there is not a major difference between the tangent
plane analysis and PNS. This seems to be because for many data sets
SRVFs do not tend to spread very broadly around the sphere (so the
tangent plane approximation is frequently adequate). This may be due
to the fact that for $\gamma(t)$ increasing, $\psi_{\gamma}(t)>0$,
and hence $\gamma$ lies in the positive orthant of the sphere, which
in finite dimensions is $\left\{ \boldsymbol{x}\in\mathbb{R}^{d}:\left\Vert \boldsymbol{x}\right\Vert _{2}=1,x_{1}>0,\cdots,x_{d}>0\right\} $.
However a toy example showing the tangent plane approximation can
be a serious problem is shown in Figure \ref{figH.L:ToyNoSphere}.

\begin{figure}[H]
\includegraphics[scale=0.35]{ChapterH/OODAbookChpHFigLToyNoSphere}

\caption{Far left panel shows a family of piece-wise linear warps $\gamma(t)$.
Center left shows representations of each warp in terms of heights
at $t=\frac{1}{3},\frac{2}{3}$, together with projections onto Euclidean
PCA. Center right is reconstruction of Euclidean projections with
cyan curves highlighting invalid warping functions. Far left shows
PNS projections, which are all valid warps. \label{figH.L:ToyNoSphere} }
\end{figure}
Figure \ref{figH.L:ToyNoSphere} first appeared in the PhD dissertation
of \citet{Lu2013Dissertation}, and is Figure 5 of \citet{YuLuMarron2017}.
The far left panel contains a set of warping functions $\gamma(t)$,
which are piece-wise linear, all having just two breakpoints at $t=\frac{1}{3}$
and $\frac{2}{3}$. Thus each warp is entirely characterized by $\left(\gamma(\frac{1}{3}),\gamma(\frac{2}{3})\right)$,
and is shown as a small black circle in the center left panel. Because
warps must be increasing, i.e. $\gamma(\frac{1}{3})<\gamma(\frac{2}{3})$,
all of these circles are above the diagonal black line that represents
$\gamma(\frac{1}{3})=\gamma(\frac{2}{3})$. The results of a direct
Euclidean PCA on the family of warps is also overlaid on the center
left panel. The red line shows the first eigen direction, and the
symbols are the projections of the warps onto the first mode of variation.
Red plus signs are those projections that are actually warps (i.e.
are increasing functions), while the cyan circles indicate those of
the projections that are actually outside of the set of warps (because
they are no longer increasing functions). The impact of replacing
the Euclidean PCA of the warping\index{warping|)} functions with
a PNS of the corresponding SRVF\index{Square Root Velocity Function|)}s
is shown in the right panel. Note that each projection onto the first
PNS component remains a valid warp. Another view of the effect of
PNS is the blue curve shown in the center left panel. That is a plot
of $\left(\gamma(\frac{1}{3}),\gamma(\frac{2}{3})\right)$ for a dense
grid of points lying along the PNS one dimensional representation.
It shows how appropriately managing the curvature of the space keeps
that analysis in the range of functions of interest. A similar issue
occurred when studying covariance matrices in Section \ref{subsec:Covariance-matrix-distances}
where linear extrapolation using the Euclidean \index{metric|)}metric\index{metric!Euclidean}
can lead to negative eigenvalues, hence leaving the space of covariance
matrices (which are all positive semi-definite). Recall that challenge
was similarly addressed by using a curved manifold\index{manifold}
for representing covariance matrix \index{data objects!covariance matrices}data
objects. 

Figure \ref{figH.M:ToyWarpPNS} shows another example of the improvement
available from PNS of the Fisher-Rao warp function\index{warp function|)}s
over the simple tangent plane PCA shown in Figure \ref{figH.J:ToyWarpPCA}.
That is based on the Shifted Betas data\index{Shifted Betas data|)}
in the top left panel of Figure \ref{figH.A:ToyBadPCA} which was
used with the same coloring to illustrate many points in Section \ref{secH.A:Fisher-Rao-Reg}.
The input to the PNS is the set of Fisher-Rao warps, $\gamma_{1,K},\cdots,\gamma_{n,K}$,
shown in the top left panel of Figure \ref{figH.J:ToyWarpPCA}. Figure
\ref{figH.M:ToyWarpPNS} shows various projections onto the rank 2
Euclidean PCA subspace. The circles are the (rank 2) projections of
the warps onto this space, with the horizontal axis showing the first
PC eigen direction, and the vertical axis the second. The x signs
show the first PCA mode of variation (shown in a different way in
the middle left panel of Figure \ref{figH.J:ToyWarpPCA}), whose constraint
of being linear makes them a poor approximation of the clearly nonlinear
data set. That strong distortion is addressed by the plus signs, which
are the projections of the first PNS approximation of this data set.
That mode of variation is essentially nonlinear with respect to this
data view. It clearly provides a far better approximation of the data.
A careful look at the horizontal and vertical scales shows this effect
in not terribly strong (also reflected by the rather small variation
in the center panel of Figure \ref{figH.K:ToyHorizPCA}) although
still considerable. A more dramatic example (which is based on even
stronger phase variation\index{phase!variation|)}) appears in Section
2 of \citet{YuLuMarron2017}.\index{modes of variation!amplitude and phase|)}

\begin{figure}[H]
\includegraphics[scale=0.5]{ChapterH/OODAbookChpHFigMToyWarpPNS}

\caption{Shows benefits of a PNS analysis of the warps in the top left panel
of Figure \ref{figH.J:ToyWarpPCA} whose projections onto the 2-d
PCA subspace appear as circles. The x signs are the rank 1 PCA approximations.
Much better approximation comes from the projections of the rank 1
PNS approximation shown as plus signs. \label{figH.M:ToyWarpPNS}}
\end{figure}
A real data example, where PNS gives a big benefit over tangent plane
PCA\index{Principal Component Analysis|)} can be found in Section
3 of \citet{YuLuMarron2017}, in the context of real time monitoring
of blood glucose. Additional examples appear in \citet{StaicuLu2014AneuRisk}
and in the analysis of the Juggling data\index{Juggling data} by
\citet{LuEJS2014juggling} shown in Figure \ref{fig:Juggling}.\index{data objects!amplitude - phase|)}
\index{data objects!curves|)} \index{functional data analysis|)}
