
\chapter{PCA Details and Variants\label{chapP:PCA-Details}}

PCA\index{Principal Component Analysis|(} (recall from Section \ref{sec:Curves-as-Data},
this is an acronym for Principal Component Analysis) has been very
useful at many points in preceding chapters. A good source for many
aspects of PCA (including the history, such as the name itself going
back to \citet{Hotelling1933PCA}) is \citet{Jolliffe2002}. Common
terminology for PCA is ``\index{dimension reduction}dimension reduction'',
which appears to be rooted in the era of pre-computational statistics,
when indeed dimension reduction was a critical task. In the current
era of abundant computational capabilities, more relevant and specific
ways of thinking about PCA are:
\begin{itemize}
\item \emph{Data \index{visualization|(}visualization}. For example where
scores distribution plots\index{scores!distribution plot} clearly
show the \emph{relationships} between\index{data objects!relationships between}
data objects. This was done in preceding chapters at many points,
starting with the Spanish Mortality data\index{Spanish Mortality data}
in Figure \ref{fig:Mortality-Scatterplot} and the \index{Twin Arches data}Twin
Arches toy data in Figure \ref{fig:50d-ToyPCAScatMat}. Further insights
came from decomposition into easily interpretable one dimensional
modes of variation\index{modes of variation|(} (each revealing both
loadings\index{loadings} and scores), starting again with Spanish
Mortality in Figures \ref{fig:Mortality-PC1} - \ref{fig:Mortality-PC2},
and with the 2-d Toy data\index{2-d Toy data} in Figures \ref{fig:2dToy-PC1}
- \ref{fig:2dToy-Decomp}. 
\item \emph{Denoising\index{denoising}}. Useful in situations where a low
rank approximation of the data can retain most of the signal and thus
PCA may reduce noise in the data.
\item \emph{Efficient data representation\index{data representation|(}}.
Which comes up in HDLSS\index{High Dimension Low Sample Size} situations
where $d>n$ as discussed in Section \ref{secM-B:HDLSS}. The $(n-1)\times n$
matrix that contains all nonzero PC \index{scores|(}scores is essentially
a rotation of the mean centered\index{mean!center|(} data into the
lower dimensional space $\mathbb{R}^{n}$, giving an $(n-1)$ dimensional
(exact) representation of the data, as detailed in Section \ref{subsecP-A-D:PCA-Comput}. 
\end{itemize}
Several quite different but useful ways of thinking about PCA are
given in Section \ref{secP-A:PCA-Views}. Section \ref{secP-B:TwoBlockDecomp}
discusses methods that parallel PCA in the sense of shifting the focus
from variation in a single data block to analyzing \emph{joint} variation
between two different but \emph{linked} blocks of data.

Often really good ideas are discovered several times independently,
usually under a different name. PCA is a clear example of this, with
alternate names appearing in other fields as listed in Table \ref{tabP-A:AlternatePCAnames}.

\begin{table}[H]
\begin{tabular}{|c|c|}
\hline 
Name & Field\tabularnewline
\hline 
\hline 
Factor Analysis & Social Sciences\tabularnewline
\hline 
Karhunen Loeve Expansion & Probability Theory / Electrical Engineering\tabularnewline
\hline 
\index{Empirical Orthogonal Functions}Empirical Orthogonal Functions
(EOF) & Geo and Climate Sciences\tabularnewline
\hline 
\index{Proper Orthogonal Decomposition}Proper Orthogonal Decomposition
(POD) & Applied Mathematics\tabularnewline
\hline 
\end{tabular}

\caption{\index{PCA!alternate names}Alternate names for concepts essentially
similar to PCA in other fields. \label{tabP-A:AlternatePCAnames}}
\end{table}
Not all of these are direct synonyms for PCA. For example, as discussed
in Section \ref{subsecP-A-C:Gauss-Likelihood}, in psychometrics factor
analysis\index{factor analysis} is not exactly PCA, as it includes
simultaneous likelihood estimation of the residual variance. However,
the name factor analysis is also used in a number of other fields
where it actually is a synonym for PCA. Also, the \index{Karhunen-Loeve expansion}Karhunen-Loeve
expansion is the same idea, but typically applied to probability distributions
instead of to data sets. Of course these can be related by thinking
of a data set in terms of its empirical discrete probability distribution,
using the idea of \emph{empirical probability measure}\index{empirical probability measure},
that puts probability mass $n^{-1}$ on each data object. 

\section{Viewpoints of PCA\label{secP-A:PCA-Views}}

As discussed above there are many ways of considering PCA, most of
which provide differing useful insights. As noted in Section \ref{secF.E:ComboView-Centering}
\emph{data centering\index{centering|(}} issues tend to be downplayed,
but actually have a surprisingly important and sometimes non-intuitive
role (see for example the combined view\index{combined view|(} analysis
of the Twin Arches data in Figures \ref{figF.F1:ComboViewInputData}
and \ref{figF.F2:ComboViewOAMeanResid}). This is true for both discussions
about rows or columns of the data matrix\index{data matrix} being
viewed as the data objects (recall the discussion in Section \ref{secBA.1:Foundations_Terminology}),
and also in data visualizations, as studied below in Section \ref{subsecP-A-A:Data-Centering}.
Good insight into many aspects of PCA, including the decomposition
of the column object centered data matrix into insightful modes of
variation, comes from studying the Singular Value Decomposition\index{Singular Value Decomposition|(}
(SVD) in Section \ref{subsecP-A-B:SVD}. A traditional and important
approach to PCA is through a low rank factor model whose parameters
are estimated by Gaussian likelihood estimation (which can be deceptive)
as reviewed in Section \ref{subsecP-A-C:Gauss-Likelihood}. Relevant
computational and graphical issues are discussed in Section \ref{subsecP-A-D:PCA-Comput}.

Recalling notation from Section \ref{secBA.2:OverviewObjFeat}, and
the use of tildes to indicate random quantities from Table \ref{tab7-1:TildeNotation},
write a $d\times n$ data matrix as
\begin{equation}
\widetilde{\boldsymbol{X}}=\left[\begin{array}{ccc}
\widetilde{x}_{1,1} & \cdots & \widetilde{x}_{1,n}\\
\vdots & \ddots & \vdots\\
\widetilde{x}_{d,1} & \cdots & \widetilde{x}_{d,n}
\end{array}\right]=\left[\begin{array}{ccc}
\widetilde{\boldsymbol{x}}_{1} & \cdots & \widetilde{\boldsymbol{x}}_{n}\end{array}\right],\label{eqP-A:DataMatrix}
\end{equation}
where the data object column vectors (sometimes called ``cases''
or even ``samples'') are
\[
\widetilde{\boldsymbol{x}}_{j}=\left[\begin{array}{c}
\widetilde{x}_{1,j}\\
\vdots\\
\widetilde{x}_{d,j}
\end{array}\right]\in\mathbb{R}^{d},
\]
and where each number $\widetilde{x}_{i,j}$ is the trait value (sometimes
called ``variable'' or ``feature'') for $j=1,\cdots,n$, and $i=1,\cdots,d$.
Linear algebra provides a useful mathematical backbone for understanding
many aspects of OODA, and an important issue is the space in which
this is done. Many mathematical treatments focus on the vector spaces
$\mathbb{R}^{d}$ and/or $\mathbb{R}^{n}$, but in this chapter we
take the unconventional approach of developing the ideas in the vector
space $\mathbb{R}^{d\times n}$. Recall that space from (\ref{eqB.E:Matrices})
where the symbol $\times$ in the exponent is used to indicate the
vector space of matrices, in contrast with $\mathbb{R}^{dn}$, which
denotes the corresponding vector space of long vectorize\index{vectorize}d
column vectors of length $dn$. Subspaces and projections in $\mathbb{R}^{d\times n}$
provide a simple framework for understanding the relationships between
many OODA operations. This approach is a simplification of the ideas
that \citet{Mammen2001SmoothingAsProjection} developed to understand
how kernel smoothing\index{smoothing} (discussed in Sections \ref{secN-A:NoHisto}
and \ref{secN-B:SmoothBasics}) can be viewed as a projection.

Recall from (\ref{eqB.F:IdentityMatrix}) that $\boldsymbol{I}_{d}$
denotes the $d\times d$ identity matrix\index{identity matrix},
from (\ref{eq:ones matrix}) that $\boldsymbol{1}_{d,n}$ is the $d\times n$
matrix of ones, and from (\ref{eqB.G:0sMatrix}) that $\boldsymbol{0}_{d,n}$
indicates the $d\times n$ matrix of $0$s.

The operation of \emph{projection\index{projection|(}} is fundamental
here. Given any metric\index{metric!space} space with distance\index{distance}
$\delta$, a subset $\mathcal{S}$, and an element $\boldsymbol{x}$,
the projection of $\boldsymbol{x}$ onto $\mathcal{S}$ is the closest
point in $\mathcal{S}$ to $\boldsymbol{x}$, i.e.
\[
P_{\mathcal{S}}\left(\boldsymbol{x}\right)=\arg\min_{\boldsymbol{s}\in\mathcal{S}}\delta\left(\boldsymbol{s},\boldsymbol{x}\right).
\]
Centering operations are studied in Section \ref{subsecP-A-A:Data-Centering},
using projections in $\mathbb{R}^{d\times n}$, with respect to the
Frobenius norm\index{Frobenius norm} $\left\Vert \cdot\right\Vert _{F}$,
defined in (\ref{eq:Frobeniusnorm}). Additional useful linear algebra
comes from the corresponding Frobenius inner product\index{Frobenius inner product}
\begin{equation}
\left\langle \boldsymbol{M},\boldsymbol{N}\right\rangle _{F}=\sum_{i=1}^{d}\sum_{j=1}^{n}M_{i,j}N_{i,j},\label{eqP.AB:FrobeniusInnerProd}
\end{equation}
on $\mathbb{R}^{d\times n}$.

Two relevant subspaces of $\mathbb{R}^{d\times n}$ are based on the
concept of \emph{flat vector\index{flat vector}s} whose entries are
all the same. For example in $\mathbb{R}^{d}$, the set of flat vectors
is the subspace 
\[
\left\{ \left(\begin{array}{c}
u\\
\vdots\\
u
\end{array}\right):u\in\mathbb{R}\right\} ,
\]
which is sometimes called ``the 45 degree line'', although that
only makes sense in the case $d=2$. The flat vectors are also sometimes
called ``constant vectors'', but the ``flat'' terminology seems
easier to keep in mind in the context of the matrix space $\mathbb{R}^{d\times n}$. 

The concept of flat vectors allows interpretation of conventional
univariate means as projections. First note that a unit basis vector
(i.e. with norm 1) of the subspace of flat vectors in $\mathbb{R}^{n}$
is $\frac{1}{\sqrt{n}}\boldsymbol{1}_{n,1}$. Hence given $\boldsymbol{x}=\left(\begin{array}{c}
x_{1}\\
\vdots\\
x_{n}
\end{array}\right)\in\mathbb{R}^{n}$, the $\mathbb{R}^{n}$ projection coefficient of $\boldsymbol{x}$
onto the flat subspace is the inner product $\frac{1}{\sqrt{n}}\boldsymbol{1}_{1,n}\boldsymbol{x}=\sqrt{n}\overline{x}$.
Multiplying that coefficient by the basis vector results in the $\mathbb{R}^{n}$
projection\index{projection|)} of $\boldsymbol{x}$ onto the subspace
of flat vectors as 
\begin{equation}
\frac{1}{\sqrt{n}}\boldsymbol{1}_{n,1}\left(\frac{1}{\sqrt{n}}\boldsymbol{1}_{1,n}\boldsymbol{x}\right)=\frac{1}{n}\boldsymbol{1}_{n,n}\boldsymbol{x}=\boldsymbol{1}_{n,1}\overline{x},\label{eqP.AA:UniMeanProj}
\end{equation}
which is the flat vector whose common entry is $\overline{x}$.

To understand data centering as \index{mean!as projection}projection
in the matrix space $\mathbb{R}^{d\times n}$, useful notation is
$\mathcal{S}_{FT}$ for the subspace (of $\mathbb{R}^{d\times n}$)
in which all of the row trait\index{row trait!vector|(}\index{trait|(}
vectors are flat, where $FT$ stands for \emph{Flat Traits}\index{flat traits},
i.e. 
\[
\mathcal{S}_{FT}=\left\{ \boldsymbol{u}\boldsymbol{1}_{1,n}:\boldsymbol{u}\in\mathbb{R}^{d}\right\} .
\]
Similarly define the subspace (of $\mathbb{R}^{d\times n}$) consisting
of matrices composed of flat column object vectors (\emph{Flat Objects}\index{flat objects})
to be 
\[
\mathcal{S}_{FO}=\left\{ \boldsymbol{1}_{d,1}\boldsymbol{v}^{t}:\boldsymbol{v}\in\mathbb{R}^{n}\right\} .
\]
Applying the transpose of (\ref{eqP.AA:UniMeanProj}) to each row
of the data matrix $\widetilde{\boldsymbol{X}}$ gives the projection
\[
P_{\mathcal{S}_{FT}}\left(\widetilde{\boldsymbol{X}}\right)=\widetilde{\boldsymbol{X}}\left(\frac{1}{\sqrt{n}}\boldsymbol{1}_{n,1}\right)\left(\frac{1}{\sqrt{n}}\boldsymbol{1}_{1,n}\boldsymbol{x}\right)=\widetilde{\boldsymbol{X}}\left(\frac{1}{n}\boldsymbol{1}_{n,n}\right)=\overline{\boldsymbol{x}}_{CO}\boldsymbol{1}_{1,n},
\]
which essentially extends the $d\times1$ column object mean\index{column object mean|(}
$\overline{\boldsymbol{x}}_{CO}$ from (\ref{eqF1:ColumnObjectMean})
into a $d\times n$ element of $\mathcal{S}_{FT}$. The column object
centered version of the data matrix can also be written as a projection
onto the orthogonal (with respect to the Frobenius inner product (\ref{eqP.AB:FrobeniusInnerProd}))
complementary subspace $\mathcal{S}_{FT}^{\perp}=\left\{ \boldsymbol{M}\in\mathbb{R}^{d\times n}:\boldsymbol{M}\perp\mathcal{S}_{FT}\right\} $
as, 
\[
P_{\mathcal{S}_{FT}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)=\widetilde{\boldsymbol{X}}-\overline{\boldsymbol{x}}_{CO}\boldsymbol{1}_{1,n}=\widetilde{\boldsymbol{X}}\left(\boldsymbol{I}_{n}-\frac{1}{n}\boldsymbol{1}_{n,n}\right).
\]
Similarly, row trait\index{row trait!mean|(} mean centering as defined
in (\ref{eqF2:RowFeatureMean}) can be studied in terms of projections
onto the flat object subspace $\mathcal{S}_{FO}$ 
\[
P_{\mathcal{S}_{FO}}\left(\widetilde{\boldsymbol{X}}\right)=\boldsymbol{1}_{d,1}\overline{\boldsymbol{x}}_{RT}^{t}=\left(\frac{1}{d}\boldsymbol{1}_{d,d}\right)\widetilde{\boldsymbol{X}},
\]
and its orthogonal complement $\mathcal{S}_{FO}^{\perp}$ 
\[
P_{\mathcal{S}_{FO}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)=\widetilde{\boldsymbol{X}}-\boldsymbol{1}_{d,1}\overline{\boldsymbol{x}}_{RT}^{t}=\left(\boldsymbol{I}_{d}-\frac{1}{d}\boldsymbol{1}_{d,d}\right)\widetilde{\boldsymbol{X}}.
\]


\subsection{Data Centering\label{subsecP-A-A:Data-Centering}}

Some mean centering issues were discussed in Sections \ref{secF.E:ComboView-Centering}
and \ref{subsecI-C-C:Procrustes}. As noted in \citet{Prothero2021Centering},
while centering seems like a simple issue, it can be surprisingly
slippery. Here a deeper view based on the above ideas of projections
onto $\mathbb{R}^{d\times n}$ is based on the Sine Wave\index{Sine Wave data|(}
toy data set, with $d=20$ rows and $n=10$ columns, shown in Figure
\ref{figP-A1:Centering-Raw}. The data were generated (for $i=1,\cdots,20$
and $j=1,\cdots,10$) as 
\begin{equation}
x_{i,j}=T_{1}+T_{2}+T_{3}-3+N(0,10^{-6})\label{eqP-B:CenteringToyDataGen}
\end{equation}
where 
\[
T_{1}=\sin\left(5\cdot\pi\cdot\left(i-1\right)/19\right),
\]
\[
T_{2}=0.3\cdot\left(j-5.8\right)^{2},
\]
\[
T_{3}=0.005\cdot\left(i-10.5\right)\cdot\left(j-5.5\right).
\]
 Let $\boldsymbol{T}_{1}$, $\boldsymbol{T}_{2}$ and $\boldsymbol{T}_{3}$,
denote the matrix versions of the first three terms on the right hand
side of (\ref{eqP-B:CenteringToyDataGen}). Note that each of these
is a \emph{mode of variation} as defined in Section \ref{secBA.1:Foundations_Terminology}.
These modes illustrate how column object and row trait mean centering
are understood through projections in $\mathbb{R}^{d\times n}$. In
particular, note that $\boldsymbol{T}_{1}=\boldsymbol{u}\boldsymbol{1}_{1,n}$
($\boldsymbol{u}\in\mathbb{R}^{d}$ generates the sine wave) has flat
rows, so $\boldsymbol{T}_{1}\in\mathcal{S}_{FT}$ (i.e. is a flat
trait mode) and that the variation across rows follows a sin wave.
Similarly $\boldsymbol{T}_{2}=\boldsymbol{1}_{d,1}\boldsymbol{v}^{t}$
($\boldsymbol{v}\in\mathbb{R}^{n}$ gives the parabola) has flat columns,
so $\boldsymbol{T}_{2}\in\mathcal{S}_{FO}$ and its (flat object)
variation across columns follows a parabola. Both of these modes are
clearly visible in the views of this Sine Wave data set shown in Figure
\ref{figP-A1:Centering-Raw}. The third mode $\boldsymbol{T}_{3}$
is a product of two (both mean $\boldsymbol{0}$) linear vectors and
has a much smaller coefficient to make it a much smaller contribution
to the overall variation. The linear factors have both been carefully
centered, which entails $\boldsymbol{T}_{3}\in\mathcal{S}_{FT}^{\perp}\cap\mathcal{S}_{FO}^{\perp}$
which will leave $\boldsymbol{T}_{3}$ remaining after both centering
operations. Finally a very low level of independent Gaussian variation
is also added. 

Figure \ref{figP-A1:Centering-Raw} explores the Sine Wave data set
of (\ref{eqP-B:CenteringToyDataGen}) using a different formatting
of the combined \index{data matrix!combined view}data\index{combined view|)}
matrix view introduced in Section \ref{secF.E:ComboView-Centering}.
A \index{heat-map|(}heat-map view (as discussed in Sections \ref{secF.A:Heatmap-Views})
is shown in the left panel. As expected from the form of (\ref{eqP-B:CenteringToyDataGen}),
the heat-map contains both vertical sine wave and horizontal parabolic
patterns. 

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterP/OODAbookChpPFigAA1ToyRaw}

\caption{Combined view of the Sine Wave raw data matrix. Heat-map in the left
panel. Curve views of same data set, with columns as data objects
in the center panel and row trait vectors on the right. Visual connections
between the three different views takes some effort to make. \label{figP-A1:Centering-Raw}}
\end{figure}
Another view of the Sine Wave data matrix (in the spirit of the top
right panels of Figures \ref{figF.F1:ComboViewInputData}-\ref{figF.F4:ComboViewMode2})
is of the column data objects (essentially sine waves), shown as $n=10$
curves in the center panel. That view suggests that vertical shift
of the (common) sine pattern is the dominant mode of variation. A
third view highlights the $d=20$ row trait vectors in $\mathbb{R}^{n}$
shown as curves in the right panel. That graphic is in the spirit
of the lower left panels of Figures \ref{figF.F1:ComboViewInputData}-\ref{figF.F4:ComboViewMode2}.
It gives the quite different impression that the dominant mode of
variation is instead about vertical shift of the parabolic curve.
A careful look at the matrix view in the left panel reveals what is
happening: the variation in the height of the sine waves follows the
parabola, and the variation in the heights of the parabola follows
the sine wave pattern. An important point is \emph{neither} is a single
mode of variation, as defined in Section \ref{subsecBA.1.1:FormalizingModes},
as they are not rank one matrices. In fact both are the sum of \emph{two}
modes shown as $\boldsymbol{T}_{1}+\boldsymbol{T}_{2}$ in (\ref{eqP-E:OverallMeanCentData}).
The first mode $\boldsymbol{T}_{1}=\boldsymbol{u}\boldsymbol{1}_{1,n}$
can be thought of as either extending the sine wave $\boldsymbol{u}$
in a flat horizontal way (giving the common sine wave in the center
panel) or else providing sinusoidal variation to the flat line $\boldsymbol{1}_{1,n}$.
The second mode $\boldsymbol{T}_{2}=\boldsymbol{1}_{d,1}\boldsymbol{v}^{t}$
similarly plays a dual role. It both provides parabolic variation
(modeled by $\boldsymbol{v}^{t}$) to the vertical flat line $\boldsymbol{1}_{d,1}$
to drive the variation in the center panel, and also gives a vertical
flat extension of the parabola $\boldsymbol{v}^{t}$ apparent in the
right panel. 

The effect of column object mean centering, i.e. $P_{\mathcal{S}_{FT}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)$,
is studied in Figure \ref{figP-A2:Centering-Col}. Because $\boldsymbol{T}_{1}\in\mathcal{S}_{FT}$
that subtraction of the column data object mean essentially removes
the vertical sine wave from the heat map (left panel of Figure \ref{figP-A2:Centering-Col})
making the horizontal parabolic structure of the mode $\boldsymbol{T}_{2}$
even more apparent. The conventional column data object view in the
middle panel shows that the common structure of the curves in the
middle of Figure \ref{figP-A1:Centering-Raw}, i.e. the sine wave
structure, has been removed, leaving mostly just the flat vertical
shift mode of variation $\boldsymbol{T}_{2}$ (with heights determined
by the parabola). Also note that the linear $\boldsymbol{T}_{3}$
mode of variation begins to be visible in these curves as well. This
same type of removal of mean structure was very useful in the Spanish
Mortality data\index{Spanish Mortality data} analysis of Figure \ref{fig:Mortality-Mean-Resid}
where the column object mean nicely captured the mortality structure
that is common over time, while the mean residuals contained the variation
across time. It also was demonstrated using the Tilted Parabolas data\index{Tilted Parabolas data}
in the top row of Figure \ref{fig:10d-ToyPCA}. The alternate row
trait view (right panel of Figure \ref{figP-A2:Centering-Col}) instead
shows that column object mean centering removes most of the $\boldsymbol{T}_{1}$
driven variation visible in the right panel of Figure \ref{figP-A1:Centering-Raw},
again since that variation was driven by the sinusoidal structure
in the data, ending up with the set of parabolas that are nearly all
the same in $\boldsymbol{T}_{2}$.

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterP/OODAbookChpPFigAA2ToyColObjCent}

\caption{Effect of column object centering for the Sine Wave data of Figure
\ref{figP-A1:Centering-Raw}, using the same panel format. Shows removal
of vertical sine wave component, leaving nearly constant columns whose
heights follow the parabola. \label{figP-A2:Centering-Col}}
\end{figure}
Figure \ref{figP-A3:Centering-Row} shows the effect of complementary
row trait mean centering, i.e. $P_{\mathcal{S}_{FO}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)$,
on the Sine Wave data set. The heatmap view in the left panel shows
that this time the horizontal parabola effect $\boldsymbol{T}_{2}$
(defined in (\ref{eqP-B:CenteringToyDataGen})) has been removed,
leaving mostly horizontal stripes following the $\boldsymbol{T}_{1}$
vertical sine wave pattern. In the conventional column data object
curve view in the center panel, this projection has the effect of
removing the dominant and important $\boldsymbol{T}_{2}$ vertical
shift mode of variation (since $\boldsymbol{T}_{2}\in\mathcal{S}_{FO}$).
The vertical axes in the center and right panels of Figures \ref{figP-A1:Centering-Raw},
\ref{figP-A2:Centering-Col} and \ref{figP-A3:Centering-Row} all
use a common scale to show the $\boldsymbol{T}_{2}$ parabolic variation
has the larger magnitude relative to the $\boldsymbol{T}_{1}$ sine
wave component. A major point illustrated here is that the type of
centering used can have a strong visual impact, although frequently
too little attention is paid to this point. Finally the right panel
shows that this operation removes the $\boldsymbol{T}_{2}$ parabolic
structure from the corresponding view of the original data on the
right of Figure \ref{figP-A1:Centering-Raw}. Note that the $\boldsymbol{T}_{3}$
tilted linear component is now a little more apparent (actually in
both plots).

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterP/OODAbookChpPFigAA3ToyRowObjCent}

\caption{Row trait centering of Sine Wave data from Figure \ref{figP-A1:Centering-Raw},
with the same panel format. Shows removal of $\boldsymbol{T}_{2}$
horizontal parabolic component.\label{figP-A3:Centering-Row}}
\end{figure}
Another data centering option is \emph{double centering\index{double centering}\index{centering!double}}
shown in Figure \ref{figP-A4:Centering-Overall}. Double centering
can be viewed in a number of ways. The simplest is through simultaneous
removal of both types of mean via a composition of projections $P_{\mathcal{S}_{FT}^{\perp}}\left(P_{\mathcal{S}_{FO}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)\right)$.
Other ways of thinking about (or even recipes for computation of)
double centering include 
\[
P_{\mathcal{S}_{FT}^{\perp}}\left(P_{\mathcal{S}_{FO}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)\right)=\left(\widetilde{\boldsymbol{X}}-P_{\mathcal{S}_{FO}}\left(\widetilde{\boldsymbol{X}}\right)\right)-P_{\mathcal{S}_{FT}}\left(\widetilde{\boldsymbol{X}}-P_{\mathcal{S}_{FO}}\left(\widetilde{\boldsymbol{X}}\right)\right)=
\]
\[
=\widetilde{\boldsymbol{X}}-P_{\mathcal{S}_{FO}}\left(\widetilde{\boldsymbol{X}}\right)-P_{\mathcal{S}_{FT}}\left(\widetilde{\boldsymbol{X}}\right)+P_{\mathcal{S}_{FT}\cap\mathcal{S}_{FO}}\left(\widetilde{\boldsymbol{X}}\right).
\]
Writing this in matrix form gives the double centered version of the
data 
\begin{equation}
\overline{\overline{\boldsymbol{X}}}=\widetilde{\boldsymbol{X}}-\overline{\boldsymbol{x}}_{CO}\cdot\boldsymbol{1}_{1,n}-\boldsymbol{1}_{d,1}\cdot\overline{\boldsymbol{x}}_{RT}^{t}+\overline{x}_{AA}\cdot\boldsymbol{1}_{d,n}=\label{eqP-E:OverallMeanCentData}
\end{equation}
\[
=\left(\boldsymbol{I}_{d}-\frac{\boldsymbol{1}_{d,d}}{d}\right)\widetilde{\boldsymbol{X}}\left(\boldsymbol{I}_{n}-\frac{\boldsymbol{1}_{n,n}}{n}\right).
\]
where the \emph{grand mean\index{grand mean}} (i.e. mean over all
of the entries of the data matrix) is the scalar
\[
\overline{x}_{AA}=d^{-1}\sum_{i=1}^{d}\overline{x}_{i,A}=n^{-1}\sum_{j=1}^{n}\overline{x}_{A,j}=(nd)^{-1}\sum_{i=1}^{d}\sum_{j=1}^{n}\widetilde{x}_{i,j}=(nd)^{-1}\boldsymbol{1}_{1,d}\boldsymbol{X}\boldsymbol{1}_{n,1}.
\]
One way of understanding why the grand mean $\overline{x}_{AA}$ should
be added back in is that an overall vertical shift in the entire data
matrix (e.g. the fourth term in (\ref{eqP-B:CenteringToyDataGen}))
will be subtracted twice from the previous terms, so it should be
added back in to give the correct overall impact. 

\begin{figure}[H]
\includegraphics[scale=0.8]{ChapterP/OODAbookChpPFigAA4ToyOverallCent}

\caption{Double centered version of the Sine Wave data in Figure \ref{figP-A1:Centering-Raw}.
Shows most structure has been removed leaving just the product of
sloping lines.\label{figP-A4:Centering-Overall}}
\end{figure}
Note that in Figure \ref{figP-A4:Centering-Overall} both the horizontal
and vertical effects that dominated Figure \ref{figP-A1:Centering-Raw}
have been removed. That leaves mostly just the product of linear components
generated by $T_{3}\in\mathcal{S}_{FT}^{\perp}\cap\mathcal{S}_{FO}^{\perp}$
showing up as a distinctive pattern that results from a set of lines
whose slopes change linearly from positive to negative. Note that
in the heat-map\index{heat-map|)} in the left panel, the color bar
has been changed to reflect a finer color scale (than in Figure \ref{figP-A1:Centering-Raw}),
since use of the original scale results in very light colors that
are harder to interpret. However, the vertical axes in the center
and right panels have been kept the same, so one can still visually
gauge the relatively small magnitude of the tilted lines that are
a consequence of the $\boldsymbol{T}_{3}$ product linear mode of
variation. Further discussion of these different notions of mean centering
of a data matrix, and their impact on PCA type decompositions of variation
can be found in \citet{ZhangMarronShen2006Centering}. 

Next, this Sine Wave data set defined in (\ref{eqP-B:CenteringToyDataGen})
is used to understand how these centering operations interact with
classical functional PCAs, such as illustrated using the \index{Tilted Parabolas data}Tilted
Parabolas example in Figure \ref{fig:10d-ToyPCA}. In particular,
Figure \ref{figP-A5:Centering-ColCurvDat} shows an analysis of the
Sine Wave data, using the format of Figure \ref{fig:10d-ToyPCA}.
The top left panel shows the column object raw data curves using a
rainbow color scheme (starting with magenta for the left hand column
of the heat map through red for the right hand column, i.e. over indices
$j=1,\cdots,10$). These are colored versions of the curves in the
top center panel of Figure \ref{figP-A1:Centering-Raw}, where the
color order is best understood using the scores plot in the bottom
right panel of Figure \ref{figP-A5:Centering-ColCurvDat}. The top
center panel shows the data object mean curve, $\overline{\boldsymbol{x}}_{CO}$,
which can be thought of in an $\mathbb{R}^{d\times n}$ way as an
overlay of the $n$ identical columns (as curves) of $\overline{\boldsymbol{x}}_{CO}\boldsymbol{1}_{1,n}=P_{\mathcal{S}_{FT}}\left(\widetilde{\boldsymbol{X}}\right).$
The top right panel shows the residuals from subtracting that mean,
i.e. $P_{\mathcal{S}_{FT}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)$.
The latter are colored versions of those in the center panel of Figure
\ref{figP-A2:Centering-Col}. This makes it clear that the classical
FDA\index{functional data analysis} mean centering operation should
be thought of as removing a mode of variation from the data (in this
case the $\boldsymbol{T}_{1}$ mode). The PC1 mode of variation in
the bottom left panel of Figure \ref{figP-A5:Centering-ColCurvDat}
is essentially $\boldsymbol{T}_{2}=\boldsymbol{1}_{d,1}\boldsymbol{v}^{t}$,
i.e. a product of the flat \index{loadings|(}loadings vector $\boldsymbol{1}_{d,1}$
with parabolic scores $\boldsymbol{v}^{t}$. The latter are highlighted
in the scores plot in the bottom right, where both the rainbow colors
and the parabolic pattern are clear. An important point is that the
appearance of the flat trait mode $\boldsymbol{T}_{2}$ as the PC1
mode of variation should not be thought of as typical, and instead
is just a consequence of the design of this particular example. However,
first PC modes frequently have at least some influence\index{influence|(}
from the flat trait mode $P_{\mathcal{S}_{FO}}\left(\widetilde{\boldsymbol{X}}\right)$,
see for example the Spanish Mortality data\index{Spanish Mortality data}
in Figure \ref{fig:Mortality-PC1}. The residuals in the bottom center
plot of Figure \ref{figP-A5:Centering-ColCurvDat} are linear with
slopes again following the rainbow colors, which capture the small
variation in the slopes of the line segments in the top right panel,
thus displaying the $\boldsymbol{T}_{3}$ component of the variation
that is a colored version of the curves in the center panel of Figure
\ref{figP-A4:Centering-Overall}.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterP/OODAbookChpPFigAA5ToyColObjPCA1}

\includegraphics[scale=0.85]{ChapterP/OODAbookChpPFigAA5ToyColObjPCA2}

\caption{Functional PCA of the $n=10$ column objects in the Sine Wave data.
Shows the sine wave component is common structure appearing in the
column object mean (top center panel). PC1 mode of variation is essentially
a constant function (flat mode), as seen in the bottom left loadings
plot. Residuals from this mode show up as lines with ordered slopes
(bottom center). Rainbow color scheme follows the ordering of the
columns, clearly understood in the scores plot (lower right). \label{figP-A5:Centering-ColCurvDat}}
\end{figure}
The dual version of FDA, based on row trait vectors as data objects,
is shown in Figure \ref{figP-A6:Centering-RowCurvDat}. This time
a very different heat scale color scheme is used, running from black
through red to yellow following the indices $i=1,\cdots,20$, to indicate
that now color indicates rows of the data matrix $\widetilde{\boldsymbol{X}}$.
The top left panel shows the colored parabolic version of the curves
in the left panel of Figure \ref{figP-A1:Centering-Raw}. This time
it is the $\boldsymbol{T}_{2}$ parabolic mode of variation that shows
up as common structure captured by the row trait mean (top center).
The sine wave component $\boldsymbol{T}_{1}$ is now much harder to
see in the first PC mode in the bottom left panel, which again are
essentially flat lines. This is because the sine wave appears in the
ordering of the curves, as reflected in the scores shown in the lower
right panel. The residuals in the bottom center panel are again lines
with linearly varying slope driven by $\boldsymbol{T}_{3}$, where
the indices $i$ and $j$ have swapped roles relative to the bottom
right panel of Figure \ref{figP-A5:Centering-ColCurvDat}. 

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterP/OODAbookChpPFigAA6ToyRowObjPCA1} 

\includegraphics[scale=0.85]{ChapterP/OODAbookChpPFigAA6ToyRowObjPCA2}

\caption{Functional PCA of the $d=20$ row trait vectors in the Sine Wave data.
Heat color scheme is used to highlight row ordering. From this viewpoint
the parabola appears as mean structure, and the sine wave component
generates the variation of the flat curves in the first mode of variation.
Again the residuals are lines with ordered slopes.\label{figP-A6:Centering-RowCurvDat} }
\end{figure}
It can be challenging to connect the structure seen in the column
object PCA in Figure \ref{figP-A5:Centering-ColCurvDat} in an intuitive
fashion with the row trait PCA\index{row trait PCA}\index{PCA!row trait}
shown in Figure \ref{figP-A6:Centering-RowCurvDat}. Attention to
where the terms of (\ref{eqP-B:CenteringToyDataGen}) appear in Figures
\ref{figP-A5:Centering-ColCurvDat} and \ref{figP-A6:Centering-RowCurvDat},
makes clear the differences and also connections between these PCAs
based on column data objects and on \index{row trait!vector}row trait
vectors. The main point is that loadings in one correspond to scores
in the other. In this sense the two types of PCA are duals to each
other. That point is generally somewhat obscured by centering issues,
but can be clarified by focusing on modes of variation. This time
the PC1 mode of variation in the bottom left panel of Figure \ref{figP-A6:Centering-RowCurvDat}
is essentially $\boldsymbol{T}_{1}=\boldsymbol{u}\boldsymbol{1}_{1,n}$,
which is a product of the constant direction vector $\boldsymbol{1}_{1,n}\in\mathbb{R}^{n}$
with sine wave scores $\boldsymbol{u}$. The scores plot in the lower
right reveals the sine wave pattern following the circles from cool
colors at the bottom through hotter colors at the top. Also this time
the bi-linear PC1 residuals in the bottom center plot have a form
similar to those in the bottom center plot of Figure $\ref{figP-A4:Centering-Overall}$,
but with a more informative coloring.

A special property of this Sine Wave\index{Sine Wave data} toy data
set is that the mean mode of each is the first PC mode of the other.
Hence Figures \ref{figP-A5:Centering-ColCurvDat} and \ref{figP-A6:Centering-RowCurvDat}
are able to illustrate the important general principle that for each
\emph{common} mode of variation, the scores for one are the loadings
for the other. However, in general for column object PCA, after mean
object centering, the constant direction vector $\boldsymbol{1}_{d,1}\in\mathbb{R}^{d}$
will not point in any single PC direction, so the modes of variation
of the mean residuals will be not be the same as for row trait object
PCA. Thus exact duality will usually not hold, except for modes where
both loadings and scores are orthogonal to their respective flat directions,
i.e. sum to 0. When this type of duality is desirable (i.e. when it
is important that the scores on one analysis be the loadings of the
other), it can be guaranteed by starting the PCA at the double centered
version of the data $\overline{\overline{\boldsymbol{X}}}$ from (\ref{eqP-E:OverallMeanCentData}).
This amounts to separately accounting for both of the mean-based modes
of variation, before finding the PC modes.

Another viewpoint of centering ideas is \emph{data filtering}\index{data filtering}.
Filtering a data set by row trait mean centering $P_{\mathcal{S}_{FO}^{\perp}}\left(\widetilde{\boldsymbol{X}}\right)$
removes the flat mode of variation in (conventional) column object
PCA (sometimes quite desirable, e.g. for removing library effects
in gene expression). Similarly filtering by column object mean centering
removes the flat mode of variation in row trait (sometimes called
dual) PCA. Other common types of filtering include scaling (e.g. dividing
either rows or columns by their standard deviation) as discussed in
Section \ref{secD4:Standardization}, and transformation\index{transformation|(}
as discussed in Section \ref{secD2:Preprocess-Trans}. 

\subsection{Singular Value Decomposition\label{subsecP-A-B:SVD}}

A typical conception of PCA is sequential in nature through the components
(modes of variation), where at each step one searches for a direction
of maximal variation in the subspace orthogonal to the previously
found directions, and then uses projections of the data objects on
these directions to define the modes. For example this viewpoint was
taken in discussing the analysis of the Spanish Mortality data\index{Spanish Mortality data}
in Section \ref{sec:Curves-as-Data} and the Tilted Parabolas\index{Tilted Parabolas data}
(Twin Arches) example in Section \ref{subsecB.B.1:ExampleParabolas}
(\ref{subsecB.B.2:ExampleTwinArches} respectively). But an interesting
property of PCA is that finding the full set of modes can be usefully
considered to be a single matrix operation.

As noted for example in Section 3.5 of \citet{Jolliffe2002}, \emph{Singular
Value Decomposition} (SVD) provides a direct and simple approach to
computing PCA, and to understanding its properties. The full matrix
version of the SVD of a $d\times n$ data matrix $\widetilde{\boldsymbol{X}}$
is 
\begin{equation}
\widetilde{\boldsymbol{X}}=\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^{t},\label{eqP-E:SVD}
\end{equation}
where $\boldsymbol{U}\in O(d)$ (the set of orthogonal matrices as
defined in Section \ref{subsec:Covariance-matrix-distances}) is a
$d\times d$ orthonormal basis matrix of $\mathbb{R}^{d}$, where
$\boldsymbol{D}$ is a $d\times n$ diagonal matrix (i.e. all matrix
entries are $0$ except on the main diagonal) of \emph{\index{singular values|(}singular
values} that are nonnegative and assumed to be sorted in decreasing
order, and where $\boldsymbol{V}\in O(n)$ is an $n\times n$ orthonormal
basis matrix of $\mathbb{R}^{n}$. 

An interesting interpretation of the SVD comes from thinking of $\widetilde{\boldsymbol{X}}$
as the matrix of coefficients of a linear transformation\index{transformation!linear}
from $\mathtt{\mathbb{R}}^{n}$ to $\mathbb{R}^{d}$, when multiplying
from the left, and also of linear transformation from $\mathbb{R}^{d}$
to $\mathtt{\mathbb{R}}^{n}$, under right multiplication applied
to the transpose of a vector. In both cases, note that the transformations
corresponding to orthonormal matrices are \emph{isometries\index{isometries}.}
These are elements of $O(d)$ (or $O(n)$), which are a little more
general than rigid rotations (elements of $SO(d)$ or $SO(n)$ as
defined in Section \ref{subsecI-C-C:Procrustes}), since various axis
flippings are also included. Note that they preserve Euclidean distance\index{distance!Euclidean}s.
Furthermore linear transformation by a diagonal matrix simply rescales
(possibly zeroing out) each vector entry. Thus the SVD in (\ref{eqP-E:SVD})
represents both (i.e. left and right) linear transformations implied
by $\widetilde{\boldsymbol{X}}$ as essentially a rotation, followed
by a coordinate by coordinate rescaling (or zeroing), which is followed
by another rotation.

While the full matrix version of the SVD defined above lends itself
to linear transformation interpretation\index{SVD!linear transformation interpretation},
it is usually numerically inefficient, in the sense that it can be
written in terms of smaller matrices. When $d\neq n$ the diagonal
matrix $\boldsymbol{D}$ has at most only $d\wedge n$ (where $\wedge$
denotes the minimum) nonzero elements, so there are either rows or
columns of $\boldsymbol{D}$ that are all zero and thus are not actually
used in the product reconstruction of $\widetilde{\boldsymbol{X}}$.
In particular, replacing by appropriate sub-matrices, the SVD in (\ref{eqP-E:SVD})
can have the same form where $\boldsymbol{U}$ is $n\times(d\wedge n)$,
$\boldsymbol{D}$ is $(d\wedge n)\times(d\wedge n)$, and $\boldsymbol{V}$
is $(d\wedge n)\times n$.

Further reduction is available in the lower rank case when $\widetilde{\boldsymbol{X}}$
has rank $r<(d\wedge n)$, because then some of the singular values
are $0$, so again there are parts of the matrices that have no effect
on the product. In this case, (\ref{eqP-E:SVD}) still has the same
form but now $\boldsymbol{U}$ is $n\times r$, $\boldsymbol{D}$
is $r\times r$, and $\boldsymbol{V}$ is $r\times n$. This is the
most efficient \emph{loss-less} version of SVD which gives an exact
representation of $\widetilde{\boldsymbol{X}}$. The columns of 
\begin{equation}
\boldsymbol{U}=\left[\boldsymbol{u}_{1}\cdots\boldsymbol{u}_{r}\right],\quad\boldsymbol{V}=\left[\boldsymbol{v}_{1}\cdots\boldsymbol{v}_{r}\right]\label{eqP-I:SingularVectors}
\end{equation}
 are called the \emph{left} and \emph{right singular vectors\index{singular vectors}}
respectively. 

SVD can be thought of as providing the solution to a number of optimization
problems. Two of these are in terms of optimal projections onto direction
vectors in $\mathbb{R}^{d}$ and in $\mathbb{R}^{n}$. Recalling the
notation (\ref{eqB.C:Spheres}) for high dimensional spheres, the
set of all direction vectors in $\mathbb{R}^{d}$ is the sphere $S^{d-1}=\left\{ \boldsymbol{u}\in\mathbb{R}^{d}:\left\Vert \boldsymbol{u}\right\Vert =1\right\} $.
Note that $S^{d-1}$ is one way to represent the Grassmannian manifold\index{manifold!Grassmannian},
$\mathbf{Gr}(1,d)$, i.e. the set of all 1-d subspaces of $\mathbb{R}^{d}$.
The $1\times n$ row vector whose entries are projection coefficients
(scores) of each data object (column of $\widetilde{\boldsymbol{X}}$)
onto a given direction vector $\boldsymbol{u}\in S^{d-1}$ is given
as $\boldsymbol{u}^{t}\widetilde{\boldsymbol{X}}$. A useful measure
of \emph{signal power in the direction} $\boldsymbol{u}$ is the sum
of squared scores, $\boldsymbol{u}^{t}\widetilde{\boldsymbol{X}}\left(\boldsymbol{u}^{t}\widetilde{\boldsymbol{X}}\right)^{t}=\boldsymbol{u}^{t}\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{X}}{}^{t}\boldsymbol{u}$.
The left singular vectors $\boldsymbol{u}_{1},\cdots,\boldsymbol{u}_{d}$
give (sequential) \emph{\index{maximal signal power}maximal signal
power} of possible projections in $\mathbb{R}^{d}$ in the sense that
\begin{equation}
\boldsymbol{u}_{l}=\arg\max_{\boldsymbol{u}\in S^{d-1},\boldsymbol{u}\perp\boldsymbol{u}_{1},\cdots,\boldsymbol{u}_{l-1}}\boldsymbol{u}^{t}\widetilde{\boldsymbol{X}}\widetilde{\boldsymbol{X}}{}^{t}\boldsymbol{u}\label{eqP-J:SVDDirnOpt}
\end{equation}
for $l=1,\cdots,d$. Furthermore the maximum signal powers are the
squares of the \emph{singular values} $s_{1},\cdots,s_{d}$, the diagonal
elements of $\boldsymbol{D}$ (using the convention $s_{l}=0$ for
$l>\left(d\wedge n\right)$). In the case of the appropriately centered
and scaled data matrix in (\ref{eqP-G:ColObjCentData}) below, these
squared singular values drive the relative heights of the red circles
in the scree plot\index{scree plot}s introduced in Figure \ref{figB-Y:Screeplots}.
Thus SVD simultaneously provides the full PCA decomposition without
the need for sequential computation. 

The right singular vectors (columns of $\boldsymbol{V}$) solve parallel
optimization problems in $\mathbb{R}^{n}$, in particular, given a
direction vector $\boldsymbol{v}\in S^{n-1}\subset\mathbb{R}^{n}$,
the $d\times1$ column vector of projection coefficients (loadings)
of the rows of $\widetilde{\boldsymbol{X}}$ onto $\boldsymbol{v}$
is $\widetilde{\boldsymbol{X}}\boldsymbol{v}$. Signal power in the
direction $\boldsymbol{v}$ is the sum of squares $\left(\widetilde{\boldsymbol{X}}\boldsymbol{v}\right)^{t}\widetilde{\boldsymbol{X}}\boldsymbol{v}=\boldsymbol{v}^{t}\widetilde{\boldsymbol{X}}^{t}\widetilde{\boldsymbol{X}}\boldsymbol{v}$,
which are sequentially maximized by the right singular vectors $\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{n}$,
\[
\boldsymbol{v}_{l}=\arg\max_{\boldsymbol{v}\in S^{n-1},\boldsymbol{v}\perp\boldsymbol{v}_{1},\cdots,\boldsymbol{v}_{l-1}}\boldsymbol{v}^{t}\widetilde{\boldsymbol{X}}{}^{t}\widetilde{\boldsymbol{X}}\boldsymbol{v}
\]
for $l=1,\cdots,n$, where the maximum values are the \emph{same}
singular values\emph{\index{singular values|)}} $s_{1},\ldots,s_{n}$
using appropriate $0$s as needed to fill in discrepancies due to
$d\neq n$. 

Another important optimization problem solved by SVD (key to deep
understanding of the modes of variation generated by PCA) can be seen
by writing it as a sum of rank 1 matrices:\index{}
\begin{equation}
\widetilde{\boldsymbol{X}}=\sum_{l=1}^{r}s_{l}\boldsymbol{u}_{l}\boldsymbol{v}_{l}^{t},\label{eqP-L:SVDcomps}
\end{equation}
where $s_{l}$ is the $l$-th singular value and $\boldsymbol{u}_{l}$
and $\boldsymbol{v}_{l}$ are defined in (\ref{eqP-I:SingularVectors}).
This representation gives easy insight into the good matrix approximation
properties of SVD. In particular, for $k\leq r$ (the rank of $\widetilde{\boldsymbol{X}}$),
from the ordering of the singular values $s_{1}\geq\cdots\geq s_{r}>0$,
it follows that $\widetilde{\boldsymbol{X}}_{k}=\sum_{l=1}^{k}s_{l}\boldsymbol{u}_{l}\boldsymbol{v}_{l}^{t}$
is the \index{best rank $k$ approximation}best rank $k$ approximation
of $\widetilde{\boldsymbol{X}}$ in the sense that 
\[
\widetilde{\boldsymbol{X}}_{k}=\arg\min_{\boldsymbol{M}\in\mathcal{R}_{k}}\left\Vert \widetilde{\boldsymbol{X}}-\boldsymbol{M}\right\Vert _{F},
\]
where $\mathcal{R}_{k}$ is the set (not a subspace) of matrices of
rank $\leq k$. Note that $\widetilde{\boldsymbol{X}}_{k}$ has a
representation of the form (\ref{eqP-E:SVD}), $\widetilde{\boldsymbol{X}}_{k}=\boldsymbol{U}_{k}\boldsymbol{D}_{k}\boldsymbol{V}_{k}^{t}$
where $\boldsymbol{\boldsymbol{U}}_{k}$, $\boldsymbol{D}_{k}$, $\boldsymbol{V}_{k}^{t}$
are the first $k$ columns of $\boldsymbol{U}$, the upper $k\times k$
sub-diagonal of $\boldsymbol{D}$, and the first $k$ rows of $\boldsymbol{V}^{t}$,
respectively. This also shows how each rank $k$ SVD can be viewed
as a projection onto $\mathcal{R}_{k}$ in $\mathbb{R}^{d\times n}$,
the space of $d\times n$ matrices, as was done in Section \ref{subsecP-A-A:Data-Centering}.
In addition, each rank one matrix $s_{l}\boldsymbol{u}_{l}\boldsymbol{v}_{l}^{t}$
is the projection of the data matrix $\widetilde{\boldsymbol{X}}$
onto the one dimensional subspace generated by $\boldsymbol{u}_{l}\boldsymbol{v}_{l}^{t}$,
and the rank $k$ approximation is the projection onto the $k$ dimensional
subspace generated by $\boldsymbol{u}_{1}\boldsymbol{v}_{1}^{t},\cdots,\boldsymbol{u}_{k}\boldsymbol{v}_{k}^{t}$.
Finally note that the representation (\ref{eqP-L:SVDcomps}) connects
with the discussion of modes of variation in Section \ref{subsecBA.1.1:FormalizingModes}.
In particular, each rank 1 matrix $s_{l}\boldsymbol{u}_{l}\boldsymbol{v}_{l}^{t}$
(for $l=1,\cdots,k$) is a mode of variation, containing both the
loadings ($\boldsymbol{u}_{l}$) and scores ($s_{l}\boldsymbol{v}_{l}^{t}$).

The representation (\ref{eqP-L:SVDcomps}) also provides a very useful
connection between the SVD loadings and scores. In particular, the
scores are projections of the data onto the (subspace generated by
the) respective loadings vectors (in $\mathbb{R}^{d}$) in the sense
that for $l=1,\cdots,r$ 
\begin{equation}
\boldsymbol{u}_{l}^{t}\widetilde{\boldsymbol{X}}=\boldsymbol{u}_{l}^{t}\sum_{l'=1}^{r}s_{l'}\boldsymbol{u}_{l'}\boldsymbol{v}_{l'}^{t}=s_{l}\boldsymbol{v}_{l}^{t}.\label{eqP.U:Loadings2Scores}
\end{equation}
Similarly, loadings have a parallel simple representation in terms
of scores as 
\begin{equation}
\widetilde{\boldsymbol{X}}\boldsymbol{v}_{l}=\sum_{l'=1}^{r}s_{l'}\boldsymbol{u}_{l'}\boldsymbol{v}_{l'}^{t}\boldsymbol{v}_{l}=s_{l}\boldsymbol{u}_{l}.\label{eqP.V:Scores2Loadings}
\end{equation}
The matrix version of (\ref{eqP.U:Loadings2Scores}) is the projection
of each data object (column of $\widetilde{\boldsymbol{X}}$) onto
each of the singular vectors, through the matrix inner product calculation
\begin{equation}
\boldsymbol{U}^{t}\widetilde{\boldsymbol{X}}=\boldsymbol{U}^{t}\boldsymbol{U}\boldsymbol{D}\boldsymbol{V}^{t}=\boldsymbol{D}\boldsymbol{V}^{t}.\label{eqP-K:SVDscores}
\end{equation}
These inner products give the coefficients of the projections that
are called SVD scores, which are the basis of many useful visualizations
(e.g. scatterplot matrices\index{scores!scatterplot|(} as discussed
in Section \ref{secF.C:Scatterplot-Matrix-Views}) in this book. 

Perhaps the most direct \index{PCA!mathematical understanding}mathematical
understanding of PCA comes from viewing it as the SVD of the (sample
size scaled) column object centered version of the data, 
\begin{equation}
\breve{\boldsymbol{X}}=n^{-1/2}\left(\widetilde{\boldsymbol{X}}-\overline{\boldsymbol{x}}_{CO}\boldsymbol{1}_{1,n}\right)=n^{-1/2}\widetilde{\boldsymbol{X}}\left(\boldsymbol{I}_{n}-\frac{1}{n}\boldsymbol{1}_{n,n}\right),\label{eqP-G:ColObjCentData}
\end{equation}
written as 
\begin{equation}
\breve{\boldsymbol{X}}=\breve{\boldsymbol{\boldsymbol{U}}}\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}.\label{eqP-H:SVDofData}
\end{equation}
For $l=1,\cdots,r$, the $l$-th set of PCA \emph{scores}\index{PCA!scores},
i.e. the $l$-th \emph{principal components}\index{principal components},
are in the $l$-th row of the $r\times n$ matrix $\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}$.
Recall these scores provide insights into how the data objects relate
to each other through scatterplots as seen in Section \ref{secF.C:Scatterplot-Matrix-Views}
and other places. As noted starting in Section \ref{sec:Curves-as-Data},
insight into the drivers of these relationships comes from the \emph{loadings}\index{PCA!loadings}.
For $l=1,\cdots,r$, the $l$-th column of $\breve{\boldsymbol{U}}$
is the $l$-th loadings vector which is the direction (recall a vector
of norm $1$) in $\mathbb{R}^{d}$ of $l$-th largest variation in
the data. For $i=1,\cdots,d$, the $i$-th entry of the loadings vector
reflects direction and magnitude of the influence of the $i$-th variable
on the $l$-th direction. The PC scores and loadings are intimately
related to each other via projections. In particular, since 
\begin{equation}
\breve{\boldsymbol{U}}^{t}\breve{\boldsymbol{X}}=\breve{\boldsymbol{U}}^{t}\breve{\boldsymbol{U}}\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}=\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t},\label{eqP.AI:ScoresAsProd}
\end{equation}
the $l$-th row of the scores matrix $\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}$
is the $n\times1$ vector of inner products (essentially projection
coefficients) of the $l$-th loadings vector with the data objects
(columns of $\breve{\boldsymbol{X}}$). Similarly 
\begin{equation}
\breve{\boldsymbol{X}}\left(\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}\right)^{t}\breve{\boldsymbol{D}}^{-2}=\breve{\boldsymbol{U}}\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}\breve{\boldsymbol{V}}\breve{\boldsymbol{D}}\breve{\boldsymbol{D}}^{-2}=\breve{\boldsymbol{U}}\label{eqP.AJ:LoadingsAsProd}
\end{equation}
 shows that the loadings can be represented as a normalization\index{normalization}
(rescaling by the inverse variances) of the inner products of the
centered data matrix and the scores. In summary, (\ref{eqP.AI:ScoresAsProd})
and (\ref{eqP.AJ:LoadingsAsProd}) show that both loadings and scores
can be derived from the other by an appropriate product with the scaled,
centered data matrix.

The column centering operation in (\ref{eqP-G:ColObjCentData}) provides
enhanced interpretability in several ways, including connection with
the sequential view of PCA. First, using $\breve{\boldsymbol{X}}$
instead of $\widetilde{\boldsymbol{X}}$ in (\ref{eqP-J:SVDDirnOpt})
gives the interpretation of PC directions first discussed in Section
\ref{sec:Curves-as-Data}: $\breve{\boldsymbol{u}}_{1}$ is the direction
vector (based at the sample mean) that \emph{maximizes sample variance}
of the projections. Note that it is the centering operation and division
by $n^{1/2}$ that changes the criterion from sum of squares to sample
variance. Since the columns of $\breve{\boldsymbol{U}}$ are orthonormal,
each $\breve{\boldsymbol{u}}_{l}$ (for $l=2,\cdots,r$) is the direction
vector in the subspace orthogonal to $\breve{\boldsymbol{u}}_{1},\cdots,\breve{\boldsymbol{u}}_{l-1}$
that maximizes the sample variance of the data projections. A second
enhancement resulting from column object centering was illustrated
in the bottom panels of Figure \ref{figP-B:SVDvsPCA}, comparing uncentered
SVD with PCA in a toy example. In particular, the column object mean
centering used in PCA (together with the orthogonality of the rows
of $\breve{\boldsymbol{V}}^{t}$, i.e. the columns of $\breve{\boldsymbol{V}}$)
results in more easily interpretable \emph{uncorrelated} scores scatterplots.
As noted in Section \ref{secF.E:ComboView-Centering}, PCA scores
scatterplots all have correlation $0$ because otherwise the direction
vector could be rotated to result in a larger sample variance of the
scores. Another view of this phenomenon comes from calculating the
$r\times r$ matrix of inner products of the scores vectors, again
using the orthonormality of the columns of $\breve{\boldsymbol{V}}$
and the fact that $\breve{\boldsymbol{D}}$ is diagonal, 
\begin{equation}
\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}\left(\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}\right)^{t}=\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}\breve{\boldsymbol{V}}\breve{\boldsymbol{D}}^{t}=\breve{\boldsymbol{D}}^{2}.\label{eqP.AE:OrthoScores}
\end{equation}
Orthogonality of the scores vectors (recall rows of $\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}$),
and hence $0$ sample correlation of the entries of those vectors
as seen in scores scatterplots, follows from the fact that $\breve{\boldsymbol{D}}^{2}$
is also a diagonal matrix, whose off diagonal entries are thus all
$0$. 

Without column object mean centering, there is generally some (i.e.
non-zero) sample correlation of SVD scores plots, again as demonstrated
in the lower left panel of Figure \ref{figP-B:SVDvsPCA}. However,
in many situations, this correlation is not distracting, because frequently
one of the singular vectors (usually the first) happens to point approximately
in the direction of the sample column object mean. In particular,
many data sets tend to have an important flat mode of variation, for
example $\boldsymbol{T}_{2}$ as illustrated in (\ref{eqP-B:CenteringToyDataGen}).
That flat mode is often roughly in the first SVD direction, so the
SVD modes of variation are frequently similar to those of a standard
PCA.

Again a major lesson of Figure \ref{figP-B:SVDvsPCA} is that column
mean object centering results in uncorrelated PC scores. The symmetry
of the SVD in terms of rows and columns, readily apparent from (\ref{eqP-E:SVD})
and (\ref{eqP-H:SVDofData}), suggests a parallel relationship: row
trait mean centering implies that such pairwise scatterplots of the
PCA loadings will be similarly uncorrelated. While column object mean
centering is routine for PCA, row trait mean centering (or double
mean centering as defined in (\ref{eqP-E:OverallMeanCentData})) is
not. As shown in Figure \ref{figP-B:SVDvsPCA} appropriate centering
will ensure that scores are uncorrelated. This does not seem to have
been widely noticed for two reasons. First, with the noticeable exception
of \emph{biplot} data views, as defined by \citet{Gabriel1971Biplots},
and more completely developed in \citet{GowerHand1996Biplots,Gower1966DS},
there has not been strong interest in visualizing PC loadings. While
biplots can be insightful in low dimensional contexts, a version adapted
to high dimensions has not yet been developed. The second reason is
that in typical analyses, the column mean object centering removes
the grand mean $\overline{x}_{AA}$, which often is a large factor
in the lack of correlation in the scores. An exception is shown in
Figure \ref{figP-C:SpanMortLoadCorr}, based on the Spanish Mortality
data\index{Spanish Mortality data} studied in Section \ref{sec:Curves-as-Data}.

\begin{figure}[H]
\includegraphics[scale=0.85]{ChapterP/OODAbookChpPFigCMortLoadCorr}

\caption{PCA loadings for Spanish Mortality data studied in Section \ref{sec:Curves-as-Data}.
Heat coded coloring indicates age from 0 (black) to 98 (yellow). Shows
correlations in loadings caused by the lack of row trait\index{row trait|)}
mean centering. \label{figP-C:SpanMortLoadCorr}}
\end{figure}
Note that the correlations between the first and second as well as
the second and third loadings vectors are small, but the first and
third are strongly negatively correlated. Furthermore the first loadings
do not have mean $0$. Thus, in situations where it is important to
simultaneously consider both scores and loadings, such as the multi-block
analyses studied in Section \ref{secP-B:TwoBlockDecomp}, it can be
important to use an SVD of the double centered version of the data
$\overline{\overline{\boldsymbol{X}}}$ defined in (\ref{eqP-E:OverallMeanCentData}).

Also note that there are strong connections with the column object
oriented analysis in Section \ref{sec:Curves-as-Data}. In particular
the pattern in the first loadings can be seen in the top center panel
by tracing (through the ages on the horizontal axis) of the upper
magenta curve in the left panel of Figure \ref{fig:Mortality-PC1}.
There is an initial relatively small value (black) with an initial
jump followed by slow decrease with later wiggles due to the age rounding
(red and yellow). In that same panel of Figure \ref{figP-C:SpanMortLoadCorr}
the pattern in the second loadings (left panel of Figure \ref{fig:Mortality-PC2})
also follows the visible magenta curve, initially negative swinging
to positive (for the 20 to 45 year olds) and then back to negative.
This is natural because as noted in Section \ref{subsecP-A-A:Data-Centering}
the usual PCA (with columns as data objects) scores and loadings defined
above essentially (modulo centering\index{centering|)}) become the
loadings and scores respectively when the rows are taken as data objects,
which is now easy to understand from taking the transpose of (\ref{eqP-H:SVDofData}).

It is also insightful to formulate these quantities in terms of the
sample covariance matrix $\widehat{\boldsymbol{\varSigma}}$ defined
at (\ref{eqB.F:SampleCovariance}). Given a data matrix $\widetilde{\boldsymbol{X}}$,
the entries of $\widehat{\boldsymbol{\varSigma}}$ are the sample
variances 
\begin{equation}
\widehat{var}_{i}=\frac{1}{n}\sum_{j=1}^{n}\left(\widetilde{x}_{i,j}-\overline{x}_{i,A}\right)^{2}\label{eqP.N:SampleVar}
\end{equation}
and the sample covariances 
\begin{equation}
\widehat{cov}_{i,i'}=\frac{1}{n}\sum_{j=1}^{n}\left(\widetilde{x}_{i,j}-\overline{x}_{i,A}\right)\left(\widetilde{x}_{i',j}-\overline{x}_{i',A}\right)\label{eqP.M:SampleCov}
\end{equation}
 for $i,i'=1,\cdots,d$. The sample covariance matrix is related to
PCA through representing it as the outer product 
\begin{equation}
\widehat{\boldsymbol{\varSigma}}=\breve{\boldsymbol{X}}\breve{\boldsymbol{X}}^{t}.\label{eqP.W:CovMatOuterProd}
\end{equation}
Using the singular value representation (\ref{eqP-H:SVDofData}) gives
\begin{equation}
\widehat{\boldsymbol{\varSigma}}=\breve{\boldsymbol{\boldsymbol{U}}}\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}\breve{\boldsymbol{V}}\breve{\boldsymbol{D}}^{t}\breve{\boldsymbol{\boldsymbol{U}}}^{t}=\breve{\boldsymbol{\boldsymbol{U}}}\breve{\boldsymbol{\varLambda}}\breve{\boldsymbol{\boldsymbol{U}}}^{t}\label{eqP-O:EigenAnalofCov}
\end{equation}
where $\breve{\boldsymbol{\varLambda}}=\breve{\boldsymbol{D}}\breve{\boldsymbol{D}}^{t}$is
a diagonal matrix whose entries are the squares of the entries of
$\breve{\boldsymbol{D}}$, i.e. the squares of the singular values
$\breve{s}_{l}$. The representation (\ref{eqP-O:EigenAnalofCov})
is called an \emph{eigenvalue decomposition} or \emph{eigen analysis\index{eigen analysis}}
at many points above. This reveals that in addition to the (\ref{eqP-H:SVDofData})
representation of PCA as a Singular Value Decomposition of the column
object centered data matrix, it can also be thought of as an eigenvalue
analysis of the sample covariance matrix. The entries of $\breve{\boldsymbol{\varLambda}}$
are the eigenvalues, which divided by their sum (i.e. the total energy
in the centered data matrix $\breve{\boldsymbol{X}}$) are the heights
of the scree plot\index{scree plot}s shown e.g. in Figure \ref{figB-Y:Screeplots}.
Such eigenvalues were also fundamental to the \emph{random matrix
theory} studied in Section \ref{secM-A:RandomMatrix}.

\subsection{Gaussian Likelihood View\label{subsecP-A-C:Gauss-Likelihood}}

A common theme throughout this book has been a \emph{nonparametric},
or data centric, development of analytical methods for OODA. This
is a contrast to many introductions to multivariate analysis where
the focus is instead on underlying parametric statistical models and
resulting likelihood approaches. Low rank \emph{latent factor} models
provide a traditional approach of this latter type to PCA. In particular,
assume the $d\times n$ data matrix can be written in the form 
\[
\widetilde{\boldsymbol{X}}=\boldsymbol{\mu}\boldsymbol{1}_{1,n}+\boldsymbol{L}\boldsymbol{S}+\widetilde{\boldsymbol{E}},
\]
where for some low rank $r$, $\boldsymbol{L}$ is a $d\times r$
orthonormal loadings matrix, $\boldsymbol{S}$ is an $r\times n$
scores matrix whose rows are orthogonal and sum to $0$ (for identifiability)
and $\widetilde{\boldsymbol{E}}$ is a $d\times n$ matrix of errors
whose entries are often assumed to be independent $N\left(0,\sigma^{2}\right)$.
Rows of $\boldsymbol{S}$ are often thought of as unobserved latent
factors, and their estimation is the primary goal of \emph{factor
analysis\index{factor analysis}}. Standard calculations show that
for a given $r$, using notation from (\ref{eqF1:ColumnObjectMean})
and (\ref{eqP-H:SVDofData}), the maximum likelihood estimates of
$\boldsymbol{\mu}$, $\boldsymbol{L}$ and $\boldsymbol{S}$ are $\overline{\boldsymbol{x}}_{CO}$,
$\breve{\boldsymbol{\boldsymbol{U}}}$ and $\breve{\boldsymbol{D}}\breve{\boldsymbol{V}}^{t}$
respectively. This is a sense in which standard PCA can be thought
of from the likelihood viewpoint. 

While this approach is mathematically elegant, as noted in Section
\ref{secB2:Data-Visualization} it has the serious downside that it
has fed a common misconception that PCA is only useful when the data
approximately follow a multivariate Gaussian distribution. As seen
at many points in this book, PCA provides insightful data visualization\index{visualization|)}s
in many non-Gaussian situations. Important real data examples of this
appear in Figures \ref{fig:LungCancerPCABrushed} (Lung Cancer\index{Lung Cancer}
data) and \ref{figB-Q:PanCanPCA} (Pan Cancer). 

As noted at the beginning of this chapter, properly defined factor
analysis is closely related to PCA (as opposed to being exactly the
same). The difference is that factor analysis includes estimation
of the residual variance $\sigma^{2}$ as part of the likelihood calculation.
This results in the same loadings (column vectors of $\breve{\boldsymbol{\boldsymbol{U}}}$),
i.e. directions for visualizations, as PCA. However the scores are
somewhat impacted by the simultaneous estimation of $\sigma^{2}$. 

\subsection{PCA Computational Issues\label{subsecP-A-D:PCA-Comput}}

Good discussion of various approaches to computing PCA can be found
in Section A.1 of \citet{Jolliffe2002}. Using most software packages,
which contain well-optimized SVD functions, the SVD approach to computation
described in Section \ref{subsecP-A-B:SVD} is both insightful and
generally computationally useful.

As noted at the beginning of this chapter, PCA can give an efficient
(exact) low dimensional representation of a $d\times n$ data matrix.
In particular, it follows from (\ref{eqP-G:ColObjCentData}) and (\ref{eqP-H:SVDofData})
that 
\begin{equation}
\widetilde{\boldsymbol{X}}=\overline{\boldsymbol{x}}_{CO}\boldsymbol{1}_{1,n}+\breve{\boldsymbol{\boldsymbol{U}}}_{r}\left(n^{1/2}\breve{\boldsymbol{D}}_{r}\breve{\boldsymbol{V}}_{r}^{t}\right),\label{eqP.P:PCADataRep}
\end{equation}
where $r$ is the rank of $\widetilde{\boldsymbol{X}}-\overline{\boldsymbol{x}}_{CO}\boldsymbol{1}_{1,n}$
(which is usually the rank of $\widetilde{\boldsymbol{X}}$ when $n>d$
and one smaller otherwise), where the \emph{loadings} matrix $\breve{\boldsymbol{\boldsymbol{U}}}_{r}$
contains the first $r$ columns of $\breve{\boldsymbol{\boldsymbol{U}}}$
(these are the PCA directions of projection in $\mathbb{R}^{d}$)
and where $\left(n^{1/2}\breve{\boldsymbol{D}}_{r}\breve{\boldsymbol{V}}_{r}^{t}\right)$
is the $r\times n$ matrix of PCA scores (projection coefficients),
plotted at many points in this book starting with Figure \ref{fig:Mortality-PC1}.
Thus the scores are a translated and orthogonal transformed (rotation
and/or sign flip) version of the data that are a lower dimensional
isometric\index{isometric} (hence loss-less) representation in $\mathbb{R}^{r}$.

The data representation\index{data representation|)} (\ref{eqP.P:PCADataRep})
is closely related to the formula for simulating a $d\times n$ random
matrix $\widetilde{\boldsymbol{X}}\sim N_{d}\left(\boldsymbol{\mu},\boldsymbol{\varSigma}\right)$
as 
\begin{equation}
\widetilde{\boldsymbol{X}}=\boldsymbol{\mu}\boldsymbol{1}_{1,n}+\boldsymbol{\varSigma}^{1/2}\boldsymbol{Z}=\boldsymbol{\mu}\boldsymbol{1}_{1,n}+\boldsymbol{U}\boldsymbol{\Lambda}^{1/2}\boldsymbol{U}^{t}\boldsymbol{Z}=\boldsymbol{\mu}\boldsymbol{1}_{1,n}+\boldsymbol{U}\boldsymbol{\Lambda}^{1/2}\boldsymbol{Z}^{*}\label{eqP.T:SimulateGaussians}
\end{equation}
where both $\boldsymbol{Z}$ and $\boldsymbol{Z}^{*}$ are $d\times n$
matrices of independent standard normals since $\boldsymbol{U}^{t}$
is an orthogonal matrix. To simulate Gaussians that follow the data
set in any matrix $\widetilde{\boldsymbol{X}}$ as well as possible,
use (\ref{eqP.T:SimulateGaussians}) with $\boldsymbol{\mu}$ (and
$\boldsymbol{U}\boldsymbol{\Lambda}^{1/2}$) replaced by the sample
version $\overline{\boldsymbol{x}}_{CO}$ (and $n^{1/2}\breve{\boldsymbol{\boldsymbol{U}}}_{r}\breve{\boldsymbol{D}}_{r}$
respectively).

\section{Two Block Decompositions\label{secP-B:TwoBlockDecomp}}

As noted at many points above, a set of data objects that is usefully
viewed as columns of a matrix is a common format (recall the OODA
convention of cases as columns of the data matrix and traits, i.e.
variables as rows). However, an increasing number of OODA contexts
have been appearing where data objects naturally consist of several
such column vectors of traits that are grouped in meaningful ways.
A canonical example is the TCGA data set, introduced in Sections \ref{subsecB.2.A:LungCancerData}
and \ref{subsecB.2.B:PanCancerExample}. In addition to the gene expression
studied those sections, around a dozen other types of measurements
were made for each patient, including protein expression and mutation
status. Most measurement types result in long vectors of numbers,
so the resulting data objects are naturally considered to be sets
of vectors. For clear study of not just the measurements but also
the relationships between the groups, it is convenient to group the
vectors into several matrices, called data \emph{blocks} here. In
particular, for $b=1,\cdots,B$ let the $d_{b}\times n$ data block
be 
\begin{equation}
\widetilde{\boldsymbol{X}}^{(b)}=\left[\widetilde{\boldsymbol{x}}_{1}^{(b)}\cdots\widetilde{\boldsymbol{x}}_{n}^{(b)}\right].\label{eqP.R:Multi-BlockNotation}
\end{equation}
An important aspect of multi-block data is that for each $j=1,\cdots,n$
the corresponding measurement vectors $\widetilde{\boldsymbol{x}}_{j}^{(1)},\cdots,\widetilde{\boldsymbol{x}}_{j}^{(B)}$
are all made on the same experimental unit, e.g. tissue sample in
the case of the TCGA data first discussed in Section \ref{subsecB.2.A:LungCancerData}.
There are several synonyms for ``multi-block'', such as \emph{multi-view}
which is common in the machine learning\index{machine learning} literature.
An accessible overview of multi-block methods can be found in \citet{DeBie2005Eigenproblems}.

Time honored approaches to multi-block data, also called \index{data integration}data
integration, are \emph{Partial Least Squares\index{Partial Least Squares|(}}
(PLS) studied in Section \ref{subsecP-B-A:PLS} and \emph{Canonical
Correlation Analysis\index{Canonical Correlation Analysis}} (CCA)
discussed in Section \ref{subsecP-B-B:CCA}. As discussed there, these
powerful methods can be approached from several viewpoints, and have
a number of variants. Some of these variants are usefully understood
and contrasted through their corresponding modes of variation.

A simple starting point for this is the case $B=2$ and the joint
covariance matrix 
\begin{equation}
\left[\begin{array}{cc}
\widehat{\boldsymbol{\varSigma}}^{(1)} & \widehat{\boldsymbol{\varSigma}}^{(1,2)}\\
\widehat{\boldsymbol{\varSigma}}^{(1,2)} & \widehat{\boldsymbol{\varSigma}}^{(2)}
\end{array}\right],\label{eqP.S:JointCovMat}
\end{equation}
where $\widehat{\boldsymbol{\varSigma}}^{(1)}$ and $\widehat{\boldsymbol{\varSigma}}^{(2)}$
are the respective within block sample covariance matrices (defined
at (\ref{eqB.F:SampleCovariance}) and usefully represented as an
outer product in (\ref{eqP.W:CovMatOuterProd})) of the two blocks,
and where $\widehat{\boldsymbol{\varSigma}}^{(1,2)}$ is the \emph{\index{cross covariance matrix}cross
covariance matrix} whose entries are sample covariances of the form
(\ref{eqP.M:SampleCov}) with $i=1,\cdots,d_{1}$ (and $i'=1,\cdots,d_{2}$)
representing the index of a trait (variable) in $\widetilde{\boldsymbol{X}}^{(1)}$
(and $\widetilde{\boldsymbol{X}}^{(2)}$ respectively). In the spirit
of (\ref{eqP.W:CovMatOuterProd}), using block indexed versions of
the column object centered data matrices (\ref{eqP-G:ColObjCentData}),
the $d_{1}\times d_{2}$ cross covariance matrix has an outer product
representation as $\widehat{\boldsymbol{\varSigma}}^{(1,2)}=\breve{\boldsymbol{X}}^{(1)}\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}$,
where $\breve{\boldsymbol{X}}^{(b)}$ is the column object centered
version of the data in block $b$ for $b=1,2$. 

Additional extensions of PLS and CCA ideas include kernel (in the
classification sense of Section \ref{secJ-B:KernelMethods}) versions
in \citet{BachJordan2002KernelICA}, the sparsity approaches of \citet{Gao2015SparseCCA}
and the deep learning implementations proposed by \citet{Andrew2013DeepCCA}.
For more than two blocks there are dozens of extensions as seen in
\citet{Tenenhaus2011RegularizedCCA}. For example, CCA alone has 20
generalizations to multiple blocks: \citet{Kettenring1971CCA,Nielsen2002MultisetCCA,Asendorf2015BeyondCCA}.

\subsection{Partial Least Squares\label{subsecP-B-A:PLS}}

The terminology Partial Least Squares (PLS) has been used for a large
number of methods for multi-block analysis. The most important of
these are a set of methods for deriving joint modes of variation\index{modes of variation!PLS}
from the cross covariance matrix $\widehat{\boldsymbol{\varSigma}}^{(1,2)}$.
Good overview appears in the survey paper \citet{Wegelin2000SurveyPLS}.
The name PLS was coined by \citet{Wold1975PLS,Wold1985EncyclopediaPLS}.

A straightforward version, labeled PLS-SVD by \citet{Wegelin2000SurveyPLS},
develops joint modes of variation using ideas parallel to PCA (as
described in Section \ref{subsecP-A-B:SVD}). It starts with the SVD
of the cross covariance matrix 
\begin{equation}
\widehat{\boldsymbol{\varSigma}}^{(1,2)}=\widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{D}}\widetilde{\boldsymbol{V}}^{t}.\label{eqP.Z:PLS-SVD}
\end{equation}
For $l=1,\cdots,d_{1}\wedge d_{2}$ the columns $\widetilde{\boldsymbol{u}}_{l}$
of $\widetilde{\boldsymbol{U}}$ and $\widetilde{\boldsymbol{v}}_{l}$
of $\widetilde{\boldsymbol{V}}$ are direction vectors in $\mathbb{R}^{d_{1}}$
and $\mathbb{R}^{d_{2}}$ (respectively) that sequentially maximize
the covariance of projections of $\widetilde{\boldsymbol{X}}^{(1)}$
onto $\widetilde{\boldsymbol{u}}_{l}$ and of $\widetilde{\boldsymbol{X}}^{(2)}$
onto $\widetilde{\boldsymbol{v}}_{l}$, subject to orthogonality of
columns of both $\widetilde{\boldsymbol{U}}$ and $\widetilde{\boldsymbol{V}}$.
Each pair of direction vectors generates \emph{joint modes of variation\index{modes of variation!joint|(}}
in $\mathbb{R}^{d_{1}\times n}$ and $\mathbb{R}^{d_{2}\times n}$
by computing corresponding scores vectors using the data matrix multiplication
principle of (\ref{eqP.AI:ScoresAsProd}). The resulting \emph{cross
covariance scores} are projection coefficients of the data objects
onto the direction loading vectors as in (\ref{eqP.U:Loadings2Scores}).
Hence these are computed as $\left(d_{1}\wedge d_{2}\right)\times n$
\emph{PLS scores\index{PLS scores} }matrices 
\[
\widetilde{\boldsymbol{S}}^{(1)}=\widetilde{\boldsymbol{U}}^{t}\left(\widetilde{\boldsymbol{X}}^{(1)}-\overline{\boldsymbol{x}}_{CO}^{(1)}\boldsymbol{1}_{1,n}\right),\;\widetilde{\boldsymbol{S}}^{(2)}=\widetilde{\boldsymbol{U}}^{t}\left(\widetilde{\boldsymbol{X}}^{(2)}-\overline{\boldsymbol{x}}_{CO}^{(2)}\boldsymbol{1}_{1,n}\right).
\]
Scatterplot matrix plots, as in e.g. Figure \ref{fig:50d-ToyPCAScatMat}
(Twin Arches data\index{Twin Arches data}), of scores (projections
on the loadings vectors) within each block gives insights (parallel
to PCA \index{scatterplot matrix}scatterplot matrix plots) into object
relationships that are maximally related to the other block. Also
of interest is to study the amount of dependence between the blocks
through, for $l=1,\cdots,d_{1}\wedge d_{2}$, viewing the scatterplots
of the scores $\widetilde{\boldsymbol{S}}_{l}^{(1)}$  vs. $\widetilde{\boldsymbol{S}}_{l}^{(2)}$.
As for PCA, the entries of the vectors $\widetilde{\boldsymbol{u}}_{1}$
and $\widetilde{\boldsymbol{v}}_{1}$ are both called \emph{\index{PLS loadings}loadings}
in this context. These loadings can be plotted as in Figure \ref{figB-X:PanCanDWDLoadings}
to give insights as to the drivers of the \emph{joint variation} between
data blocks. When the loadings are interpretable as smooth curves,
such as for the Spanish Mortality data\index{Spanish Mortality data}
studied in Section \ref{sec:Curves-as-Data}, useful insights come
from multiplying the loadings vectors by scores to produce mode of
variation plots\index{modes of variation!plot}, such as those in
Figures \ref{fig:Mortality-PC1} and \ref{fig:Mortality-PC2}, and
in many other places.

While the PLS-SVD approach is quite straightforward to compute and
gives practically useful modes of variation, as demonstrated in \citet{Singh2010PLS}
and other places noted in \citet{Wegelin2000SurveyPLS}, it does have
some limitations. In particular these modes of variation do not share
the \emph{bi-orthogonality\index{bi-orthogonal|(}} property of PCA
modes, which was that \emph{both} the columns of the matrix $\breve{\boldsymbol{\boldsymbol{U}}}$
and also the columns of $\breve{\boldsymbol{V}}$ (i.e. the rows of
$\breve{\boldsymbol{V}}^{t}$) are orthonormal bases (of $\mathbb{R}^{d}$
and $\mathbb{R}^{n}$ respectively). In particular, the rows of the
block 1 scores matrix $\widetilde{\boldsymbol{S}}^{(1)}$ are generally
\emph{not} orthogonal to each other, and similarly for the rows of
$\widetilde{\boldsymbol{S}}^{(2)}$.Because each vector of scores
has sample mean 0 of its entries, these \index{non-orthogonal axes}non-orthogonalities
entail that the sample correlations between vector entries are not
0, i.e. there is some correlation between these vectors, which can
be quite apparent and even disconcerting in scores scatterplot matrix
views. However, there is some orthogonality \emph{across} the blocks,
that follows from a calculation similar to (\ref{eqP.AE:OrthoScores}).
In particular the parallel $\left(d_{1}\wedge d_{2}\right)\times\left(d_{1}\wedge d_{2}\right)$
inner product matrix is 
\[
\widetilde{\boldsymbol{S}}^{(1)}\left(\widetilde{\boldsymbol{S}}^{(2)}\right)^{t}=\widetilde{\boldsymbol{U}}^{t}\left(\widetilde{\boldsymbol{X}}^{(1)}-\overline{\boldsymbol{x}}_{CO}^{(1)}\boldsymbol{1}_{1,n}\right)\left(\widetilde{\boldsymbol{V}}^{t}\left(\widetilde{\boldsymbol{X}}^{(2)}-\overline{\boldsymbol{x}}_{CO}^{(2)}\boldsymbol{1}_{1,n}\right)\right)^{t}=
\]
\[
=\widetilde{\boldsymbol{U}}^{t}\left(\widetilde{\boldsymbol{X}}^{(1)}-\overline{\boldsymbol{x}}_{CO}^{(1)}\boldsymbol{1}_{1,n}\right)\left(\widetilde{\boldsymbol{X}}^{(2)}-\overline{\boldsymbol{x}}_{CO}^{(2)}\boldsymbol{1}_{1,n}\right)^{t}\widetilde{\boldsymbol{V}}=
\]
\begin{equation}
=\widetilde{\boldsymbol{U}}^{t}\widehat{\boldsymbol{\varSigma}}^{(1,2)}\widetilde{\boldsymbol{V}}=\widetilde{\boldsymbol{U}}^{t}\left(\widetilde{\boldsymbol{U}}\widetilde{\boldsymbol{D}}\widetilde{\boldsymbol{V}}^{t}\right)\widetilde{\boldsymbol{V}}=\widetilde{\boldsymbol{D}}.\label{eqPF:PLScrossBlockOrtho}
\end{equation}
Hence the diagonality of the singular value matrix $\widetilde{\boldsymbol{D}}$
shows that the scores in each block are orthogonal to the non-corresponding
scores in the other block. This does not seem to be widely acknowledged,
perhaps because there does not appear to be motivation to construct
such cross block scores scatterplots\index{scores!scatterplot|)}.

The above (within blocks) lack of orthogonality of the PLS-SVD scores
can be viewed as motivating a number of variants. The best known of
these comes from \citet{Wold1975PLS,Wold1985EncyclopediaPLS} who
proposed \emph{\index{PLS regression}PLS regression} as an improvement
of \emph{PCA regression}\index{PCA regression}. The latter is useful
in multiple linear regression contexts with many and/or collinear
predictors, which tend to cause numerical instability in classical
least squares algorithms. The main idea of PCA regression is to replace
those predictors with a reduced set of PCA scores, which often provide
a good low rank approximation (thus containing most of the relevant
information of the original predictors while hopefully eliminating
noise), sometimes called \emph{\index{latent variables}latent variables}.
That also avoids numerical instability in the regression formulas
because of the orthogonality of the PCA scores vectors. An obvious
drawback is that PCA is driven completely by the variation in the
predictors, which has no explicit connection to the responses. The
intuitively attractive improvement of PLS regression is to replace
the low rank approximation based on SVD of (\ref{eqP-H:SVDofData})
inherent to PCA with a more directly relevant approximation based
on SVD of (\ref{eqP.S:JointCovMat}). That gives modes of variation
that focus on the connection between predictors and responses, instead
of merely maximizing variation among the predictors. However, as discussed
above, the scores from the SVD of (\ref{eqP.S:JointCovMat}) are not
entirely suitable as predictors because of their non-orthogonality.
This can be thought of as the motivation for Wold's PLS which does
have uncorrelated scores. Those give both enhanced numerical performance
of the regression algorithm, and also more interpretable scores in
latent variable terms. 

There are various approaches to Wold's version of PLS. An algorithm-centered
approach (focused on sequentially finding appropriate directions)
can be found in \citet{Wegelin2000SurveyPLS}. The idea is to find
a first mode of variation using (\ref{eqP.Z:PLS-SVD}) as above. Additional
modes are found by iteratively solving (\ref{eqP.Z:PLS-SVD}) with
the added constraint of orthogonality of the scores vectors. Many
authors prefer a model-based latent factor way of thinking along the
lines of that approach to PCA described in Section \ref{subsecP-A-C:Gauss-Likelihood},
with the same type of sequential implementation. Other approaches
are more optimization based, see for example \citet{Rosipal2005PLS}
and \citet{Xu2013MultiViewSurvey}. 

The orthogonality of the scores in Wold's PLS comes at a price of
loss of orthogonality of the loadings vectors. The orthogonality of
scores across blocks shown in (\ref{eqPF:PLScrossBlockOrtho}) is
similarly lost in Wold's PLS. There do not seem to be two block decompositions
into modes of variation that share the bi-orthogonality property of
PCA.

A question that does not appear to be addressed in the literature
is how far are the scores resulting from PLS-SVD from being orthogonal,
i.e. how far are the PLS-SVD modes of variation from being bi-orthogonal?
This is investigated in Table \ref{tabP.B:OrthoPLSscores} which summarizes
results from a small  simulation study based on $1000$ realizations
of pairs of matrices with i.i.d. Gaussian entries of dimensions $d\times n$.
Each pair was double centered as described in Section \ref{subsecP-A-A:Data-Centering}
and scores from an SVD of the cross covariance matrix $\widehat{\boldsymbol{\varSigma}}^{(1,2)}$
were computed. The angles between the first and second scores vectors
for both matrices were computed, and the minimums and averages have
been summarized in Table \ref{tabP.B:OrthoPLSscores}. 

\begin{table}[H]
\begin{tabular}{|c|c|c|c|}
\hline 
$d$ & $n$ & Min. Angle & Max. Angle\tabularnewline
\hline 
\hline 
10000  & 10 & 88.70 & 89.67\tabularnewline
\hline 
1000 & 100 & 86.59 & 89.19\tabularnewline
\hline 
100 & 1000 & 84.65 & 88.83\tabularnewline
\hline 
10 & 10000 & 87.85 & 89.55\tabularnewline
\hline 
\end{tabular}

\caption{Summary of simulation study, exploring angles between PLS-SVD scores
vectors, in the case of i.i.d. Gaussian data. Shows scores vectors
tend to be orthogonal across a wide range of settings. \label{tabP.B:OrthoPLSscores}}
\end{table}
Note that broadly over this wide range of contexts, these angles are
all perhaps surprisingly close to orthogonal, i.e. $90^{\circ}$.
In the HDLSS\index{High Dimension Low Sample Size} case ($d\gg n$)
this is not surprising, given the fundamental ideas from Section \ref{secM-B:HDLSS}
that random directions tend to be orthogonal in those cases. Less
obvious is why the angles are also very close to $90^{\circ}$ in
the other cases. An open theoretical problem seems to be an asymptotic
analysis\index{asymptotic analysis} that quantifies these observations.
Of course, the results of Table \ref{tabP.B:OrthoPLSscores} only
apply to pure noise contexts, and some types of low rank signal can
easily lead to much smaller angles between scores vectors. 

An important limitation of PLS is that because it is driven by covariance
it does not always focus completely on joint structure but instead
also feels variation within each block, as demonstrated in Figure
\ref{figP.D:PCA-PLS-CCA}. For example when one variable has a much
larger variance than the others, it will tend to unduly influence
PLS. A solution to this problem is given in Section \ref{subsecP-B-B:CCA}. 

\subsection{Canonical Correlations\label{subsecP-B-B:CCA}}

As noted in Section \ref{subsecP-B-A:PLS}, PLS-SVD provides modes
of joint variation that maximize covariance between scores. As discussed
in the context of correlation PCA in Section \ref{secD4:Standardization},
a potential drawback of covariance is that it feels the units of the
data. This can create challenges in several contexts, an important
example being data blocks whose trait\index{trait|)}s are measured
in non-commensurate\index{commensurate|(} units. In those situations
it is very sensible to replace the covariance criterion with the unit
free Pearson's correlation coefficient. That goal results in CCA\index{Canonical Correlation Analysis|(},
which was named by \citet{Hotelling1936CCA}, who as noted at the
beginning of this chapter also coined the term PCA. 

CCA starts with directions $\widetilde{\widetilde{\boldsymbol{u}}}\in\mathbb{R}^{d_{1}}$
and $\widetilde{\widetilde{\boldsymbol{v}}}\in\mathbb{R}^{d_{2}}$.
Corresponding vectors of projections of the ($\frac{1}{n}$ rescaled)
column object centered data (\ref{eqP-G:ColObjCentData}) are $\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}}$,
$\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}\in\mathbb{R}^{n}$.
Their sample correlation is 
\[
\widehat{\rho}\left(\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}},\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}\right)=\frac{\widehat{cov}\left(\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}},\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}\right)}{\left(\widehat{var}\left(\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}}\right)\widehat{var}\left(\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}\right)\right)^{1/2}}=
\]
\[
=\frac{\left(\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}}\right)^{t}\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}}{\left(\left(\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}}\right)^{t}\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{u}}}\left(\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}\right)^{t}\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}\widetilde{\widetilde{\boldsymbol{v}}}\right)^{1/2}}=
\]
\[
=\frac{\widetilde{\widetilde{\boldsymbol{u}}}^{t}\widehat{\boldsymbol{\varSigma}}^{(1,2)}\widetilde{\widetilde{\boldsymbol{v}}}}{\left(\widetilde{\widetilde{\boldsymbol{u}}}^{t}\widehat{\boldsymbol{\varSigma}}^{(1)}\widetilde{\widetilde{\boldsymbol{u}}}\widetilde{\widetilde{\boldsymbol{v}}}^{t}\widehat{\boldsymbol{\varSigma}}^{(2)}\widetilde{\widetilde{\boldsymbol{v}}}\right)^{1/2}}.
\]
Assuming both sample covariance matrices $\widehat{\boldsymbol{\varSigma}}^{(1)}$
and $\widehat{\boldsymbol{\varSigma}}^{(2)}$ are full rank and hence
invertible, this sample correlation can be rewritten in a form allowing
calculation by a direct SVD using the change of variables 
\begin{equation}
\check{\boldsymbol{u}}=\left(\widehat{\boldsymbol{\varSigma}}^{(1)}\right)^{1/2}\widetilde{\widetilde{\boldsymbol{u}}},\;\check{\boldsymbol{v}}=\left(\widehat{\boldsymbol{\varSigma}}^{(2)}\right)^{1/2}\widetilde{\widetilde{\boldsymbol{v}}}.\label{eqP.X:CCAvarChange}
\end{equation}
Note that to arrive at directions that maximize correlation it is
not enough to just make each variable commensurate (i.e. to standardize
each variable by dividing by its standard deviation as done to compute
the correlation matrix discussed in Section \ref{secD4:Standardization}).
Instead standardization\index{standardization} by full root inverse
covariance matrices (sometimes called ``sphering'', as done in Figure
\ref{figJ.B:FLDIntro}) is needed, which essentially results in projected
correlations over all directions in $\mathbb{R}^{d_{1}}$ and $\mathbb{R}^{d_{2}}$. 

The change of variables (\ref{eqP.X:CCAvarChange}) results in the
(reparameterized) sample correlation becoming 
\[
\widehat{\rho}\left(\check{\boldsymbol{u}},\check{\boldsymbol{v}}\right)=\frac{\check{\boldsymbol{u}}^{t}\left(\widehat{\boldsymbol{\varSigma}}^{(1)}\right)^{-1/2}\widehat{\boldsymbol{\varSigma}}^{(1,2)}\left(\widehat{\boldsymbol{\varSigma}}^{(2)}\right)^{-1/2}\check{\boldsymbol{v}}}{\left(\check{\boldsymbol{u}}^{t}\check{\boldsymbol{u}}\check{\boldsymbol{v}}^{t}\check{\boldsymbol{v}}\right)^{1/2}}.
\]
But the SVD optimization seeks unit (i.e. direction) vectors which
have norm 1. Hence to maximize $\widehat{\rho}\left(\check{\boldsymbol{u}},\check{\boldsymbol{v}}\right)$,
it is enough to consider $\check{\boldsymbol{u}}$ and $\check{\boldsymbol{v}}$
with $\left\Vert \check{\boldsymbol{u}}\right\Vert =\left\Vert \check{\boldsymbol{v}}\right\Vert =1$,
and thus to compute the SVD 
\begin{equation}
\left(\widehat{\boldsymbol{\varSigma}}^{(1)}\right)^{-1/2}\widehat{\boldsymbol{\varSigma}}^{(1,2)}\left(\widehat{\boldsymbol{\varSigma}}^{(2)}\right)^{-1/2}=\check{\boldsymbol{U}}\check{\boldsymbol{D}}\check{\boldsymbol{V}}^{t},\label{eqP.AC:SVD4CCA}
\end{equation}
where as above $\check{\boldsymbol{U}}$ and $\check{\boldsymbol{V}}$
are orthonormal matrices and $\check{\boldsymbol{D}}$ is diagonal.
Furthermore, their columns are sets of orthogonal maximizers of $\widehat{\rho}\left(\check{\boldsymbol{u}},\check{\boldsymbol{v}}\right)$.
Inverting the variable change (\ref{eqP.X:CCAvarChange}) returns
those sets of coefficients to the original data scale. In particular,
the matrices of \emph{CCA loadings} are 
\begin{equation}
\widetilde{\widetilde{\boldsymbol{U}}}=\left(\widehat{\boldsymbol{\varSigma}}^{(1)}\right)^{-1/2}\check{\boldsymbol{U}},\;\widetilde{\widetilde{V}}=\left(\widehat{\boldsymbol{\varSigma}}^{(2)}\right)^{-1/2}\check{\boldsymbol{V}}.\label{eqP.AD:CCAloadings}
\end{equation}
Note that while the columns of $\check{\boldsymbol{U}}$ and $\check{\boldsymbol{V}}$
are orthonormal, this generally does not hold for the CCA loadings
$\widetilde{\widetilde{\boldsymbol{U}}}$ and $\widetilde{\widetilde{\boldsymbol{V}}}$.
However, the latter do provide very useful \index{modes of variation!CCA}modes
of variation. In particular, treating the columns of $\widetilde{\widetilde{\boldsymbol{U}}}$
and $\widetilde{\widetilde{V}}$ as direction vectors and finding
the corresponding projection coefficient matrices (again in the spirit
of multiplying by the data blocks as in (\ref{eqP.AI:ScoresAsProd}))
give the \emph{CCA scores} 
\[
\widetilde{\widetilde{\boldsymbol{S}}}^{(1)}=\widetilde{\widetilde{\boldsymbol{U}}}^{t}\breve{\boldsymbol{X}}^{(1)},\;\boldsymbol{S}^{(2)}=\widetilde{\widetilde{\boldsymbol{V}}}^{t}\breve{\boldsymbol{X}}^{(2)}.
\]
An interesting property of these scores is that 
\begin{equation}
\widetilde{\widetilde{\boldsymbol{S}}}^{(1)}\left(\widetilde{\widetilde{\boldsymbol{S}}}^{(1)}\right)^{t}=\widetilde{\widetilde{\boldsymbol{U}}}^{t}\breve{\boldsymbol{X}}^{(1)}\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}\widetilde{\widetilde{\boldsymbol{U}}}=\check{\boldsymbol{U}}^{t}\left(\widehat{\boldsymbol{\varSigma}}^{(1)}\right)^{-1/2}\widehat{\boldsymbol{\varSigma}}^{(1)}\left(\widehat{\boldsymbol{\varSigma}}^{(1)}\right)^{-1/2}\check{\boldsymbol{U}}=I_{n},\label{eqP.Y:CCAsocoresOrtho}
\end{equation}
i.e. the rows of the CCA scores matrix turn out to be orthonormal
vectors in $\mathbb{R}^{n}$. A very parallel calculation shows that
the rows of $\widetilde{\widetilde{\boldsymbol{S}}}^{(2)}$ are also
orthonormal. Thus the transformation\index{transformation|)} (\ref{eqP.X:CCAvarChange})
results in an interesting contrast between CCA and the PLS-SVD introduced
at the beginning of Section \ref{subsecP-B-A:PLS}. In particular,
PLS-SVD has orthonormal loadings (not scores), while CCA has orthonormal
scores and not loadings. This again suggests that the bi-orthogonality
property of PCA modes of variation is indeed a special property. An
interesting variation of CCA that deliberately produces bi-orthogonal\index{bi-orthogonal|)}
(in a particular sense) \index{modes of variation!biorthogonal}modes
of variation has been proposed by \citet{Shu2020D-CCA}.

Various rescalings of the CCA scores appear in different software
packages. A common variation multiplies both sides (\ref{eqP.Y:CCAsocoresOrtho})
by $n$, which gives a scores\index{scores|)} matrix that is interpretable
as having an identity covariance matrix. When variables are not commensurate\index{commensurate|)},
as discussed in Section \ref{secD4:Standardization}, it can be useful
to replace the covariance matrices that are input to CCA with corresponding
correlation matrices. 

The case of covariance matrices not of full rank is sensibly handled
by an application of the methodology called \emph{generalized SVD}.
\index{Singular Value Decomposition!generalized} That also conveniently
allows for lower rank versions of CCA. However, it should be kept
in mind that in high dimensional cases, such methods generally lead
to drastic overfitting\index{overfitting}. This is caused by the
existence of pairs of directions that spuriously give perfect correlation
of the projections, which happens with probability one in high dimensional
cases. Note that such pairs will appear in the analysis ahead of pairs
of directions representing important underlying joint variation. In
such cases PLS is a simple approach that can give results that are
much less sensitive to spurious noise artifacts than CCA. There are
also many other approaches to this problem that use sparsity\index{sparsity}
ideas and other types of regularization. In some elegant unpublished
work, Iain Carmichael uses the generalized SVD approach to formulate
a general framework that encompasses, PLS, CCA and many more. That
paper also includes access to a number of extensions of PLS to the
multi-block (multi-view) case of $B>2$. An earlier survey of multi-block
methods can be found in \citet{Kettenring1971CCA}. 

A special case of CCA, useful in classification\index{classification!relation to CCA}
as discussed in Chapter \ref{chapJ:Classification} is Canonical Variate
Analysis\index{Canonical Variate Analysis}. This generalizes LDA
(see Section \ref{secJ-A:ClassicalMethods}) to the multi-class case,
by computing CCA with one data block consisting of row indicators
(zeros and ones) for each class.

A toy example highlighting the difference between CCA and PLS-SVD
is given in Figure \ref{figP.D:PCA-PLS-CCA}. This is a two block
example, with the first data block containing data object vectors
of the form $\left(\begin{array}{cc}
x_{1} & x_{2}\end{array}\right)^{t}\in\mathbb{R}^{2}$. The second block consists of scalars $y\in\mathbb{R}$. The top
panel displays the relationship between the data blocks using two
rotated views, where each of the $n=100$ data objects appears as
a point of the form $\left(\begin{array}{ccc}
x_{1} & x_{2} & y\end{array}\right)^{t}\in\mathbb{R}^{3}$ indicated with a plus sign. There is a very strong relationship between
data blocks in the sense that the points lie entirely in the cyan
plane. As shown in the lower right panel, the data are distributed
within the cyan plane to have a strong first PC direction shown in
red in all panels. That red PC1 direction has been chosen to also
lie within the first data block $\mathbb{R}^{2}$ space, shown as
the yellow plane. Both PLS and CCA seek optimal pairs of directions,
one in each space. Because the second data block lies in $\mathbb{R}$
there is only one possible direction, so the optimizing direction
is just the vertical axis. More informative are the PLS and CCA directions
in the second block (which thus lie in the yellow plane), shown in
red and green respectively. The CCA direction points in the direction
of steepest increase of the cyan plane, which is natural because that
will maximize the projected correlation (the target of CCA). On the
other hand, the PLS direction is a compromise between the CCA and
PC 1 directions, because it optimizes covariance which responds also
to the very strong variance in the PC 1 direction.

\begin{figure}[H]
\includegraphics[scale=0.75]{ChapterP/OODAbookChpPFigDPLSvCCA1A}\includegraphics[scale=0.75]{ChapterP/OODAbookChpPFigDPLSvCCA1B}

\includegraphics[scale=0.85]{ChapterP/OODAbookChpPFigDPLSvCCA2}

\caption{Toy example highlighting the difference between PLS (green) and CCA
(blue) directions. The concatenated ($\mathbb{R}^{3}$) data points
all lie in the cyan plane. The $\breve{\boldsymbol{X}}^{(1)}$ data
block lies in the yellow plane, as do the PLS-SVD and CCA directions.
This shows how CCA maximizes correlation with the $\breve{\boldsymbol{X}}^{(2)}$
data that lie in the vertical axis, while the PLS direction is strongly
influence\index{influence|)}d by the red PC1 direction. \label{figP.D:PCA-PLS-CCA} }
\end{figure}
A range of compromise directions between PCA\index{Principal Component Analysis|)},
PLS\index{Partial Least Squares|)} and CCA have been proposed in
the unpublished PhD dissertation of \citet{Lee2007ContinuumCCA},
based on the \emph{continuum regression\index{continuum regression}}
ideas of \citet{StoneBrooks1990ContReg}.

CCA also has a strong relationship to classical linear regression.
Using a small circle above symbols to distinguish that notation, the
linear regression model is usually written as 
\[
\mathring{\boldsymbol{y}}=\mathring{\boldsymbol{X}}\mathring{\boldsymbol{\beta}}+\mathring{\boldsymbol{\epsilon}},
\]
where $\mathring{\boldsymbol{y}}$ is an $n\times1$ vector of responses,
$\mathring{\boldsymbol{X}}$ is an $n\times d$ design matrix whose
columns are thought of as multivariate predictors, $\mathring{\boldsymbol{\beta}}$
is a $d\times1$ vector of regression coefficients, and $\mathring{\boldsymbol{\epsilon}}$
is an $n\times1$ vector of residuals. Direct connection with the
above CCA context comes from assuming the data have been centered
in the sense that the entries of $\mathring{\boldsymbol{y}}$ have
sample mean 0, as do the columns of $\text{\ensuremath{\mathring{\boldsymbol{X}}}}$.
That allows using the column object centered notation (\ref{eqP-G:ColObjCentData})
to write $\breve{\boldsymbol{X}}^{(1)}=\frac{1}{\sqrt{n}}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}$
and $\breve{\boldsymbol{X}}^{(2)}=\frac{1}{\sqrt{n}}\mathring{\boldsymbol{y}}^{t}$.
This gives inputs to CCA as 
\[
\widehat{\boldsymbol{\varSigma}}^{(1)}=\breve{\boldsymbol{X}}^{(1)}\left(\breve{\boldsymbol{X}}^{(1)}\right)^{t}=\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{X},}
\]
\[
\widehat{\boldsymbol{\varSigma}}^{(1,2)}=\breve{\boldsymbol{X}}^{(1)}\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}=\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{y}},
\]
\[
\widehat{\boldsymbol{\varSigma}}^{(2)}=\breve{\boldsymbol{X}}^{(2)}\left(\breve{\boldsymbol{X}}^{(2)}\right)^{t}=\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{y}}}}^{t}\mathring{\boldsymbol{y}}=\hat{\sigma}_{Y}^{2},
\]
where the scalar $\hat{\sigma}_{Y}^{2}$ is an estimate of the standard
deviation of the entries of $\mathring{\boldsymbol{y}}$. The corresponding
SVD to compute CCA as in (\ref{eqP.AC:SVD4CCA}) is thus applied to
the $d\times1$ matrix 
\[
\left(\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{X}}\right)^{-1/2}\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{y}}\hat{\sigma}_{Y}^{-1}.
\]
That SVD\index{Singular Value Decomposition} results in a $d\times1$
vector of norm 1 
\[
\check{\boldsymbol{U}}=\frac{\left(\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{X}}\right)^{-1/2}\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{y}}\hat{\sigma}_{Y}^{-1}}{\left\Vert \left(\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{X}}\right)^{-1/2}\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{y}}\hat{\sigma}_{Y}^{-1}\right\Vert },
\]
the single singular value 
\[
\check{\boldsymbol{D}}=\left\Vert \left(\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{X}}\right)^{-1/2}\frac{1}{n}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{y}}\hat{\sigma}_{Y}^{-1}\right\Vert 
\]
 and just the scalar $\check{\boldsymbol{V}}=1$. Inverting the change
of variables as in (\ref{eqP.AD:CCAloadings}) and dividing by the
length of the resulting vector gives a CCA loadings\index{loadings|)}
vector which points in the same direction in $\mathbb{R}^{d}$ as
the familiar least squares regression estimate $\hat{\boldsymbol{\beta}}=\left(\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{X}}\right)^{-1}\text{\ensuremath{\mathring{\boldsymbol{X}}}}^{t}\mathring{\boldsymbol{y}}$
of the coefficient vector $\mathring{\boldsymbol{\beta}}$. An important
consequence is that the direction in the object space\index{object space}
$\mathbb{R}^{d}$, where a set of data objects $\widetilde{\boldsymbol{X}}$
best correspond to a set of $n$ values in the vector $\widetilde{\boldsymbol{y}}$,
is proportional to the least squares estimate $\hat{\boldsymbol{\beta}}$
(assuming $\widetilde{\boldsymbol{X}}$ has been column object mean\index{column object mean|)}
centered\index{mean!center|)} and that the mean of the $\widetilde{\boldsymbol{y}}$
values is 0). \index{Canonical Correlation Analysis|)}

\subsection{Joint and Individual Variation Explained\label{subsecP-B-C:JIVE}}

The methods described in Sections \ref{subsecP-B-A:PLS} and \ref{subsecP-B-B:CCA}
all focus on the investigation of \emph{joint variation}\index{joint variation},
i.e. finding \index{modes of variation!joint|)}modes of variation\index{modes of variation!individual}
that reflect how data blocks vary together with each other. But in
many applications there is also keen interest in \emph{individual
variation}\index{individual variation}, i.e. variation that is block
specific and in some sense independent (either stochastically, or
else in terms of appropriate orthogonality of modes of variation)
of the other blocks. That idea was formalized as a method called Joint
and Individual Variation Explained (JIVE) by \citet{Lock2013JIVE}.
An approach based on the simpler and more direct linear algebra method
of \emph{\index{Principal Angle Analysis}Principal Angle Analysis}
is the Angle-based JIVE of \citet{Feng2018JIVE}. \citet{Yu2017JIVE}
and \citet{Carmichael2019JIVEimage} used that method to demonstrate
the value of careful investigation of both joint and individual variation
in neuroscience and cancer genomics/imaging applications. An approach
to the case of \emph{\index{partially shared blocks}partially shared
blocks} (allowing joint variation between various subsets of data
blocks) was proposed by \citet{GaynanovaLi2019MultiView}.\index{modes of variation|)\index{}}
