# Regularised Canonical Correlation Analysis: graphical lasso, biplots and
beyond

Lennie Wells$\hskip 0.85358pt{}^{1}$, Kumar Thurimella$\hskip 0.85358pt{}^{2}$
& Sergio Bacallado$\hskip 0.85358pt{}^{1}$

###### Abstract

Recent developments in regularized Canonical Correlation Analysis (CCA)
promise powerful methods for high-dimensional, multiview data analysis.
However, justifying the structural assumptions behind many popular approaches
remains a challenge, and features of realistic biological datasets pose
practical difficulties that are seldom discussed. We propose a novel CCA
estimator rooted in an assumption of conditional independencies and based on
the Graphical Lasso. Our method has desirable theoretical guarantees and good
empirical performance, demonstrated through extensive simulations and real-
world biological datasets. Recognizing the difficulties of model selection in
high dimensions and other practical challenges of applying CCA in real-world
settings, we introduce a novel framework for evaluating and interpreting
regularized CCA models in the context of Exploratory Data Analysis (EDA),
which we hope will empower researchers and pave the way for wider adoption.

Keywords: Canonical Correlation Analysis; Graphical Lasso; Multi-Omics; Multi-
view data.

††1Department of Pure Mathematics and Mathematical Statistics, University of
Cambridge††2Department of Chemical Engineering and Biotechnology, University
of Cambridge

##  1 Introduction

Canonical correlation analysis (CCA) is a powerful tool in multivariate
statistics and multiview learning. It provides a principled approach to
understand the correlation structure between two datasets or views of data.
CCA is closely related to the better known technique of PCA (Jolliffe, 1986;
Jolliffe and Cadima, 2016), and has similar linguistic ambiguities: the term
CCA can refer to a population canonical decomposition (a mathematical object
determined by the law of a pair of random vectors), to a related sample CCA
object (determined by a pair of data matrices), or to associated frameworks
for data analysis. Moreover, there are many different interpretations of CCA;
leading to a great number of extensions and alternative methods (Roweis and
Ghahramani, 1999; Bach and Jordan, 2005; Klami et al., 2013; Tenenhaus et al.,
2014; Chapman and Wang, 2021; Chapman et al., 2023). A more complete review of
this wider context is given in Appendix F.

This article takes a frequentist perspective of high-dimensional CCA, where we
assume independent and identically distributed observations from some unknown
distribution and wish to estimate certain aspects of the population canonical
decomposition. Unfortunately, the classical notion of sample CCA breaks down
when there are more dimensions than samples, and in response, many regularised
extensions of CCA have been proposed. Some of these extensions are now widely
used, particularly in the fields of genetics (Witten et al., 2009) and
neuroscience (Mihalik et al., 2022).

Even in this narrower context, CCA can be used to accomplish many different
goals. One very popular goal is dimension reduction of high-dimensional data:
to obtain low dimensional representations of the samples that preserve salient
information and can be used for downstream tasks, such as classification,
regression, and visualisation. Though our results can be applied to this use-
case of dimension reduction, our statistical approach to CCA places more
emphasis on the variables than the samples. Indeed, our main goal is to
understand the relationship between variables in the population distribution.
In this context, CCA can be a powerful tool for Exploratory Data Analysis
(EDA) and provide interpretable conclusions to guide scientific enquiry.

Before we can outline our contributions in more detail we need to define CCA
and related terminology. Please note that our choice of language and notation
may differ from any given source, due to a lack of consistency in the
literature where terms such as _loading_ can refer to various objects.

###  1.1 Terminology for CCA

To formalise the population CCA problem, suppose we have two random variables,
$X,Y$ taking values in $\mathbb{R}^{p},\mathbb{R}^{q}$ respectively with joint
covariance partitioned as

$\Sigma=\left(\begin{array}[]{cc}\Sigma_{xx}&\Sigma_{xy}\\\
\Sigma_{yx}&\Sigma_{yy}\end{array}\right)\;.$ (1)

CCA defines successive vectors $u_{k}\in\mathbb{R}^{p},v_{k}\in\mathbb{R}^{q}$
by the programs

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\,&u^{\top}\Sigma_{xy}v\\\
\text{subject to }&u^{\top}\Sigma_{xx}u\leq 1,\>v^{\top}\Sigma_{yy}v\leq 1,\\\
&u^{\top}\Sigma_{xx}u_{j}=v^{\top}\Sigma_{yy}v_{j}=0\text{ for }1\leq j\leq
k-1\;.\end{split}$ (2)

We call the optimal value $\rho_{k}$ the $k^{\text{th}}$ canonical
correlation, call $u_{k},v_{k}$ the $k^{\text{th}}$ pair of canonical
directions, or simply weights and call the projections
$u_{k}^{\top}X,v_{k}^{\top}Y$ the first pair of canonical variates. We shall
always refer to the components of the original $X,Y$ as variables, and to the
transformed variables as (canonical) variates. The quantities
$\Sigma_{xx}u_{k},\Sigma_{yy}v_{k}$ will also be important; we will refer to
these as canonical loading vectors, or simply loadings. In particular the
$k^{\text{th}}$ pair maximises
$\operatorname{Corr}(u_{k}^{\top}X,v_{k}^{\top}Y)$ subject to the
orthogonality constraints. When $\Sigma_{xx},\Sigma_{yy}$ are both full rank,
we can define $\min\\{p,q\\}$ canonical direction pairs by Equation 2 and then
extend this to give a pair of bases $(u_{k})_{k},(v_{k})_{k}$, for
$\mathbb{R}^{p},\mathbb{R}^{q}$, respectively. The situation is more subtle
when $\Sigma_{xx},\Sigma_{yy}$ are not of full rank; we discuss this case in
detail in Appendix B. We shall loosely use canonical decomposition to refer to
this full collection of vectors, or to the associate collection of canonical
variates.

In practice, a statistician will not have access to the population covariance
matrices, but only samples. We assume these data are realisations of
independent and identically distributed copies of $(X,Y)$. The sample CCA
problem is to estimate certain aspects of the canonical decomposition, for
example the first few canonical direction pairs, from these samples. The
classical method for the sample problem replaces the true covariance matrices
in (2) with sample covariance matrices; see Section 2.2.

###  1.2 Initial contribution: graphical CCA (gCCA)

Existing regularised CCA methods assume regularity of the canonical directions
and apply penalties to enforce this. The main alternatives are ridge CCA
(rCCA) (Vinod, 1976), corresponding to an $\ell_{{2}}$ penalty on the
directions, and a family of sparse CCA (sCCA) methods (Gao et al., 2016; Mai
and Zhang, 2019), which penalise the $\ell_{{1}}$-norm of the canonical
directions; $\ell_{{1}}$ penalties are widely used in high-dimensional
statistics to promote sparsity (Hastie et al., 2015; Wainwright, 2019). A
family of methods based on the earlier work of Witten et al. (2009) is also
often branded as sparse CCA, even though the orthogonality constraints are
defined rather differently; following Mihalik et al. (2022) we call these
sparse Partial Least Squares (sPLS) methods. It is important to consider sPLS
methods because they are widely applied, with appealing open-source
implementations. We define versions of rCCA, sCCA and sPLS in Section 3, with
further discussion and a wider literature review in Section F.1.

Theoretical work has shown that certain methods can effectively estimate the
canonical directions when the true underlying directions are sparse (Gao et
al., 2016; Mai and Zhang, 2019). However, it is hard to justify why sparse
directions should be assumed in real world datasets. We posit that sparsity of
the population precision matrix can be a more natural structural assumption,
derived from a high degree of conditional independence in a Gaussian model;
see Section 4.1. Thus, we propose graphical CCA (gCCA), a plug-in method
exploiting the powerful Graphical Lasso estimator (Yuan and Lin, 2007;
Friedman et al., 2008). Full details are given in Section 4, which includes
non-asymptotic upper bounds on the estimation error. However, these
theoretical results are only as informative as their assumptions permit; to
understand whether gCCA is better in practice, we need to compare it to
existing regularisation CCA methods on real data.

###  1.3 Main contribution: frameworks for comparison and interpretation

In real-world applications, regularised CCA presents several hurdles to the
practitioner. First, an ever-growing toolbox of estimators demands model
comparison. Tuning penalty parameters and choosing the number of canonical
directions (essentially, the variate subspace dimension) add further layers of
complexity. There is little guidance in the literature on how to choose
appropriate criteria for evaluation, which will depend on the number of
canonical correlations under consideration. The real challenge lies in
realistic datasets, which deviate from idealized models used for theoretical
analysis and rarely exhibit clear gaps in the canonical correlations. Most
commonly, practitioners choose the number of directions arbitrarily or only
consider the top pair, leading to potential inconsistencies between methods
and unreliable conclusions, especially in data-scarce settings. This
underscores the urgent need for systematic analysis frameworks which allow the
practitioner to evaluate and contrast a range of regularised CCA estimators.

The relevance of different evaluation criteria depends on the downstream task
of interest. To explore this issue, we consider a broad family of evaluation
criteria, which we introduce in Section 5. This family includes two main
classes: criteria for correlation signal captured, and criteria for estimation
accuracy (for directions, variates and loadings). For each such class, one can
consider multiple successive directions individually or as successive
subspaces. Moreover, each class of criteria can be applied in the oracle case
where the true distribution is known (relevant for synthetic data) or in the
empirical case where the true distribution is unknown and sample splitting is
required (relevant for real data).

Our motivation for graphical CCA originated from a project studying the human
gut microbiome with the hope of generating insights to help understand and
treat Irritable Bowel Disease (IBD). We wished to apply CCA to explore the
relationship between metabolites and enzymes in the human gut. The structure
of biochemical pathways provides strong motivation for expecting conditional
independencies in the data. We use a dataset related to this problem for
illustrative purposes, both in its original form, and to provide meaningful
constructions for synthetic data using a parametric bootstrap approach. We
shall refer to this as the ‘Microbiome’ dataset in the sequel, and give
further details in section 8.1.

In Section 6 we systematically compare different CCA methods on such synthetic
data using this broad family of criteria. The parametric bootstrap model
derived from our Microbiome dataset exhibits markedly different
characteristics to illustrative models typically studied in the regularised
CCA literature. We make the following main observations:

  * •

gCCA works well. It outperforms rCCA and sCCA across a wide range criteria.
Additionally, as might be expected, sPLS performs significantly worse than the
three CCA methods across both correlation and estimation criteria.

  * •

Cross Validation (CV) works well. The CV-selected tuning parameter is usually
near the oracle-optimal one for any given correlation objective; moreover such
CV parameters often give near-oracle-optimal performance across a range of
correlation and estimation criteria.

  * •

Correlation is easier than estimation. Many successive directions with large
correlation signal are recovered, even in high dimensions, but these are not
always near the true canonical directions.

  * •

Variates and loadings are easier to estimate than weights and successive
subspaces are easier to estimate than successive pairs. In fact, very
different sequences of weights can lead to very similar variate subspaces.
However, it is important to use comparable accuracy metrics.

These observations have significant practical implications. The former two
observations reassure us that CCA is feasible in certain high-dimensional
practical settings. When using CCA for EDA, we would ideally only wish to
spend time closely inspecting aspects of the estimated canonical decomposition
which have good statistical properties. The latter two observations suggest
that there may be situations where it is sensible to closely inspect the
loading vectors but not the weight vectors. These observations are consistent
with the geometry of CCA, previous theoretical results (Gao et al., 2016; Ma
and Li, 2020) and previous empirical studies using CCA (Gu and Wu, 2018).
Evaluating estimation accuracy on real data is challenging, but accurate
estimators should certainly be stable to sample splitting and small changes in
hyperparameter choice. Stability also to the choice of regularisation gives
confidence that the observations are due to the signal in the data rather than
the chosen structural assumptions.

In Section 8 we propose a novel framework for applying CCA in practice,
leveraging these observations. This framework allows a practitioner to fit
many different CCA methods with a range of penalty parameters and then compare
the different estimates. It includes visualisations to assess which estimates
capture significant correlation and determine what aspects of those estimates
are stable to sample splitting, hyperparameter, or choice of regularisation.
We suggest that a practitioner only focus on stable aspects of the estimates
for closer inspection.

To aid such closer inspection we propose a powerful ‘biplot’ visualisation
method extending the proposal of ter Braak (1990). These biplots show
correlations between the original variables and the canonical variates, but
can also be interpreted through a latent variable formulation of CCA. They
geometry of the biplots reflects the subspace nature of the CCA problem with a
natural correspondence between registration of the estimates and
transformations in the biplot. The biplots are particularly insightful in
three dimensions, and when analysed interactively with modern plotting
libraries. All our code is available on GitHub at
https://github.com/W-L-W/RegularisedCCA, and we recommend the reader
experiment with the Jupyter notebooks provided. We hope that the framework we
introduce may lead to high-dimensional CCA becoming more widely used for EDA.

###  1.4 How to read this paper

Our initial contribution is in Section 4; to put this into context we present
existing methods in Section 3. The main numerical experiments on synthetic
data are in Section 6 with supporting definitions in Section 5; experiments on
real data are in Section 8 with supporting definitions in Section 7. However,
we first present essential background and notation for the CCA problem in
Section 2.

We include only minimal mathematical statements in the main text; full details
are given in the supplement. In particular, Appendices A and B provides a
self-contained treatment of the mathematical foundations of CCA, including
technical results which are often stated without proof or citation. This
includes versions of all the results we present in the following Section 2. We
hope this may provide a useful reference for future work.

##  2 CCA Background and Notation

###  2.1 Formulation via Singular Value Decomposition

CCA is well known to be equivalent to a certain Singular Value Decomposition
(SVD). The SVD is a very well studied object in matrix analysis (Stewart and
Sun, 1990), with many attractive geometrical properties, which can be
translated to CCA. In particular, suppose the within-view covariance matrices
are invertible, and define

$T\vcentcolon=\Sigma_{xx}^{-1/2}\Sigma_{xy}\Sigma_{yy}^{-1/2}\;.$ (3)

Then the singular vectors of $T$ are $\Sigma_{xx}^{1/2}u_{k}$ and
$\Sigma_{yy}^{1/2}v_{k}$ and the singular values are precisely $\rho_{k}$
where $u_{k},v_{k},\rho_{k}$ are the $k^{\text{th}}$ canonical pair and
canonical correlation respectively. This follows easily from a standard
formulation of SVD, as outlined in Appendix B. An immediate consequence of
this equivalence is the identity

$\Sigma_{xy}=\Sigma_{xx}^{1/2}\,T\,\Sigma_{yy}^{1/2}=\Sigma_{xx}\left(\sum_{k=1}^{K}\rho_{k}u_{k}v_{k}^{\top}\right)\Sigma_{yy}$
(4)

for any $K\leq\min(p,q)$ such that $\rho_{k}=0$ for each $k>K$.

###  2.2 Sample CCA

Write $\mathbf{X}=(x_{1},\dots,x_{N})^{\top}\in\mathbb{R}^{N\times
p},\mathbf{Y}=(y_{1},\dots,y_{N})^{\top}\in\mathbb{R}^{N\times q}$ for our
data-matrices. We follow the convention, inspired by Hastie et al. (2015),
that vectors and matrices of samples (with leading dimension $N$) all appear
in bold. This includes the columns of the data matrices, which we write as
$\mathbf{X}_{i}=((x_{n})_{i})_{n=1}^{N}\in\mathbb{R}^{N}$ for $i\in 1,\dots,p$
and $\mathbf{Y}_{j}$ for $j\in 1,\dots,q$. We assume the data is mean-centred.
When implementing CCA algorithms in practice, there are many different ways
one might apply this mean-centering, but this usually has a only a small
impact on the results, so we will not discuss the issue further. With mean-
centred variables, one may use the shorthand for sample covariance matrices:

$\mathbf{C}=\left(\begin{array}[]{ll}\mathbf{C}_{xx}&\mathbf{C}_{xy}\\\
\mathbf{C}_{yx}&\mathbf{C}_{yy}\end{array}\right)=\left(\begin{array}[]{ll}\tfrac{1}{N}\mathbf{X}^{\top}\mathbf{X}&\tfrac{1}{N}\mathbf{X}^{\top}\mathbf{Y}\\\
\tfrac{1}{N}\mathbf{Y}^{\top}\mathbf{X}&\tfrac{1}{N}\mathbf{Y}^{\top}\mathbf{Y}\end{array}\right)\;.$
(5)

Hence the classical CCA estimator solves

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\,&u^{\top}\mathbf{C}_{xy}v\\\
\text{subject to }&u^{\top}\mathbf{C}_{xx}u\leq
1,\>v^{\top}\mathbf{C}_{yy}v\leq 1,\\\
&u^{\top}\mathbf{C}_{xx}u_{j}=v^{\top}\mathbf{C}_{yy}v_{j}=0\text{ for }1\leq
j\leq k-1\;.\end{split}$ (6)

###  2.3 Variate-based formulation of CCA

Rewriting the program from Equation 2 in terms of covariances between (the
canonical) variates rather than covariance matrices yields

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\,&\operatorname{Cov}(u^{\top}X,v^{\top}Y)\\\
\text{subject to }&\operatorname{Var}(u^{\top}X)\leq
1,\>\operatorname{Var}(v^{\top}Y)\leq 1,\\\
&\operatorname{Cov}(u^{\top}X,u_{j}^{\top}X)=\operatorname{Cov}(v^{\top}Y,v_{j}^{\top}Y)=0\text{
for }1\leq j\leq k-1\;.\end{split}$ (7)

Before giving a natural sample analogue, we need to introduce notation for
sample covariances and correlations. Let
$\mathbf{e},\mathbf{f}\in\mathbb{R}^{N}$ be vectors with mean zero. Then we
define

$\displaystyle\widehat{\operatorname{Cov}}(\mathbf{e},\mathbf{f})\vcentcolon=\tfrac{1}{N}\mathbf{e}^{\top}\mathbf{f},\qquad\widehat{\operatorname{Var}}(\mathbf{e})\vcentcolon=\tfrac{1}{N}\mathbf{e}^{\top}\mathbf{e}=\widehat{\operatorname{Cov}}(\mathbf{e},\mathbf{e}),\qquad\widehat{\operatorname{Corr}}(\mathbf{e},\mathbf{f})\vcentcolon=\frac{\mathbf{e}^{\top}\mathbf{f}}{\lVert\mathbf{e}\rVert\lVert\mathbf{f}\rVert}=\frac{\widehat{\operatorname{Cov}}(\mathbf{e},\mathbf{f})}{\sqrt{\widehat{\operatorname{Var}}(\mathbf{e})\widehat{\operatorname{Var}}(\mathbf{f})}}\;.$

Then the sample analogue of Equation 7 is

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\,&\widehat{\operatorname{Cov}}(u^{\top}\mathbf{X},v^{\top}\mathbf{Y})\\\
\text{subject to }&\widehat{\operatorname{Var}}(u^{\top}\mathbf{X})\leq
1,\>\widehat{\operatorname{Var}}(v^{\top}\mathbf{Y})\leq 1,\\\
&\widehat{\operatorname{Cov}}(u^{\top}\mathbf{X},u_{j}^{\top}\mathbf{X})=\widehat{\operatorname{Cov}}(v^{\top}\mathbf{Y},v_{j}^{\top}\mathbf{Y})=0\text{
for }1\leq j\leq k-1\;.\end{split}$ (8)

######  Remark 1.

It is now easy to see how CCA breaks down when $q\geq n$. Indeed, suppose
$\mathbf{Y}$ is of full row rank; then for any $u\in\mathbb{R}^{p}$ there is
some $v\in\mathbb{R}^{q}$ such that $\mathbf{Y}v=\mathbf{X}u$, and therefore
$\widehat{\operatorname{Corr}}(\mathbf{X}u,\mathbf{Y}v)=1$. Therefore all the
canonical correlations are in fact 1.

######  Remark 2 (Random variables to vectors of samples).

Here we saw that the sample CCA version could be formulated by replacing
random variables with sample vectors and corresponding population covariance
operators with sample covariance operators. We will see this repeatedly in the
rest of the article.

###  2.4 Further notation and background

It will often be helpful to use the shorthand

$U_{K}=\left(u_{1},\dots,u_{K}\right)\in\mathbb{R}^{p\times K},\quad
V_{K}=\left(v_{1},\dots,v_{K}\right)\in\mathbb{R}^{q\times K}$ (9)

for matrices formed by concatenation of successive canonical directions.

We shall also use the notion of canonical angles. These are a sequence of
angles that reflect how well aligned two subspaces of a vector space are, and
are well-studied in matrix analysis (Stewart and Sun, 1990). They are used to
define the $\sin\Theta$ distance and $\cos\Theta$ similarity between
subspaces, which we will use to evaluate how close estimated subspaces are
from a ground truth or between each other. Definitions and an extended
background may be found in Section A.3.

##  3 Existing Regularised Sample CCA Methods

A large number of methods for regularised CCA have been proposed. These
methods generally add some regularisation to (6) to obtain an objective
function whose optimum defines an estimator. For brevity we only describe
three methods that are representative of most of these approaches; we will
compare our method with these representative methods in Sections 6 and 8. We
give a more thorough review and discussion of the many variants and extensions
of these representative methods in Appendix F.

###  3.1 Ridge-regularised CCA (rCCA)

Ridge regularisation for CCA (Vinod, 1976) is the oldest and most established
regularised CCA method (Uurtio et al., 2017). For our purposes, it is most
convenient to define estimated canonical directions through the successive
programs

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\,&u^{\top}\mathbf{C}_{xy}v\\\
\text{subject to
}&u^{\top}\left((1-c_{x})\mathbf{C}_{xx}+c_{x}I_{p}\right)u\leq
1\;,\>v^{\top}\left((1-c_{y})\mathbf{C}_{yy}+c_{y}I\right)v\leq 1\;,\\\
&u^{\top}\left((1-c_{x})\mathbf{C}_{xx}+c_{x}I_{p}\right)u_{j}=v^{\top}\left((1-c_{y})\mathbf{C}_{yy}+c_{y}I\right)v_{j}=0\text{
for }1\leq j\leq k-1\;,\end{split}$ (10)

where $c_{x},c_{y}\in[0,1]$ are tuning parameters. The estimator can therefore
be viewed as a plug-in estimator, corresponding to regularised covariance
estimate

$\displaystyle\hat{\Sigma}^{\text{ridge}}=\left(\begin{array}[]{ll}(1-c_{x})\mathbf{C}_{xx}+c_{x}I_{p}&\mathbf{C}_{xy}\\\
\mathbf{C}_{yx}&(1-c_{y})\mathbf{C}_{yy}+c_{y}I\end{array}\right)\;.$

Note that taking $c_{x}=c_{y}=0$ recovers the classical population CCA. By
contrast, taking $c_{x}=c_{y}=1$ replaces the within-view covariance matrices
by identity matrices; analogously to the discussion in Section 2.1 this
precisely to an SVD of the sample cross-covariance matrix $\mathbf{C}_{xy}$.
We follow Mihalik et al. (2022) and refer to this case as Partial Least
Squares (PLS). In fact there is a whole family of methods referred to as PLS
in the applied statistics literature; these are all closely related to the SVD
of $\mathbf{C}_{xy}$ but often defined by iterations, rather than explicit
objectives. We refer the interested reader to Rosipal and Krämer (2006) for a
broader review.

###  3.2 Sparse Partial Least Squares (sPLS)

One of the most common methods for sparse CCA used in practice is the
Penalised Matrix Decomposition (PMD) approach proposed in Witten et al. (2009)
and implemented in the R package PMA. In fact PMD is a more general method for
computing low rank approximations to matrices, which can also be applied other
problems such as sparse Principal Component Analysis (PCA). In the context of
multiview data, their method estimates the $k^{\text{th}}$ pair of canonical
directions by optimisers of

$\begin{split}\operatorname{maximize}_{u,v}\>&u^{\top}\mathbf{C}_{xy}^{(k)}v\\\
\text{ subject to}\>&\|u\|_{2}^{2}\leqslant 1,\enspace\|v\|_{2}^{2}\leqslant
1,\\\ &\|u\|_{1}\leqslant s_{x},\enspace\|v\|_{1}\leqslant
s_{y}\;,\end{split}$ (11)

where

$\mathbf{C}_{xy}^{(k)}=\mathbf{C}_{xy}-\sum_{l=1}^{k-1}d_{l}u_{l}v_{l}^{\top},\quad
d_{l}=u_{l}^{\top}\mathbf{C}_{xy}v_{l}$

are deflated versions of the cross-covariance. This deflation encourages
orthogonality between successive estimates, rather than being forced
explicitly.

An important difference between this program (11) and the sample CCA program
(8) is that the constraints
$u_{k}^{\top}\mathbf{C}_{xx}u_{k}=v_{k}\mathbf{C}_{yy}v_{k}=1$ in have been
replaced by $u_{k}^{\top}u_{k}=v_{k}^{\top}v_{k}=1$. This is mainly for
computational convenience. The justification is that the two quantities are
equivalent when $\mathbf{C}_{xx}=I_{p},\mathbf{C}_{yy}=I_{q}$; the authors
claim that this is a reasonable approximation in certain high-dimensional
genomics applications, where within view population covariances $\Sigma_{xx}$,
$\Sigma_{yy}$ may be assumed to be near identity.

However, in cases where this approximation is poor, the method does not
correspond to CCA. Indeed, on removing the $\ell_{{1}}$ penalty terms the
programs recover the SVD of $\mathbf{C}_{xy}$, and so correspond to PLS.
Following Mihalik et al. (2022), we therefore refer to this method as sparse
PLS (sPLS).

###  3.3 Sparse CCA (sCCA)

Many works have since proposed computationally feasible CCA estimators which
penalise the $\ell_{1}$-norm of canonical directions. For brevity we will only
compare to the method proposed in Suo et al. (2017); this is arguably the most
natural generalisation of sample CCA (2) and is faster or more performant than
existing alternatives. This sCCA method defines estimated canonical directions
through the successive programs

$\begin{split}\operatorname{minimize}_{u,v}\;&-u^{\top}\mathbf{C}_{xy}v+\lambda_{x}\|u\|_{1}+\lambda_{y}\|v\|_{1}\\\
\text{ subject to}\;&u^{\top}\mathbf{C}_{xx}u\leqslant 1\;,\enspace
v^{\top}\mathbf{C}_{yy}v\leqslant 1\;,\\\
&u^{\top}\mathbf{C}_{xx}u_{j}=v^{\top}\mathbf{C}_{yy}v_{j}=0\text{ for }1\leq
j\leq k-1\,.\end{split}$ (12)

###  3.4 Notes on implementation and application

Algorithmic details for solving the programs Equations 10, 11 and 12 are only
of secondary interest in this article. It is important that algorithms are
reliable and fast enough to run on a large dataset for many different tuning
parameter values in a reasonable amount of time, but beyond this point the
limiting factor is human ability to interpret the output.

We now briefly outline the algorithms we used in the forthcoming numerical
experiments. The rCCA program Equation 10 is in fact equivalent to population
CCA on the regularised joint covariance matrix so can be solved with standard
numerical algorithms for CCA; we used the efficient implementation in the
valuable cca-zoo python package (Chapman and Wang, 2021). We implemented the
remaining two methods ourselves. The sPLS program Equation 11 can be solved
efficiently by iterative soft-thresholding and rescaling (Witten et al.,
2009). The sCCA program Equation 12 is harder to solve. The original algorithm
from Suo et al. (2017) used alternating convex search (ACS), i.e. iteratively
optimise over $u$ with $v$ fixed and $v$ with $u$ fixed. We where able to
speed up this algorithm significantly by interlacing partial updates for $u$
and $v$; we give full details and further background in Appendix J.

One potential disadvantage of the sPLS and sCCA programs is that they have
non-convex objectives. Therefore there may be many local minima, and the
solution returned may depend on the initialisation used. In practice, we did
not find this to a problem: reasonable choices of initialisation found local
minima which were ‘good’ even if not globally optimal.

Another potential problem with the sPLS and sCCA programs is that they have a
potentially large number of tuning parameters. Though this provides
flexibility, it is very expensive to perform tuning parameter selection over a
high-dimensional grid of parameters; even a two dimensional parameter grid
will be significantly more expensive than a one parameter family. We therefore
tied all the tuning parameters to be the same in our numerical experiments;
this had the added advantage of easier comparability with rCCA.

We next define our proposed CCA method: graphical CCA (gCCA). As we will see,
it avoids these two potential disadvantages: it has pleasant optimisation
properties from a convex formulation and only a single tuning parameter.

##  4 Our method: Graphical CCA (gCCA)

###  4.1 Background for the Graphical Lasso

Graphical models are a well-studied class of probabilistic models that use a
graph structure to constrain the distribution of a random vector and have
applications across statistics and machine learning (Bishop, 2006). These
models associate a vector-valued random variable $Z=(Z_{1},\dots,Z_{\bar{p}})$
with the vertex set $\mathcal{V}=\\{1,\dots,\bar{p}\\}$ of some graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$, but the probabilistic constraint can
take a variety of different forms, depending on the statistical setting
(Lauritzen, 1996). A central statistical problem is to estimate the edge set
$\mathcal{E}$ from many independent observations of $Z$ and is called
structure estimation.

Gaussian graphical models have particularly pleasant properties (Lauritzen,
1996). These models are associated with undirected graphs and constrain the
probability distribution so that pairs of variables $Z_{i},Z_{j}$ are
conditionally independent given $(Z_{k})_{k\neq i,j}$ whenever
$(i,j)\notin\mathcal{E}$. It is a fundamental result (Lauritzen, 1996,
Proposition 5.2) that such conditional independencies correspond to zeros in
the precision matrix $\Omega=\operatorname{Var}(Z)^{-1}$:

######  Lemma 1.

Let $Z\in\mathbb{R}^{\bar{p}}$ be a Gaussian random vector with (invertible)
precision matrix $\Omega$. Then $Z_{i}\perp Z_{j}\mid(Z_{k})_{k\neq
i,j}\iff\Omega_{ij}=0$.

The Graphical Lasso (GLasso) is an established approach to structure
estimation, which involves adding an $\ell_{1}$ penalty to a Gaussian
likelihood to define a convex program. Formally, suppose we have a collection
of samples $z_{1},\dots,z_{N}$ from some zero-mean multivariate Gaussian with
empirical covariance matrix
$\mathbf{C}=\frac{1}{N}\sum_{n=1}^{N}z_{n}z_{n}^{\top}\in\mathbb{R}^{\bar{p}\times\bar{p}}$
(in practice, we apply this to mean-centred samples). Then our estimator for
the precision matrix is defined by

$\hat{\Omega}\in\operatorname{argmax}_{\tilde{\Omega}\succeq\mathbf{0}}\\{\log\det\tilde{\Omega}-\operatorname{trace}(\mathbf{C}\tilde{\Omega})-\lambda\rho_{1}(\tilde{\Omega})\\}\;,$
(13)

where $\lambda\in(0,\infty)$ is a penalty parameter, and $\rho_{1}(\Omega)$ is
the $\ell_{{1}}$-norm of the off-diagonal entries of $\Omega$. This objective
was originally proposed in Yuan and Lin (2007) but soon prompted a wide
variety extensions and analysis; see Hastie et al. (2015) for an overview and
further references.

The GLasso has good theoretical behaviour under the assumption that the true
precision matrix $\Omega$ is sparse. Indeed, strong theoretical guarantees can
be proven even for non-Gaussian data by interpreting the objective Equation 13
as an $\ell_{{1}}$-penalised Bregman divergence (Ravikumar et al., 2011).
However, beyond certain special cases (Loh and Wainwright, 2013), there is in
general no reason to assume that non-Gaussian data will have sparse precision
matrix, even under a sparse graphical model. Nevertheless, one may expect
approximately Gaussian data to be approximately sparse and therefore for the
GLasso to perform well on a variety of real datasets. Indeed, the GLasso is
now widely used across a variety of scientific disciplines, as testified by
the thousands of citations accrued by early papers by Yuan and Lin (2007);
Friedman et al. (2008).

###  4.2 Definition of our estimator

To apply the Glasso to CCA, we take $\bar{p}=p+q$ and define random variables
$Z\in\mathbb{R}^{\bar{p}}$ by concatenation of the variables $(X,Y)$; then the
covariance of $Z$ is precisely $\Sigma$ as defined in (1). Our key observation
is that $\Omega$ determines $\Sigma$ and therefore also determines the
canonical decomposition. We propose first estimating $\Omega$ by
$\hat{\Omega}$ using a graphical lasso algorithm, and then computing canonical
decomposition corresponding to this estimate. For clarity we present this in
pseudo-code in Algorithm 1. We note that there are various different GLasso
implementations available, and any of them could be used. Furthermore, our
plug-in procedure could be modified to accept alternative sparse precision
estimators (see Cai et al., 2011, and extensions), as well as more complex
estimators from sparse plus low-rank decompositions or latent variable models
(e.g. Chandrasekaran et al., 2012).

Algorithm 1 Graphical CCA (gCCA)

procedure Graphical-CCA($\mathbf{X},\mathbf{Y}$)

$\mathbf{C}\leftarrow\texttt{sample-covariance}(\mathbf{X},\mathbf{Y})$
$\triangleright$ As defined by (5)

$\hat{\Omega}\leftarrow\texttt{glasso}(\mathbf{C})$ $\triangleright$ Use any
GLasso variant

$\hat{\Sigma}\leftarrow\hat{\Omega}^{-1}$

$\hat{U},\hat{V},\hat{R}\leftarrow\texttt{canonical-
decomposition}(\hat{\Sigma})$ $\triangleright$ Use any numerical method

return $\hat{U},\hat{V},\hat{R}$

end procedure

###  4.3 Theoretical guarantees

Theoretical guarantees can be obtained by working through Algorithm 1 and
showing that there is ‘good’ estimation at each step. We now outline the key
ideas of this argument, present the main result of Proposition 3, and leave
full mathematical details and further discussion to Appendix E.

When working through Algorithm 1, we were able to leverage some sophisticated
results from the mathematical statistics literature. To show the precision
estimate is good, we can use the theory for the graphical lasso from Ravikumar
et al. (2011), as we explain later in this section. To show this implies good
CCA estimates, we work via the SVD form of CCA from Section 2.1, and the
target matrix $T$ for SVD from Equation 3. By using a variant of the Davis-
Kahan theorem (Yu et al., 2015), one can prove that a good estimate $\hat{T}$
of the true SVD target $T$ will lead to good estimates of its singular value
decomposition.

The main technical work remaining was to show that a good precision estimate
gives a good estimate of $\hat{T}$. We formalised this by showing that the
function $\hat{\Omega}\mapsto\hat{T}$ is uniformly Lipschitz, as stated in
Proposition 2 and proved in Section E.3. The proof amounts to matrix analysis
and tracking Lipschitz constants through each step of the formula that
expresses $\hat{T}$ in terms of $\hat{\Omega}$. This only required the
following simple assumption on the minimal and maximal eigenvalues of the
within-view covariance matrices, which is natural in the context of CCA, and
used in previous analysis of high-dimensional CCA (Gao et al., 2016).

######  Assumption A (Well-conditioned within-view variances).

There exist constants $M,m\in(0,\infty)$ such that

$\lVert\Sigma_{xx}\rVert_{\textrm{op}},\lVert\Sigma_{yy}\rVert_{\textrm{op}}\leq
M,\quad\lVert\Sigma_{xx}^{-1}\rVert_{\textrm{op}}\;,\lVert\Sigma_{yy}^{-1}\rVert_{\textrm{op}}\leq
m\;.$

######  Proposition 2 (Lipschitz plug-in).

Assume Assumption A$(M,m)$. Then there exists a $C\in(0,\infty)$ only
depending on $M,m,\rho_{1}$ such that whenever
$\lVert\Omega-\hat{\Omega}\rVert_{\textrm{op}}\leq\frac{1}{4M}$ we have

$\lVert\hat{T}-T\rVert_{\textrm{op}}\leq
C\lVert\hat{\Omega}-\Omega\rVert_{\textrm{op}}\;.$ (14)

Before we can state our main result, we need to describe the results of
Ravikumar et al. (2011), which give finite sample bounds on the error of
$\hat{\Omega}$. It is most convenient for us to use their operator-norm
analysis. Their rates critically depend on the maximal node degree and number
of edges of the graph $G$, but also depend on the following structural
quantities determined by the true covariance matrix $\Sigma$, or equivalently
its inverse $\Omega$, the precision matrix.

######  Assumption B (Structural assumptions for GLasso).

Assume that:

  * •

the full covariance matrix $\Sigma={\Omega}^{-1}$ has
$\ell_{\infty}$-operator-norm $\kappa_{\Sigma}$

  * •

the log-det objective is $\kappa_{\Gamma}$-strictly convex when restricted to
the non-zero entries of $\Omega$

  * •

mutual incoherence: no zero interaction is too well aligned with the true
interactions, quantified by a constant $\alpha$.

######  Proposition 3 (Graphical CCA guarantee).

Consider observing independent samples $z_{n}=(x_{n},y_{n})_{n=1}^{N}$ from a
distribution with sparse joint precision matrix $\Omega$. Let $d$ be the
maximal node degree of the corresponding interaction graph, and $s$ be the
total number of edges in the graph. Assume that the standardised random
variables in both views are sub-Gaussian with some common parameter. Assume
Assumptions A$\,(M,m)$ and B$\,(\kappa_{\Gamma},\kappa_{\Sigma},\alpha)$ hold.
Then there exists some constant $C$ only depending on these structural
quantities $(m,M,\kappa_{\Gamma},\kappa_{\Sigma},\alpha)$ and the top
canonical correlation $\rho_{1}$, such that for sample size $N\gtrsim
d^{2}\log p$ and appropriate choice of penalty parameter $\lambda$ we have

$\|\sin\Theta(X^{\top}\hat{U}_{K},X^{\top}U_{K})\|_{\mathrm{F}}\leq
C\frac{K^{1/2}}{\rho_{K}^{2}-\rho_{K+1}^{2}}\sqrt{\frac{\min\\{s+\bar{p},d^{2}\\}\log\bar{p}}{N}}$
(15)

with high probability.

We now note certain properties of this bound. Firstly, the bound considers
subspace error, quantified by $\sin\Theta$ distance, in variate space. One
could use Assumption A to propagate this to an error bound in weight-space,
but this would add an extra factor of $m$ to the constant, and $m$ may well be
large. Secondly, the rate $\sqrt{\frac{\log p}{N}}$ is standard for
$\ell_{{1}}$ regularised problems, is the same order as the minimax bound from
Gao et al. (2016) for sparse-weight CCA, and we conjecture that it is minimax
in this sparse-precision regime. Thirdly, the strength of the result from
Ravikumar et al. (2011) is really in the regime when $d^{2}\ll s+p$, and
bounds error in high-dimensional regimes $p/n\to\infty$ provided that
$d^{2}\log\bar{p}/n\to 0$. Also note that the eigen-gap
$\rho_{K}^{2}-\rho_{K+1}^{2}$ inherited from the Davis-Kahan bound means that
the bound becomes vacuous as the eigengap tends to zero — we should not expect
good estimation if there is small separation of correlation signal; similar
eigen-gaps are present in all previous bounds for CCA, e.g., Gao et al.
(2016); Ma and Li (2020).

Finally, we note that this general proof idea could give a number of
alternative bounds, by using different matrix perturbation results, taking
slightly different theorems from Ravikumar et al. (2011), or indeed using a
different regularised precision matrix estimator, such as the latent variable
GLasso (Chandrasekaran et al., 2012).

###  4.4 Relationship to existing regularised CCA methods

There are certain connections between sparse precision matrices and sparse
canonical directions. In fact, canonical directions can inherit block-wise
sparsity from block-wise sparsity of the precision matrix; we now state this
formally, and give a proof in Section E.1. We also give a second, alternative
proof using partial correlations in Section E.2 to provide complementary
intuition.

######  Proposition 4 (Sparse directions from sparse precision).

Let $A\subset[p]$ correspond to a subset of the $X$ variables. Suppose that
$(\Omega_{xy})_{aj}=0$ for all $a\in A,j\in[q]$. If $\rho_{k}>0$ then any
$k^{{}^{\text{th}}}$ canonical variate can be written as $X^{\top}u_{k}$ where
$\operatorname{supp}(u_{k})\subset[p]\setminus A$.

Based on this relationship, we might well expect certain similarities between
the behaviour of gCCA to sCCA because both will tend to enforce some sparsity
of the estimated directions. However, this is only a loose heuristic because
there are many cases where Proposition 4 is not relevant. On the one hand, it
is easy to construct a sparse precision matrix such that each variable in view
1 was related to a variable in view 2 and vice versa; then Proposition 4 would
be vacuous. On the other hand, is also easy to construct covariance matrices
which have very sparse directions, but have non-sparse precision matrices.

There are also close parallels between gCCA and rCCA. In particular, both can
be seen as applying population CCA to some regularised covariance estimate.
Moreover, each of these regularisations is specified by a single scalar tuning
parameter and enforces that the covariance estimate is strictly positive
definite. Correspondingly, one may expect close behaviour between gCCA and
rCCA, especially for small tuning parameters.

These heuristic parallels are reflected in our empirical results: gCCA often
appears to perform regularisation ‘somewhere between’ sCCA and rCCA. In the
remaining sections we therefore plot gCCA between sCCA and rCCA for ease of
comparison.

###  4.5 Notes on implementation and application

To implement Algorithm 1 we use the excellent Graphical Lasso implementation
in the gglasso python package (Schaipp et al., 2021). The computational
complexity of the GLasso is $\mathcal{O}((p+q)^{3})$ so becomes infeasible for
very high-dimensional problems. However, because the algorithm requires a
single convex optimisation problem followed by a SVD operation, rather than a
sequence of successive problems, it is faster in practice than the sCCA and
sPLS approaches for problems of up to a few hundred dimensions.

##  5 Evaluating CCA Estimates

What does it mean for one CCA estimate to be ‘better’ than another? Answering
this question is crucial for a meaningful comparison of different CCA
algorithms on simulated data; it is also crucial for model selection on real
data.

A good CCA estimate will capture a lot of correlation signal and also be close
to the underlying population quantity; these two ideas both lead to evaluation
approaches. We present approaches based on correlation in Section 5.1 and
approaches based on estimation in Section 5.2. In each case we first suggest
‘oracle’ quantities to be used when the population distribution is known, and
relevant for evaluation on synthetic data; then we suggest ‘empirical’
quantities for real data situations when the population distribution is
unknown, and which are based on sample splitting. Additionally, in each case
we both suggest approaches which work with canonical pairs successively and
approaches that work with subspaces.

The sample splitting approaches require additional notation. We will always be
in the setting of $V$-fold cross-validation with folds $\nu=1,\dots,V$. We
write the data-matrices restricted to the $\nu^{\text{th}}$ fold as
$\mathbf{X}^{[\nu]},\mathbf{Y}^{[\nu]}$ (validation set) and the remaining
data in the matrices $\mathbf{X}^{[-\nu]},\mathbf{Y}^{[-\nu]}$ (training set).
Given an algorithm which obtains estimates
$\hat{U}(\mathbf{X},\mathbf{Y}),\hat{V}(\mathbf{X},\mathbf{Y})$, we introduce
the shorthand
$\hat{U}^{[-\nu]}=\hat{U}(\mathbf{X}^{[-\nu]},\mathbf{Y}^{[-\nu]}),\hat{V}^{[-\nu]}=\hat{V}(\mathbf{X}^{[-\nu]},\mathbf{Y}^{[-\nu]})$
and corresponding lower case quantities for direction estimates.

###  5.1 Correlation captured

####  5.1.1 Oracle successive correlations

One obvious way to evaluate a single pair of estimated canonical directions is
through the correlations between projections of the respective random
variables onto these estimated directions (this is after all the optimisation
objective in Equation 2). We refer to these as oracle correlations because
computing them requires knowledge of the true covariance matrix, which is not
available on real data:

$\rho^{\text{oracle}}(\hat{u},\hat{v})=\operatorname{Corr}(\hat{u}^{\top}X,\hat{v}^{\top}Y)=\frac{\hat{u}^{\top}\Sigma_{xy}\hat{v}}{\sqrt{\rule{0.0pt}{9.04166pt}\hat{u}\Sigma_{xx}\hat{u}}\,\sqrt{\hat{v}^{\top}\Sigma_{yy}\hat{v}}}\;.$
(16)

If one has $K$ estimated pairs, one could consider the whole vector of
successive oracle correlations
$\left(\rho^{\text{oracle}}(\hat{u}_{k},\hat{v}_{k})\right)_{k=1}^{K}$ but
usually desires a single scalar metric to evaluate the estimates. We propose a
general family of such scalar metrics of the form

$\displaystyle\texttt{succ-cc-
agg}(f;\hat{U},\hat{V})=f\left(\left(\rho^{\text{oracle}}(\hat{u}_{k},\hat{v}_{k})\right)_{k=1}^{K}\right)\;,$
(17)

where $f:\mathbb{R}^{*}\to\mathbb{R}$ is some coordinate-wise increasing
function. Choosing $f=\|\cdot\|_{1}$ corresponds to sums of correlations;
choosing $f=\|\cdot\|_{2}^{2}$ gives sums of squared correlations, which are
particularly natural in light of Section 7.3. Other choices of $f$ may also be
interesting. For example the function
$f:\rho\mapsto-\frac{1}{2}\sum_{k}\log(1-\rho_{k}^{2})$ corresponds to mutual
information in a Gaussian setting (see Section F.3.2) and is therefore also
reasonable to consider. These choices of $f$ are justified by the following
result, which we state precisely and prove as Proposition 27 in Appendix D.

######  Proposition 5 (Valid aggregation functions).

Let $\phi:\mathbb{R}\to\mathbb{R}$ be a convex function such that
$\phi(x)\geq\phi(-x)~{}\forall x>0$. Then, top-$K$ CCA subspaces can be
characterised by the program

$\displaystyle\underset{\begin{subarray}{c}\\{u_{1},\dots,u_{K}\\},\\{v_{1},\dots,v_{K}\\}\\\
u_{k}^{\top}\Sigma_{xx}u_{l}=\delta_{kl},v_{k}^{\top}\Sigma_{yy}v_{l}=\delta_{kl}:\,\forall
k,l\in\\{1,\dots,K\\}\end{subarray}}{\operatorname{maximise}}{\vphantom{\Bigg{(}\Bigg{)}}\>\>\sum_{k=1}^{K}\phi\left((\rho^{\text{oracle}}(u_{k},v_{k})\right)}\;.$

This naive approach of combining correlations of successive pairs has one
major drawback though — and that is when the estimated directions are not
suitably orthogonal. In this case we may end up double-counting some of the
correlation signal. One way to resolve this is to consider instead subspace
correlations.

####  5.1.2 Oracle subspace correlations

Because the canonical correlations are unique, we can make the following
definition.

######  Definition 1.

Write $\mathcal{R}^{\text{CCA}}_{[{K}]}(X,Y)\in\mathbb{R}^{K}$ to be the
vector containing the top-$K$ canonical correlations of the pair of random
variables $(X,Y)$ as defined by (2). Write
$\hat{\mathcal{R}}^{\text{CCA}}_{[{K}]}(\mathbf{X},\mathbf{Y})\in\mathbb{R}^{K}$
to be vector containing the top-$K$ canonical correlations of the sample CCA
problem corresponding to $(\mathbf{X},\mathbf{Y})$ as defined by (8).

We suggest analysing the vector of canonical correlations corresponding to the
pair of $K$-dimensional random variables
$\hat{U}_{K}^{\top}X,\hat{V}_{K}^{\top}Y$, i.e.
$\mathcal{R}^{\text{CCA}}_{[{K}]}(\hat{U}_{K}^{\top}X,\hat{V}_{K}^{\top}Y)$
with the notation of Definition 1. Note from the discussion in Section B that
this only depends on the subspace spanned by the random variables
$\hat{U}_{K}^{\top}X,\hat{V}_{K}^{\top}Y$; therefore it only depends on the
subspace spanned by the estimates $\hat{U}_{K},\hat{V}_{K}$; these are the
primary objects of interest if we are using CCA for dimension reduction.
Similarly to in (17), we propose a family of metrics of the form

$\displaystyle\texttt{subsp-cc-
agg}_{K}(f;\hat{U},\hat{V})=f\left(\mathcal{R}^{\text{CCA}}_{[{K}]}(\hat{U}_{K}^{\top}X,\hat{V}_{K}^{\top}Y)\right)\;.$
(18)

####  5.1.3 Empirical correlations via cross validation

Each of the metrics above can be approximated by cross validation (CV). The CV
analogue of (16) is

$\displaystyle\rho^{\text{CV}}(\hat{u},\hat{v})=\frac{1}{V}\sum_{\nu=1}^{V}\widehat{\operatorname{Corr}}\left(\mathbf{X}^{[\nu]}\hat{u}^{[-\nu]},\mathbf{Y}^{[\nu]}\hat{v}^{[-\nu]}\right)\;.$

One CV analogue of (17) is

$\displaystyle\texttt{cv-succ-cc-
agg}_{K}(f;\hat{U},\hat{V})\vcentcolon=\frac{1}{V}\sum_{\nu=1}^{V}f\left(\left(\widehat{\operatorname{Corr}}\left(\mathbf{X}^{[\nu]}\hat{u}^{[-\nu]}_{k},\mathbf{Y}^{[\nu]}\hat{v}^{[-\nu]}_{k}\right)\right)_{k=1}^{K}\right)\;.$
(19)

Finally, a CV analogue of the subspace correlation notion from (18) is

$\displaystyle\texttt{cv-subsp-cc-
agg}_{K}(f;\hat{U},\hat{V})\vcentcolon=\frac{1}{V}\sum_{\nu=1}^{V}f\left(\hat{\mathcal{R}}^{\text{CCA}}_{[{K}]}\left(\mathbf{X}^{[\nu]}\hat{U}_{K}^{[-\nu]},\mathbf{Y}^{[\nu]}\hat{V}_{K}^{[-\nu]}\right)\right)\;.$
(20)

###  5.2 Estimation Error

Table 1: Estimation metrics Oracle version | Empirical version
---|---
Quantity | Definition | Quantity | Definition
wt-uk | $\sin^{2}\Theta(u_{k},\hat{u}_{k})$ | wt-uk-cv | $\binom{V}{2}^{-1}\sum_{\nu<\nu^{\prime}}\sin^{2}\Theta(\hat{u}_{k}^{[-\nu]},\hat{u}_{k}^{[-\nu^{\prime}]})$
vt-uk | $\sin^{2}\Theta(\Sigma_{xx}^{\nicefrac{{1}}{{2}}}u_{k},\Sigma_{xx}^{\nicefrac{{1}}{{2}}}\hat{u}_{k})$ | vt-uk-cv | $\binom{V}{2}^{-1}\sum_{\nu<\nu^{\prime}}\sin^{2}\Theta(\mathbf{X}\hat{u}_{k}^{[-\nu]},\mathbf{X}\hat{u}_{k}^{[-\nu^{\prime}]})$
wt-Uk | $\sin^{2}\Theta(U_{k},\hat{U}_{k})$ | wt-Uk-cv | $\binom{V}{2}^{-1}\sum_{\nu<\nu^{\prime}}\sin^{2}\Theta(\hat{U}_{k}^{[-\nu]},\hat{U}_{k}^{[-\nu^{\prime}]})$
vt-Uk | $\sin^{2}\Theta(\Sigma_{xx}^{1/2}U_{k},\Sigma_{xx}^{1/2}\hat{U}_{k})$ | vt-Uk-cv | $\binom{V}{2}^{-1}\sum_{\nu<\nu^{\prime}}\sin^{2}\Theta(\mathbf{X}\hat{U}_{k}^{[-\nu]},\mathbf{X}\hat{U}_{k}^{[-\nu^{\prime}]})$

####  5.2.1 Oracle

As in Section 5.1 there are many reasonable choices of metrics for evaluating
estimation. Indeed, it is again reasonable to consider estimates successively
or with a subspace-based approach. But now there are at least three different
aspects of the estimates of interest: variates, weights or loadings.

For simplicity we will only consider squared $\sin\Theta$ distances. These
have a natural sense of scale, with a maximum value of $K$ between two $K$
dimensional subspaces, this linear structure helps makes the graphs more
readable.

To apply in variate space we simply pre-multiply directions by
$\Sigma_{xx}^{\nicefrac{{1}}{{2}}}$, and to apply in loading space we pre-
multiply by $\Sigma_{xx}$ (here, we focus on estimation of variates). Recall
that in this section we are primarily concerned with numerical experiments
where $\Sigma_{xx}$ is indeed known, as are the true directions and so these
distances are easily computable. We define all relevant quantities, and their
corresponding shorthand, in the left hand column of Table 1.

####  5.2.2 Empirical: cross-validated instability

It becomes difficult to estimate these quantities empirically, but one can
define related notions of instability. Though these do not give good estimates
of the oracle error, they give approximate lower bounds and may tell us that a
CCA estimator has too high variance to be accurate.

The right hand column of Table 1 defines various notions of instability, and
gives corresponding shorthand. We propose looking at average difference
between pairs of estimates trained on different training folds. Notice that
the definitions given for instability in variate space (vt-uk-cv and vt-Uk-cv)
use the full dataset $\mathbf{X}$ to obtain sample canonical variates; though
this may raise concern about ‘re-using’ data, it has lower variance than
conceivable alternatives based on further sample splitting and gave sensible
empirical results.

##  6 Synthetic Data

Using synthetic data is essential for developing statistical methodology but
conclusions rely heavily on the choice of model and metrics used. Previous
work has used simplistic synthetic models, with very clear signal. By
contrast, in this section we will only consider synthetic data that we claim
to be biologically plausible, which is obtained via parametric bootstrap
resampling of a real dataset. We use mean-centred multivariate Gaussian
distributions in line with previous work, which can be justified by central
limit theorem arguments and does not appear restrictive. We keep this section
short, with minimal content to support our claims in Section 1 and to provide
context for our real data suggestions in Sections 7 and 8. Further synthetic
data experiments based on toy data-generating mechanisms, as well as different
parametric bootstrap models are included in Appendices G and H respectively.

Figure 1: Oracle metrics on the parametric-bootstrapped Microbiome dataset,
see Tables 1 and 2 for a glossary of the legends.

Figure 2: Correlation and instability metrics over the regularisation path for
the four methods on the parametric bootstrapped Microbiome dataset; error bars
for the aggregated quantities; oracle values are dotted for comparison. Top
row: sums of squared correlations, as defined by r2sk-cv and r2sk in Table 2
for $k=1,3,6$. Middle row: subspace instability in weight space, as defined by
wt-Uk-cv and wt-Uk in Table 1 for $k=1,3$. Bottom row: subspace instability in
variate space, as defined by vt-Uk-cvand vt-Uk in Table 1 for $k=1,3$.

In particular, this section considers a single synthetic version of the
Microbiome dataset described in Section 8.1. To construct the dataset, we
first obtained an estimated precision matrix using the GLasso algorithm; we
used the tuning parameter $\alpha=0.0059$, which maximised the sum of the
top-5 squared test correlations (r2s5-cv from Table 2). We then generated a
set of $n=500$ samples from a multivariate normal distribution with this
precision matrix. Though we use a single set of pseudo-random samples using an
arbitrary random seed, our observations are robust to the choice of this
random seed. We will see that many behaviours observed in this section closely
mirror those in Section 8 which suggests this synthetic dataset is realistic.

###  6.1 Metrics behave differently

Table 2: Correlation metrics Oracle version (error) | Empirical version (instability)
---|---
Quantity | Definition | Refs. | Quantity | Definition | Refs.
r2sk | $\texttt{succ-cc-agg}_{K}(\lVert\cdot\rVert_{2}^{2};\hat{U},\hat{V})$ | (17) | r2sk-cv | $\texttt{cv-succ-cc-agg}_{K}(\lVert\cdot\rVert_{2}^{2};\hat{U},\hat{V})$ | (19)
R2sk | $\texttt{subsp-cc-agg}_{K}(\lVert\cdot\rVert_{2}^{2};\hat{U},\hat{V})$ | (18) | R2sk-cv | $\texttt{cv-subsp-cc-agg}_{K}(\lVert\cdot\rVert_{2}^{2};\hat{U},\hat{V})$ | (20)

For our single bootstrapped dataset we show how the ‘oracle’ behaviour of our
different algorithms varies with tuning parameter. To evaluate correlation
captured we look at sums of squared oracle correlations. To evaluate
estimation error, we use $\sin^{2}\Theta$ distances, for both successive
directions and for subspaces, and both in weight space and in variate space;
these distances have convenient interpretations in terms of fraction of signal
shared. We summarise these different metrics in Table 2, with the text labels
that will be used in the legends of our plots (estimation metrics are in Table
1 for reference).

Figure 1 supports the following observations:

  * •

many directions containing signal can be recovered: the top row values are
high — for example the gCCA has squared fifth correlation of 0.6,
corresponding to correlation of 0.8, which is very significant.

  * •

variates are easier to estimate than weights: the losses for wt-uk are much
larger than for vt-uk and the losses for wt-Uk are much larger than the losses
for wt-Vk.

  * •

subspaces are easier to estimate than directions: the losses for wt-uk are
much larger than for wt-Uk and the losses for vt-uk are much larger than for
vt-Uk.

  * •

optimal tuning parameters for different metrics are often similar: for example
the gCCA optimal parameters are almost all between $5\times 10^{-3}$ and
$10^{-2}$; moreover all parameters in this range are near-optimal for most
different metrics

  * •

there is often no obvious best model or set of hyperparameters: in this
example, it seems like the gCCA with parameter $6\times 10^{-3}$ would be a
good overall choice; however sCCA has comparable correlation captured, and
rCCA has comparable variate loss.

  * •

sPLS does not perform CCA: it captures meaningfully lower correlation than CCA
methods with conventional orthogonality constraints, and produces weights with
near-maximal loss. However, it does recover some of the shared variate
subspace, because PLS and CCA subspaces can often be similar.

Because there is no obvious best model or set of hyperparameters, it may not
be sensible to perform strict model selection. Instead, it may make more sense
to consider a variety of different models with near-optimal hyperparameters,
and compare the corresponding estimates. This will give a rough measure of how
robust an observation drawn from these estimates is to the choice of
regularisation; such robustness is certainly desired if such an observation is
to merit further investigation. We provide more details of a possible approach
in Section 8.

###  6.2 Cross validation is effective

In Figure 2 we compare CV notions of various criteria with their oracle
counterparts. In each plot the CV quantities are plotted with solid lines and
error bars, while the oracle quantities are plotted with fainter dotted lines;
in Section 8 we will use the same style of plots for CV quantities, but
clearly cannot access the oracle quantities.

The important observation is that all the CV correlation criteria closely
match their oracle counterparts. Moreover, the penalty parameters optimising
the CV quantities are usually the same as the parameters optimising the oracle
quantities.

The picture is much more complicated for the notions of stability. As one
might expect, the estimates tend to be more stable for higher penalty
parameter. However, the variation in stability is not monotonic, and for sCCA
and gCCA there appear local minima of empirical instability, though these do
not align well with minima of oracle estimation error, or maxima of
correlation signals. Also, note that the CV instability versions are nearly
always below oracle estimation error counter parts — confirming our
expectation that instability does give a lower bound on estimation error.

##  7 Comparison and interpretation

We now move on to the questions of key practical importance. Given estimated
CCA directions, how should we interpret them? Moreover, if there are multiple
CCA estimates that are plausible how can we compare them? Following the
results of Section 6 and motivation of Section 1, we take a subspace-and-
variate centric approach. Our main goal is to introduce biplots for
interpretation of CCA estimates. We will also define a notion of overlap
matrices to visually compare the overlapping signal of two CCA estimates.
First, it is important to discuss registration.

###  7.1 Registration

One practical problem when comparing CCA estimates is that they are not well
aligned. Sometimes two estimates appear very different, but actually span very
similar subspaces. One way to deal with this is registration.

It is sensible to perform this registration ‘in variate space’, because the
variates are more stable than the weights; indeed, registering in ‘weight
space’ might fail to notice that two sets of very different weights can lead
to very similar variates. However, beyond this, there are many different types
of registration that might be reasonable. We now outline a reasonable family
of techniques for registration, which are closely related, but may still yield
complementary insights.

Suppose we have some reference estimates $\hat{U}_{K}\in\mathbb{R}^{p\times
K}$ and some new estimate for comparison
$\hat{U}^{\prime}_{K^{\prime}}\in\mathbb{R}^{p\times K^{\prime}}$ for some
$K^{\prime}\geq K$. We wish to find a $K$-dimensional subspace of columns of
$\hat{U}^{\prime}_{K^{\prime}}$ that most closely correspond to $\hat{U}_{K}$.
Let our data be $\mathbf{X}\in\mathbb{R}^{n\times p}$. Then our registered
direction estimates are $\hat{U}^{\prime}_{K^{\prime}}M^{*}$ where

$\displaystyle
M^{*}\in\operatorname*{arg\,min}_{M\in\mathcal{M}(K^{\prime}\times
K)}\|\mathbf{X}\hat{U}_{K}-\mathbf{X}\hat{U}^{\prime}_{K^{\prime}}M\|_{F}^{2}$

for some suitable set of matrices $\mathcal{M}(K^{\prime}\times
K)\subset\mathbb{R}^{K^{\prime}\times K}$.

By picking appropriate sets $\mathcal{M}$ we can register up to sign flips,
permutations (and sign flips), orthogonal transformations, or indeed up to
arbitrary linear transformations. In Section C.1 we show how to implement
these types of registration and discuss their relationship to subspace
distance notions. We have found each of these types of registration insightful
in different contexts, as we will see in the remaining sections.

###  7.2 Overlap matrices

We now provide background material for a visualisation tool which we will
demonstrate in Section 8.3.4.

######  Definition 2 (Overlap matrices).

Given two matrices $\textbf{Z},\textbf{W}\in\mathbb{R}^{N\times K}$ we call
$\textbf{Z}^{\top}\textbf{W}\in\mathbb{R}^{K\times K}$ the overlap matrix
between $\textbf{Z},\textbf{W}$; we define the squared overlap matrix to be
the elementwise-squared overlap matrix.

The key motivation for this definition is that if $\textbf{Z},\textbf{W}$ have
orthonormal columns $(\textbf{Z}_{k})_{k=1}^{K},(\textbf{W}_{k})_{k=1}^{K}$,
then from Section A.3, the squared cosine similarity between these subspaces
is

$\displaystyle\cos^{2}\Theta(\textbf{Z},\textbf{W})=\lVert\textbf{Z}^{\top}\textbf{W}\rVert_{F}^{2}=\sum_{k=1}^{K}\sum_{l=1}^{K}\langle\textbf{Z}_{k},\textbf{W}_{l}\rangle^{2}\;,$
(21)

which is precisely the sum of squared entries of the overlap matrix
$\textbf{Z}^{\top}\textbf{W}$. Moreover, by the same argument as in Equation
21 we see that sum of entries of any squared overlap submatrix correspond to
squared cosine similarity of the corresponding subspaces. The (squared)
overlap matrices therefore give an interpretable graphical display of how the
signal in Z overlaps with the signal in W.

In Section 8.3.4 we apply this to two matrices of estimated canonical
variates, i.e.

$\displaystyle\textbf{Z}=\mathbf{X}\hat{U}_{K},\quad\textbf{W}=\mathbf{X}\hat{U^{\prime}}_{K};\quad\text{so
}\textbf{Z}^{\top}\textbf{W}=\widehat{\operatorname{Cov}}(\mathbf{X}\hat{U}_{K},\mathbf{X}\hat{U^{\prime}}_{K})\;.$
(22)

If our estimated canonical variates are orthonormal (as they should be) then
the cosine similarity interpretation applies, but even if not, the overlap
matrix may still be informative. Moreover, by applying Gram–Schmidt one can
plot the overlap matrix corresponding to successive subspaces; by comparing to
the original overlap matrix this can give a clear visualisation of the degree
of non-orthogonality of the data. Note that we can also apply this in the
oracle settingby instead considering the matrix

$\displaystyle\operatorname{Cov}(U_{K}^{\top}X,\hat{U}_{K}^{\prime\,T}X)=(\Sigma_{xx}^{\nicefrac{{1}}{{2}}}U_{K})^{\top}(\Sigma_{xx}^{\nicefrac{{1}}{{2}}}U^{\prime}_{K})\;.$

Note it is also straightforward to consider analogous ‘validation versions’
for a sample split by taking

$\displaystyle\textbf{Z}=\mathbf{X}^{[\nu]}\hat{U}_{K}^{[-\nu]}\;,\quad\textbf{W}=\mathbf{X}^{[\nu]}\hat{U^{\prime}}_{K}^{[-\nu]}\;.$

###  7.3 Biplots

Biplots are a powerful tool for visualising canonical correlation analysis,
originally proposed in ter Braak (1990); these plot the variables from each of
the two views of data in a shared low-dimensional space. Though used in a wide
variety of contexts and presented in tutorials such as Uurtio et al. (2017),
they appear under-utilised in scientific applications of CCA. The biplot was
adapted in ter Braak (1990) from an earlier notion of biplot for PCA (Gabriel,
1971; Jolliffe, 1986); we warn the reader that the word biplot is often used
loosely in the applied literature to refer to a variety of plots based on
similar ideas.

The original biplot proposed in ter Braak (1990) was for sample CCA and breaks
down in the high-dimensional sample setting. We instead present a population
notion of biplot, whose sample version is equivalent to that of ter Braak
(1990) in low dimensions but which remains meaningful in high dimensions.
Though we define this using the $X$-canonical variates, one could equally
consider the $Y$-canonical variates; this is also true with all the related
formulations in the following sections.

####  7.3.1 Population biplots

######  Definition 3 (Population Correlation Biplot).

Let $\\{u_{k},v_{k}\\}_{k=1}^{K}$ be successive canonical directions. The
$X$-view biplot plots the $i^{\text{th}}$ variable in view 1 and
$j^{\text{th}}$ variable in view 2 at positions $\phi(X_{i}),\phi(Y_{j})$
defined by

$\phi(X_{i})\vcentcolon=\left(\operatorname{Corr}(X_{i},u_{k}^{\top}X)\right)_{k=1}^{K}\;,\quad\phi(Y_{j})\vcentcolon=\left(\operatorname{Corr}(Y_{j},u_{k}^{\top}X)\right)_{k=1}^{K}\;,$
(23)

respectively.

The vectors
$\left(\operatorname{Corr}(X_{i},u_{k}^{\top}X)\right)_{i=1}^{p},\left(\operatorname{Corr}(Y_{j},u_{k}^{\top}X)\right)_{j=1}^{q}$
are called (the $k^{\text{th}}$ pair of) structure correlations. These can be
seen as normalised versions of the canonical loadings via

$\displaystyle(\Sigma_{xx}u_{k})_{i}=\operatorname{Cov}(X_{i},u_{k}^{\top}X)=\operatorname{Var}(X_{i})^{\nicefrac{{1}}{{2}}}\operatorname{Corr}(X_{i},u_{k}^{\top}X)\;.$

The biplot can therefore be seen as a visualisation of the canonical loading
vectors and interpreted via the latent variable formulation of CCA from
Section B.2.2.

######  Proposition 6.

Let $\phi$ be the biplot mapping from Definition 3.

  * •

The squared norm of a variable’s representation is the proportion of the
variable’s variance explained by the the first $K$ canonical variates, and so
is less than one. Writing $\xi_{k}=u_{k}^{\top}X$ for the $k^{\text{th}}$
canonical variate for the $X$-view, we have

$\displaystyle\lVert\phi(X_{i})\rVert_{2}^{2}=\sum_{k=1}^{K}\operatorname{Corr}(X_{i},\xi_{k})^{2}=\frac{\sum_{k=1}^{K}\operatorname{Cov}(X_{i},\xi_{k})^{2}}{\operatorname{Var}(X_{i})}\leq
1\;.$ (24)

  * •

Inner products of representations with large squared norm approximate
correlations between the corresponding variables: for any $i,i^{\prime}\in[p]$

$\displaystyle\left\lvert\operatorname{Corr}(X_{i},X_{i^{\prime}})-\langle\phi(X_{i}),\phi(X_{i^{\prime}})\rangle_{2}\right\rvert\leq\left(1-\lVert\phi(X_{i})\rVert^{2}\right)^{\nicefrac{{1}}{{2}}}\left(1-\lVert\phi(X_{i^{\prime}})\rVert^{2}\right)^{\nicefrac{{1}}{{2}}}\;,$
(25)

and the analogous statement also holds for the variables in the $Y$-view.

  * •

For inner products between representations in different views we have the
stronger bound

$\displaystyle\left\lvert\operatorname{Corr}(X_{i},Y_{j})-\langle\phi(X_{i}),\phi(Y_{j})\rangle_{2}\right\rvert\leq\rho_{K+1}\left(1-\lVert\phi(X_{i})\rVert^{2}\right)^{\nicefrac{{1}}{{2}}}\left(1-\lVert\phi(Y_{j})\rVert^{2}\right)^{\nicefrac{{1}}{{2}}}\;.$
(26)

###### Proof.

Equations 24 and 25 follow from the orthonormality of the $\xi_{k}$, whereas
Equation 26 also uses the fact that
$\operatorname{Cov}(\xi_{k},Y_{j})=\rho_{k}\operatorname{Cov}(v_{k}^{\top}Y,Y_{j})$
due to the canonical correlation structure. ∎

This final bullet point suggests that inner products between representations
of variables in different views may be ‘good’ approximations of the
correlation between those variables. In fact, ter Braak (1990) proves that
this approximation is optimal in a certain sense, using the Eckhart–Young
inequality; we present the result with associated discussion in Section C.2.2.

####  7.3.2 Sample biplots for regularised CCA, and their interpretation

Conveniently, Definition 3 has a natural sample equivalent.

######  Definition 4 (Sample correlation biplot).

Given an estimate $\hat{U}_{K}$, respectively plot the $i^{\text{th}}$
variable in view 1 and $j^{\text{th}}$ variable in view 2 with coordinates

$(\widehat{\operatorname{Corr}}(\mathbf{X}_{i},\mathbf{X}\hat{u}_{k}))_{k=1}^{K}\;,\quad(\widehat{\operatorname{Corr}}(\mathbf{Y}_{i},\mathbf{X}\hat{u}_{k}))_{k=1}^{K}\;.$

If one uses estimates from (unregularised) sample CCA then observations
Equations 24, 25 and 26 from the previous section carry through, with
quantities replaced by their sample counterparts. Unfortunately, as outlined
in Section 1, the break-down of sample CCA means that typically all sample
canonical correlations are one ($\hat{\rho}_{k}=1$ for each $k$) and therefore
Equation 26 is no stronger than Equation 25.

The hope is that using estimates from appropriate regularised CCA methods will
give good estimators for the true population biplots in a statistical sense.
This appears justified in certain high-dimensional situations. Indeed, for the
Microbiome dataset, the synthetic experiments display accurate estimation of
the canonical variates (see Figure 1) and therefore also loading vectors and
biplots. Additionally, in Section 8.3.5 we see that variates for the different
regularised CCA methods on the Microbiome dataset are extremely similar after
appropriate registration (see Figure 5); they are also stable to sub-sampling
the data (see the stability plots in Figure 3).

However, we observed empirically that the biplots may give meaningful
visualisations of the data, even when they are poor statistical estimators for
the population quantities. Indeed, for the ‘Breastdata’ dataset that we
discuss in Section 8.4, the biplots for sCCA and rCCA look very different, so
cannot both be good statistical estimates; yet both appear to recover
meaningful chromosomal information (see Figure 22).

######  Remark 3 (CV for biplots is difficult).

One could evaluate the statistical properties of biplots with cross-
validation. One might estimate the canonical directions with some training
subset of the data, and evaluate the empirical correlations using a held out
testing subset of data; however, this leads to large variance in the plots,
particularly with small sample sizes. We therefore suggest using the same
(full) dataset to construct the estimates and visualise them in the biplot,
but we cannot expect it to be a good estimate of the population biplot in
general.

##  8 Real Data

We now evaluate the various CCA methods on real data sets and propose a
pipeline of tools for EDA. We describe this pipeline in Section 8.3, applied
to our motivating microbiome dataset, which we first describe in Section 8.1.
We also applied the pipeline to two existing datasets, the Nutrimouse dataset
(available in the CCA R package (González and Déjean, 2021), and originally
from Martin et al. (2007)), and the BreastData dataset (available in the PMA R
package (Tibshirani, 2020) and originally from Chin et al. (2006)); these
datasets have far fewer samples, leading to rather different CCA behaviour,
see Section 8.4. We present only a minimal selection of plots to demonstrate
our main points and leave more extensive supporting plots to Appendix I.

###  8.1 Microbiome dataset

This whole study was motivated by a desire to perform EDA on data concerning
the human gut microbiome. We use an illustrative dataset obtained from the
Integrative Human Microbiome Project, concerning the microbiome’s role in
Inflammatory Bowel Disease (IBD) (Lloyd-Price et al., 2019). This dataset
contains counts of microbial metabolites (C0s) and genes (K0s) collected
logitudinally from patients with two different types of IBD diagnosis,
Ulcerative Colitis or Crohn’s Disease, as well as healthy control patients. We
filtered the raw counts data to ensure that there were no zero counts, then
applied the following steps common in compositional data analysis: first
convert raw counts into proportions, then apply a log transform, finally mean-
centre the columns. After this we had $p=711$ metabolites, $q=148$ proteins
and $n=458$ samples.

###  8.2 Background and setup

Our framing is very different to most previous work on high-dimensional CCA.
In the extreme case, some of these works only consider the top pair of weight
vectors, for a single penalty parameter chosen by CV, and essentially inspect
the entries. By contrast, we would like to find subspaces containing ‘most’ of
the correlation signal, and would like our observations to be indicative of
true population quantities. We therefore consider both estimated canonical
variates and directions, and ensure any observations we make are at a minimum
robust to sample splitting. In addition, we acknowledge that all our
regularised CCA models are likely to be mis-specified; we expect that
different sorts of regularisation may reveal complementary insights and
therefore care more about model comparison than model selection. However this
leads to the following difficult situation.

Suppose that a variety of different estimators have been fitted, with a
variety of penalty parameters, both to the original dataset and to different
CV folds. We shall refer to a combination of algorithm and penalty parameter
as an _estimator_ , as is the convention in the sklearn package. We therefore
have large family of direction estimates, indexed by algorithm, penalty, and
fold. We seek visualisation tools to help compare between this family of
estimates at a high-level, and decide which estimates to analyse in more
detail.

###  8.3 A framework for applying high-dimensional CCA

As in Section 5, it is sensible both to compare correlations captured by
different estimators and the differences between the estimates themselves. In
each case, there is a trade-off between giving a high-level comparison of a
large number of estimators, and richer, more-detailed comparisons of a smaller
number of estimators. High level summaries of CV correlation and instability,
and of distances between many estimates are described in Sections 8.3.1 and
8.3.3 respectively. Richer views of CV correlation captured and comparison
between a small number of estimates are described in Sections 8.3.2 and 8.3.4
respectively. Finally, we illustrate biplots in Section 8.3.5.

####  8.3.1 Correlation and instability along trajectories

Figure 3: Top: CV sums of correlations as function of regularisation path for
the four methods on the Microbiome dataset; error bars for the aggregated
quantities. Bottom: stability both in weight space and variate space along the
same trajectories.

Figure 3 provides a way to visualise correlation and instability at a very
high-level. It is analogous to Figure 2: the top row shows cross-validated
correlations, using both direction and subspace based approaches, while the
bottom row shows subspace instability, both in direction space and variate
space (see Section 5). This gives a visual comparison of all the different
estimators (algorithm, penalty pairs) of interest, and illustrates the
difficulty of regularising hard enough to encourage stability but while still
capturing the main signal.

The choices of $k$ of 1,3,5 are reasonable, but arbitrary, default dimensions
to consider; it might be worth considering different values of $k$ after
further analysis (particularly inspection of plots in Section 8.3.2).

Figure 3 does indeed look very similar to Figure 2. We see evidence of the
claims previously alluded to: many successive directions with significant
signal are recovered; variates are much more stable than weights; gCCA
captures slightly more correlation than the other CCA methods; sPLS captures
significantly less correlation and is relatively unstable. We therefore omit
sPLS from the remaining plots in this section to prevent distraction, and
defer further consideration of sPLS to Appendix I.

####  8.3.2 Decay of test correlations

Figure 4: CV correlations (colours) and average values (black) for successive
direction estimates using sCCA, gCCA, rCCA on microbiome dataset; in each case
r2s5-cv optimal penalty parameters were used.

Figure 4 shows how the test correlations captured by three estimators
decreases over successive components on the Microbiome dataset; the estimators
correspond to r2s5-cv-optimal penalties for sCCA, gCCA, rCCA respectively.

The main aim of the plot is to determine the number of components to consider.
An ‘elbow’ in the plot would indicate a sensible number of components.
Although it can only compare a limited number of estimators, this plot gives
much richer information than in Figure 3. Firstly, we can consider all
successive directions up to a large maximum value without cramping the
picture. Plotting individual folds rather than error bars gives clear visual
guide of statistical significance (all significantly non-zero). Color-coding
the folds gives even richer information that occasionally merits further
investigation, see further discussion in • ‣ Section I.2.

The main observations from Figure 4 are that a large number of directions
carry significant correlation signal, and that there are no clear elbows.
Indeed, for gCCA, there is very little decay of correlations between the 5th
and 10th components, and there seems to be a large number of successive
correlations around 0.6. Moreover, all of the individual test correlations are
well above zero; which suggests the signal is strongly significant. We also
observe that gCCA consistently records higher correlations than sCCA or rCCA.

####  8.3.3 Trajectory comparison via subspace distance

Figure 5: $\sin^{2}\Theta$ distances for top-3 variate subspace between pairs
of full sample estimates; axes refer to three values of the penalty parameter
for each of the three regularisation methods: sCCA, gCCA, and rCCA.

To get a high-level comparison of many different estimators, we suggest
plotting a trajectory comparison matrix, as in Figure 5. This is a symmetric
matrix whose axes refer to various values of the penalty parameter for each of
the regularisation methods (sCCA, gCCA, and rCCA), and may be thought of as
‘concatenated trajectories’; the entries display distances between the
corresponding pairs of estimators. Figure 5 shows such a matrix on the
Microbiome dataset using top-3 $\sin^{2}\Theta$ distance in variate space
(i.e. vt-U3 from Table 1). As in the previous Section 8.3.2, we chose the
central penalty parameters by CV (with objective r2s5-cv from Table 2) then
chose further parameters a factor of 3 larger and smaller.

The reassuring observation from Figure 5 is that the estimates are all fairly
similar. Indeed, $\sin^{2}\Theta$ distance between subspaces of dimension 3
can be at most 3; most values plotted are far smaller than this, and even the
maximum value of 1.1 corresponds to ‘about two thirds the same subspace’.
Moreover, we can obtain a submatrix with maximal value 0.38 by taking sCCA
0.0042, gCCA 0.002 and rCCA 0.33.

This main observation of ‘closeness’ is particularly interesting because it
does not hold for all metrics. Indeed, plotting top-3 $\sin^{2}\Theta$
distance in direction space gave far higher distances, see Section I.1. This
illustrates how ‘variate distance’ is a weaker notion in the context of CCA;
different forms of regularisation may recover similar canonical variates from
rather different sets of directions. In practice, considering a variety of
different metrics with multiple values of $k$ would help give richer view.

####  8.3.4 Overlap Matrices help visualise sharing of signal across
directions

Figure 6: Squared overlap matrices comparing the sCCA, gCCA, rCCA estimators
with maximal CV correlation (r2s5-cv); see Section 8.3.4 for full details.

Figure 6 shows squared overlap matrices (as in Equation 22 from Section 7.2)
on the Microbiome dataset; these have additional rows and columns with the
corresponding row and column sums. These matrices help give a much richer
comparison of different estimates, and how their signal overlaps across
different directions.

The squared overlap matrices in Figure 6 compare the sCCA, gCCA, rCCA
estimators with maximal CV correlation (r2s5-cv); in each case the y-axis
corresponds to the gCCA estimate, while the x-axes correspond to sCCA, gCCA,
rCCA respectively. The sCCA and rCCA estimates have only been registered to
the gCCA estimates up to signs. In each case the full dataset is used to
obtain estimated canonical variates. It is therefore unsurprising that the
middle plot is almost the identity matrix; this only says that the fitted
canonical variates for the CV-optimal graphical lasso estimates are near-
orthogonal, as would be expected for small penalty parameters.

Figure 6 illustrates why it is important to take a subspace perspective when
comparing CCA estimates. Indeed, the first few row and column sums are close
to 1, despite the matrices not being close to the identity. In particular,
observe that sCCA appears to have ‘mixed-up’ the 2nd and 3rd directions
relative to gCCA and rCCA. On closer inspection, observe that there are also
many small, but significant off-diagonal entries which explain why certain
column sums are significantly higher than their corresponding diagonal
entries. It is particularly interesting to see how stable the estimated top-4
subspace is given there was no significant separation between the 4th and 5th
cross validated correlations in 4.

Figure 15 in Section I.1 plots the analogous squared overlap matrices but
where the sCCA and rCCA estimates have been registered to the gCCA estimates
up to orthogonal transformations. In this case the sCCA and rCCA overlap
matrices also become very close to the identity matrix. This again illustrates
the importance of taking a subspace perspective, and the benefits of
orthogonal registration.

Many other versions of overlap matrices can also give useful insights, and
Section I.1 illustrates a few possibilities. For example, one could use a
single algorithm and compare along the trajectory of penalty parameter (Figure
16), or indeed compare between different CV folds which gives a richer view of
the stability of variate subspaces (Figure 14).

####  8.3.5 Biplots

Figure 7: 2D biplots for the microbiome dataset, for sCCA, gCCA, rCCA methods
with respective r2s5-cv-optimal penalty parameters. The sCCA and rCCA
estimates have been registered with the central gCCA estimate up to orthogonal
transformations in variate space. Variates in question are from K0s. C0
variables are plotted with large, coloured, distinctive symbols that are
consistent between the six plots but otherwise arbitrary. K0 variables are
plotted with small dark green symbols to avoid cluttering the plot.

Having developed an understanding of how different estimates compare, we are
now in a good position to inspect certain estimates more closely. One
particularly useful way to do so is to use biplots, introduced in Section 7.3;
additional registration and thresholding can help with comparison and
interpretation of multiple biplots.

The top row of Figure 7 shows such biplots for the first two canonical
variates on the Microbiome dataset, using sample variates from the K0 view;
the second row plots the corresponding weights. The columns correspond to
sCCA, gCCA, rCCA respectively; the penalty used for each algorithm is that
maximising the CV sum of top-5 squared correlations (r2s5-cv). The following
processing steps were used to help comparison between the plots. The sCCA and
rCCA estimates were registered with the gCCA estimate up to orthogonal
transformations (see Section 7.1). The thresholding was different for biplots
and weight plots: for biplots we selected variables with sums of squared
correlation greater than 0.4 for the two gCCA directions (this same subset of
variables is plotted for each estimator); for weights we selected variables
whose first two gCCA weights had sum of absolute values greater than 0.07
(again same variables plotted for each estimator). There are many more K0
variables than C0 variables with significant structure correlations; we plot
the K0 variables with smaller symbols to highlight the C0 variables, which
better illustrate our main observation.

Firstly, Figure 7 illustrates that different methods may estimate similar
canonical variates, but with very different weights. Indeed, loadings plots
are almost identical; and, though there are some similarities in the weight
plots (such as the red-isosceles triangle in the top, corresponding to C00163)
they look much more different.

It may be of interest to further contrast the weight and loading plots. For
example, one could also compare the subset of variables selected for the
weight plot with those selected for the loading plot. Here, the subsets are
generally similar; but there are also some notable exceptions, such as the
variable C00163, which has the largest magnitude weight, but does not attain
required threshold for structure correlation plot.

These biplots give even richer information in three dimensions. We demonstrate
this in the following section, with further examples in Section I.1. It is
particularly convenient to explore such plots with modern plotting libraries
that allow one to rotate the plots interactively. The plots can also be
augmented by label information. This leads to the striking plots in Figure 22.

####  8.3.6 Biological comment

Figure 8: 3D biplot for Microbiome dataset, using CV-optimal gCCA, with
different pathways colored, see Section 8.3.6 for full details

Figure 8 shows a 3D biplot for Microbiome dataset using the r2s5-cv-optimal
penalty parameter. Our Microbiome dataset measures metabolites and different
gene pathways for healthy and diseased patients with IBD. In the 3D biplot we
examined the pathways enriched on either pole of the plot to hypothesize any
latent variables that may explain the opposing pathways. On one pole we
observe pathways and metabolites that are linked to lipid metabolism. Lipid
metabolism is a crucial component of bacterial growth and survival, and
includes pathways such as: Fatty Acid Synthesis, Lipid A Synthesis,
Phospholipid Synthesis and Lipid Desaturation (Ghazalpour et al., 2016). On
the other pole we observed many pathways/metabolites linked to amino acid
metabolism and several genes that corresponded to ABC Transporters. Amino acid
metabolism in microbes generally regards either the building or breakdown of
certain amino acids which are the building blocks of proteins, the functional
unit of cells. ABC (ATP Binding Cassette) Transporters use ATP to transport
various chemical substrates across a cell membrane (Chandel, 2021; Hollenstein
et al., 2007). We have highlighted three biological processes with large
coloured markers: Amino Acid Metabolism: blue, ABC transporters: red, Lipid
Metabolism: green.

During IBD disease progression a multitude of shifts can change in the
microbial flora of the gut (Kostic et al., 2014). A flare up in IBD is noted
to include an influx of host immune cells to aid and defend different tissue
regions within the gut (Neurath, 2019; Zhang et al., 2006). These changes in
IBD can potentially describe the driven differences between amino acid
metabolism/ABC transporters and lipid metabolism. On a more molecular level,
changes in nutrient availability for the microbes due to diet or inflammation
may lead to an availability of more amino acids and their corresponding
transporters and a down regulation of certain lipid sequestration. On the
other hand, if the environment has many lipids due to a breakdown of cell
membranes, many of those metabolic pathways may be upregulated in the
remaining taxonomy. Additional community pressures due to diet changes or
environmental changes may give rise to the growth of certain microbes and that
may include pathways of either amino acid metabolites or lipid intermediates.
With the aid of CCA, we are able to explore groups of genes and metabolites
that are maximally correlated, shedding light on how the genetic catalogue of
the perturbed IBD microbiome might affect metabolic pathways.

###  8.4 Contrast with datasets from previous CCA studies

We compare with two datasets previously used to evaluate CCA methods:

  * •

Nutrimouse: a dataset previously proposed for demonstration of various
multiview learning techniques (Rodosthenous et al., 2020), available in the
CCA R package (González and Déjean, 2021), and originally from Martin et al.
(2007). This contains the expression measure of 120 genes potentially involved
in nutritional problems and the concentrations of 21 hepatic fatty acids
(lipids) for 40 mice.

  * •

BreastData: proposed in Witten et al. (2009) for demonstration of sPLS with
data available in the PMA R package (Tibshirani, 2020) and originally from
Chin et al. (2006). This contains gene expression (RNA, $q=19672$ genes) and
DNA copy number measurements ($p=2149$ CGH spots — comparative genomic
hybridization) on a set of 89 samples.

Both of these datasets have far fewer samples than the Microbiome dataset. The
Nutrimouse dataset has smaller dimension, with ratios $(p/n,q/n)=(0.525,3)$.
By contrast, the Breastdata dataset has much larger dimension, and much larger
ratios $(p/n,q/n)=(24.1,221)$. We summarise the datasets in Table 3 for
reference.

Table 3: Summary of different real datasets analysed Dataset | n | p | q | X-variables | Y-variables
---|---|---|---|---|---
Microbiome | 458 | 711 | 142 | Metabolites (C0s) | Proteins (K0s)
Nutrimouse | 40 | 21 | 120 | Lipids | Gene expression
BreastData | 89 | 2149 | 19672 | DNA copy number | RNA gene expression

We now give a summary of conclusions from applying the techniques of Section
8.3 to these datasets. We leave all justification to Sections I.2 and I.3,
where we state observations and supporting results more systematically.

Despite the smaller sample size, most of our previous observations remain
valid. Indeed, many successive estimated canonical pairs have significant
cross validation signal; variates are significantly more stable than weights;
and the three genuine CCA methods captured similar correlation signal. Again
sPLS captured less signal than the CCA methods, and was relatively unstable
between folds, but this was somewhat less significant than for the Microbiome
set.

However there are some also subtle qualitative differences between the
Nutrimouse and the Microbiome cases. Firstly, Figure 18 shows that very small
penalty parameters are almost optimal in terms of correlation captured for
each of rCCA, sCCA, gCCA. In particular, for gCCA, taking a very small penalty
of $10^{-4}$ seems to give near-optimal correlation captured while also having
near-optimal stability. Secondly, the relevant trajectory comparison matrices
show that the ridge estimator is almost constant for small parameter values,
and only starts to change for parameters approaching 1. This behaviour also
manifests in all of the other diagnostic plots — correlations along
trajectories are very flat, correlation decay has persistent colour structure,
overlap matrices are near identity, and biplots vary only a little with
penalty parameter. Thirdly, gCCA also moves only a little along trajectories -
and is always very close to the ridge solution; again this manifests in each
different plot.

There are more pronounced qualitative differences with the Breastdata dataset.
The main difference is that the estimated variates are no longer stable
between algorithms; indeed, sCCA and rCCA now give very different estimates,
even after registration. This can be seen at a high level with the relevant
trajectory comparison matrices and at a more detailed level with the relevant
overlap matrices. In addition, though sCCA is now very sensitive to penalty
parameter (overlap matrices show that only the top direction is stable), it is
less sensitive to fold at the optimal penalty (the top-3 subspace is stable).
By contrast, rCCA is now almost invariant to its penalty parameter, while
still being fairly stable with respect to fold; this is remarkable and
suggests some mathematical phenomenon of potential interest.

To summarise these results, it is helpful to think about signal to noise
ratios (SNR). In high SNR regimes, such as with Microbiome and Nutrimouse
dataset, we may expect to consistently estimate stable canonical variates,
even if the directions are unstable. But for lower SNR regimes we may expect
to find very different canonical variate estimates — but that these may still
contain significant CV correlations signal — so give different insights about
the data. In these regimes it may be particularly useful to combine insights
from a variety of different regularised CCA methods, which will pick up
different sorts of signal in the data.

##  9 Discussion

Many regularised CCA methods are now practical for high-dimensional data due
to recent computational and theoretical advances. However, it is not clear
that the assumptions behind the regularisation are appropriate; moreover, it
is not accepted how to choose between different methods or indeed interpret
the output. Our first contribution was to propose a very different regularised
CCA method, gCCA, based on the Graphical Lasso, and to show it has reassuring
theoretically properties. Our second contribution was to investigate different
criteria for evaluating CCA methods, for comparing between different CCA
estimates, and indeed for interpreting the estimates. We observed that it is
much easier to estimate canonical variates than canonical directions, and to
consider successive subspaces rather than successive canonical pairs — which
was consistent with our theoretical results and intuition. This led to a
concrete framework for applying and interpreting different CCA methods to real
data; this framework culminates in the powerful visualisation tool of biplots.
Within this framework, our gCCA estimator often outperformed existing methods;
moreover, by using it alongside existing methods, we were able to get a much
broader perspective on high-dimensional CCA. We recommend future practitioners
use a variety of methods, and hope this leads to a wider use of CCA within the
statistics community.

Our work has opened up many directions for future research. To start with,
there are natural extensions of our gCCA method that might be worth exploring;
the plug-in graphical lasso approach can also be applied with any multiview
CCA framework (Tenenhaus and Tenenhaus, 2011); there are also extensions of
the graphical lasso to deal with (unobserved) latent variables (Chandrasekaran
et al., 2012), which have looser theoretical assumptions, are correspondingly
more flexible, and may perform better empirically. Additionally, there has
been a variety of recent work on improving the computational cost of the
graphical lasso, with some methods already implemented in gglasso (Schaipp et
al., 2021), and which may help scale up the method to higher-dimensional
situations.

There are also open theoretical questions surrounding CCA more generally. To
start with, there is little theoretical understanding of ridge regularised CCA
(rCCA); this is a simple, effective and popular high-dimensional CCA method
and therefore is important to understand. Firstly, we observed that even in
very high-dimensional regimes the ridge methods were reasonably stable and
recovered consistent signal. Secondly, in the very high-dimensional regimes,
ridge CCA was almost unchanged by varying the penalty parameter. Given recent
advances in random matrix theory applied to ridge estimators for linear
regression, we suspect both of these questions may be tractable. However, we
also observed that the sparse CCA methods also frequently recovered top
variate subspaces for CCA, even when the true directions were not sparse. Our
gCCA method was also remarkably effective at recovering variate subspaces even
in the mispecified case. More generally there appears to be room for theory
bridging the gap from Gao et al. (2016); Ma and Li (2020) to our observations
of robustness to mis-specification.

However, the most impactful area for research would be on the practical side.
To start it would be helpful to integrate the evaluation criteria and
comparison tools we have proposed into existing open source implementations of
CCA methods, such as cca-zoo (Chapman and Wang, 2021). It would also be
helpful to have a public repository of multiview datasets, that practitioners
could use to develop intuition for CCA, or indeed related multiview learning
algorithms. There is also scope for a broader systematic numerical comparison
between CCA algorithms than we have performed here, using a wider array of
synthetic and real datasets, and investigating how the various CCA criteria
interact with downstream tasks. Of particular interest would be further
understanding how CCA integrates with non-linear dimension reduction
techniques in a broader EDA pipeline. There is also scope to develop more
flexible regularised CCA techniques that can incorporate more flexible
biological priors, and to investigate more creative multiple applications of
CCA, perhaps by comparing CCA on different subgroups of data. We hope the
framework and analysis we have provided will encourage further research into
these interesting directions.

##  10 Competing interests

No competing interest is declared.

##  11 Author contributions statement

LW is largely responsible for the theoretical results about the gCCA procedure
and the mathematical presentation in the supplementary materials. He conceived
of the methodology described in Section 7. He did the bulk of the writing,
wrote all the code for the project, designed and carried out all numerical
experiments. SB conceived the core idea of applying the graphical lasso to
CCA, supervised the project closely, and edited the article. KT obtained the
Microbiome dataset analysed in Section 8 and wrote the corresponding
commentary in Section 8.3.6.

##  12 Acknowledgments

The authors thank James Chapman for valuable discussions and comments on an
early draft of this paper. LW is supported by the UK Engineering and Physical
Sciences Research Council (EPSRC) under grant number EP/V52024X/1. KT is
supported by a Gates Cambridge Scholarship.

## References

  * Anderson [2003] T. W. Anderson.  _An introduction to multivariate statistical analysis_.  Wiley series in probability and statistics. Wiley-Interscience, Hoboken, N.J, 3rd ed edition, 2003.  ISBN 978-0-471-36091-9.
  * Andrew et al. [2013] G. Andrew, R. Arora, J. Bilmes, and K. Livescu.  Deep canonical correlation analysis.  In _International conference on machine learning_ , pages 1247–1255. PMLR, 2013.
  * Bach and Jordan [2002] F. R. Bach and M. I. Jordan.  Kernel Independent Component Analysis.  _Journal of Machine Learning Research_ , 3(Jul):1–48, 2002.  ISSN ISSN 1533-7928.  URL https://www.jmlr.org/papers/v3/bach02a.html.
  * Bach and Jordan [2005] F. R. Bach and M. I. Jordan.  A Probabilistic Interpretation of Canonical Correlation Analysis.  Technical report, University of California, Berkeley, 2005.
  * Bhatia [1997] R. Bhatia.  _Matrix Analysis_ , volume 169 of _Graduate Texts in Mathematics_.  Springer, New York, NY, 1997.  ISBN 978-1-4612-6857-4 978-1-4612-0653-8.  doi: 10.1007/978-1-4612-0653-8.  URL http://link.springer.com/10.1007/978-1-4612-0653-8.
  * Bishop [2006] C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, Aug. 2006.  ISBN 978-0-387-31073-2.  Google-Books-ID: qWPwnQEACAAJ.
  * Borga [1998] M. Borga.  _Learning Multidimensional Signal Processing_.  PhD thesis, Linköping University, 1998.  URL http://urn.kb.se/resolve?urn=urn:nbn:se:liu:diva-54341.  Publisher: Linköping University Electronic Press.
  * Boyd [2010] S. Boyd.  Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers.  _Foundations and Trends® in Machine Learning_ , 3(1):1–122, 2010.  ISSN 1935-8237, 1935-8245.  doi: 10.1561/2200000016.  URL http://www.nowpublishers.com/article/Details/MAL-016.
  * Boyd and Vandenberghe [2004] S. P. Boyd and L. Vandenberghe.  _Convex optimization_.  Cambridge University Press, Cambridge, UK ; New York, 2004.  ISBN 978-0-521-83378-3.
  * Cai et al. [2011] T. Cai, W. Liu, and X. Luo.  A constrained $\ell$-1 minimization approach to sparse precision matrix estimation.  _Journal of the American Statistical Association_ , 106(494):594–607, 2011.
  * Cai and Zhang [2020] T. T. Cai and A. Zhang.  Rate-Optimal Perturbation Bounds for Singular Subspaces with Applications to High-Dimensional Statistics, June 2020.  URL http://arxiv.org/abs/1605.00353.  arXiv:1605.00353 [math, stat].
  * Carlsson [2021] M. Carlsson.  von Neumann’s trace inequality for Hilbert–Schmidt operators.  _Expositiones Mathematicae_ , 39(1):149–157, Mar. 2021.  ISSN 0723-0869.  doi: 10.1016/j.exmath.2020.05.001.  URL https://www.sciencedirect.com/science/article/pii/S0723086920300220.
  * Chandel [2021] N. S. Chandel.  Amino Acid Metabolism.  _Cold Spring Harbor Perspectives in Biology_ , 13(4):a040584, Jan. 2021.  ISSN , 1943-0264.  doi: 10.1101/cshperspect.a040584.  URL http://cshperspectives.cshlp.org/content/13/4/a040584.  Company: Cold Spring Harbor Laboratory Press Distributor: Cold Spring Harbor Laboratory Press Institution: Cold Spring Harbor Laboratory Press Label: Cold Spring Harbor Laboratory Press Publisher: Cold Spring Harbor Lab.
  * Chandrasekaran et al. [2012] V. Chandrasekaran, P. A. Parrilo, and A. S. Willsky.  Latent variable graphical model selection via convex optimization.  _The Annals of Statistics_ , 40(4):1935–1967, Aug. 2012.  ISSN 0090-5364, 2168-8966.  doi: 10.1214/11-AOS949.  URL https://projecteuclid.org/journals/annals-of-statistics/volume-40/issue-4/Latent-variable-graphical-model-selection-via-convex-optimization/10.1214/11-AOS949.full.  Publisher: Institute of Mathematical Statistics.
  * Chapman and Wang [2021] J. Chapman and H.-T. Wang.  CCA-Zoo: A collection of Regularized, Deep Learning based, Kernel, and Probabilistic CCA methods in a scikit-learn style framework.  _Journal of Open Source Software_ , 6(68):3823, Dec. 2021.  ISSN 2475-9066.  doi: 10.21105/joss.03823.  URL https://joss.theoj.org/papers/10.21105/joss.03823.
  * Chapman et al. [2023] J. Chapman, L. Wells, and A. L. Aguila.  Efficient Algorithms for the CCA Family: Unconstrained Objectives with Unbiased Gradients, Nov. 2023.  URL http://arxiv.org/abs/2310.01012.  arXiv:2310.01012 [cs, stat].
  * Chen et al. [2013] M. Chen, C. Gao, Z. Ren, and H. H. Zhou.  Sparse CCA via Precision Adjusted Iterative Thresholding.  Technical Report arXiv:1311.6186, arXiv, Nov. 2013.  URL http://arxiv.org/abs/1311.6186.  arXiv:1311.6186 [math, stat] type: article.
  * Chin et al. [2006] K. Chin, S. DeVries, J. Fridlyand, P. T. Spellman, R. Roydasgupta, W.-L. Kuo, A. Lapuk, R. M. Neve, Z. Qian, T. Ryder, F. Chen, H. Feiler, T. Tokuyasu, C. Kingsley, S. Dairkee, Z. Meng, K. Chew, D. Pinkel, A. Jain, B. M. Ljung, L. Esserman, D. G. Albertson, F. M. Waldman, and J. W. Gray.  Genomic and transcriptional aberrations linked to breast cancer pathophysiologies.  _Cancer Cell_ , 10(6):529–541, Dec. 2006.  ISSN 1535-6108.  doi: 10.1016/j.ccr.2006.10.009.  URL https://www.sciencedirect.com/science/article/pii/S1535610806003151.
  * Cover and Thomas [2012] T. M. Cover and J. A. Thomas.  _Elements of Information Theory_.  John Wiley & Sons, Nov. 2012.  ISBN 978-1-118-58577-1.  Google-Books-ID: VWq5GG6ycxMC.
  * Friedman et al. [2008] J. Friedman, T. Hastie, and R. Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, July 2008.  ISSN 1465-4644.  doi: 10.1093/biostatistics/kxm045.  URL https://doi.org/10.1093/biostatistics/kxm045.
  * Fu et al. [2017] X. Fu, K. Huang, M. Hong, N. D. Sidiropoulos, and A. M.-C. So.  Scalable and Flexible Multiview MAX-VAR Canonical Correlation Analysis.  _IEEE Transactions on Signal Processing_ , 65(16):4150–4165, Aug. 2017.  ISSN 1053-587X, 1941-0476.  doi: 10.1109/TSP.2017.2698365.  URL http://arxiv.org/abs/1605.09459.  arXiv:1605.09459 [stat].
  * Gabriel [1971] K. R. Gabriel.  The Biplot Graphic Display of Matrices with Application to Principal Component Analysis.  _Biometrika_ , 58(3):453–467, 1971.  ISSN 0006-3444.  doi: 10.2307/2334381.  URL https://www.jstor.org/stable/2334381.  Publisher: [Oxford University Press, Biometrika Trust].
  * Gao et al. [2015] C. Gao, Z. Ma, Z. Ren, and H. H. Zhou.  Minimax estimation in sparse canonical correlation analysis.  _The Annals of Statistics_ , 43(5):2168–2197, Oct. 2015.  ISSN 0090-5364, 2168-8966.  doi: 10.1214/15-AOS1332.  URL https://projecteuclid.org/journals/annals-of-statistics/volume-43/issue-5/Minimax-estimation-in-sparse-canonical-correlation-analysis/10.1214/15-AOS1332.full.  Publisher: Institute of Mathematical Statistics.
  * Gao et al. [2016] C. Gao, Z. Ma, and H. H. Zhou.  Sparse CCA: Adaptive Estimation and Computational Barriers.  _arXiv:1409.8565 [math, stat]_ , Apr. 2016.  URL http://arxiv.org/abs/1409.8565.  arXiv: 1409.8565.
  * Ghazalpour et al. [2016] A. Ghazalpour, I. Cespedes, B. J. Bennett, and H. Allayee.  Expanding Role of Gut Microbiota in Lipid Metabolism.  _Current opinion in lipidology_ , 27(2):141–147, Apr. 2016.  ISSN 0957-9672.  doi: 10.1097/MOL.0000000000000278.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5125441/.
  * González and Déjean [2021] I. González and S. Déjean.  CCA: Canonical Correlation Analysis, Mar. 2021.  URL https://CRAN.R-project.org/package=CCA.
  * Gorsuch [2014] R. L. Gorsuch.  _Factor analysis: Classic edition_.  Routledge, 2014.
  * Gu and Wu [2018] F. Gu and H. Wu.  Simultaneous canonical correlation analysis with invariant canonical loadings.  _Behaviormetrika_ , 45(1):111–132, Apr. 2018\.  ISSN 1349-6964.  doi: 10.1007/s41237-017-0042-8.  URL https://doi.org/10.1007/s41237-017-0042-8.
  * Hardoon and Shawe-Taylor [2011] D. R. Hardoon and J. Shawe-Taylor.  Sparse canonical correlation analysis.  _Machine Learning_ , 83(3):331–353, June 2011\.  ISSN 0885-6125, 1573-0565.  doi: 10.1007/s10994-010-5222-7.  URL http://link.springer.com/10.1007/s10994-010-5222-7.
  * Hastie et al. [2015] T. Hastie, R. Tibshirani, and M. J. Wainwright.  _Statistical Learning with Sparsity: the Lasso and Generalizations_.  CRC Monographs on Statistics and Applied Probability. Chapman & Hall, 2015.  URL https://hastie.su.domains/StatLearnSparsity/.
  * Hollenstein et al. [2007] K. Hollenstein, R. J. Dawson, and K. P. Locher.  Structure and mechanism of ABC transporter proteins.  _Current Opinion in Structural Biology_ , 17(4):412–418, Aug. 2007.  ISSN 0959-440X.  doi: 10.1016/j.sbi.2007.07.003.  URL https://www.sciencedirect.com/science/article/pii/S0959440X07001029.
  * Hotelling [1936] H. Hotelling.  Relations Between Two Sets of Variates.  _Biometrika_ , 28(3/4):321–377, 1936.  ISSN 0006-3444.  doi: 10.2307/2333955.  URL https://www.jstor.org/stable/2333955.  Publisher: [Oxford University Press, Biometrika Trust].
  * Jolliffe [1986] I. T. Jolliffe.  _Principal Component Analysis_.  Springer Series in Statistics. Springer, New York, NY, 1986.  ISBN 978-1-4757-1906-2 978-1-4757-1904-8.  doi: 10.1007/978-1-4757-1904-8.  URL http://link.springer.com/10.1007/978-1-4757-1904-8.
  * Jolliffe and Cadima [2016] I. T. Jolliffe and J. Cadima.  Principal component analysis: a review and recent developments.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 374(2065):20150202, Apr. 2016.  doi: 10.1098/rsta.2015.0202.  URL https://royalsocietypublishing.org/doi/10.1098/rsta.2015.0202.  Publisher: Royal Society.
  * Klami et al. [2013] A. Klami, S. Virtanen, and S. Kaski.  Bayesian Canonical Correlation Analysis.  _Journal of Machine Learning Research_ , 14(30):965–1003, 2013.  ISSN 1533-7928.  URL http://jmlr.org/papers/v14/klami13a.html.
  * Kostic et al. [2014] A. D. Kostic, R. J. Xavier, and D. Gevers.  The microbiome in inflammatory bowel disease: current status and the future ahead.  _Gastroenterology_ , 146(6):1489–1499, May 2014\.  ISSN 1528-0012.  doi: 10.1053/j.gastro.2014.02.009.
  * Laha and Mukherjee [2021] N. Laha and R. Mukherjee.  On Support Recovery with Sparse CCA: Information Theoretic and Computational Limits.  Technical Report arXiv:2108.06463, arXiv, Aug. 2021.  URL http://arxiv.org/abs/2108.06463.  arXiv:2108.06463 [math, stat] type: article.
  * Laha et al. [2022] N. Laha, N. Huey, B. Coull, and R. Mukherjee.  On Statistical Inference with High Dimensional Sparse CCA.  Technical Report arXiv:2109.11997, arXiv, Feb. 2022.  URL http://arxiv.org/abs/2109.11997.  arXiv:2109.11997 [math, stat] type: article.
  * Lauritzen [1996] S. L. Lauritzen.  _Graphical Models_.  Oxford Statistical Science Series. Oxford University Press, Oxford, New York, May 1996.  ISBN 978-0-19-852219-5.
  * Lloyd-Price et al. [2019] J. Lloyd-Price, C. Arze, A. N. Ananthakrishnan, M. Schirmer, J. Avila-Pacheco, T. W. Poon, E. Andrews, N. J. Ajami, K. S. Bonham, C. J. Brislawn, D. Casero, H. Courtney, A. Gonzalez, T. G. Graeber, A. B. Hall, K. Lake, C. J. Landers, H. Mallick, D. R. Plichta, M. Prasad, G. Rahnavard, J. Sauk, D. Shungin, Y. Vázquez-Baeza, R. A. White, J. Braun, L. A. Denson, J. K. Jansson, R. Knight, S. Kugathasan, D. P. B. McGovern, J. F. Petrosino, T. S. Stappenbeck, H. S. Winter, C. B. Clish, E. A. Franzosa, H. Vlamakis, R. J. Xavier, and C. Huttenhower.  Multi-omics of the gut microbial ecosystem in inflammatory bowel diseases.  _Nature_ , 569(7758):655–662, May 2019.  ISSN 1476-4687.  doi: 10.1038/s41586-019-1237-9.  URL https://www.nature.com/articles/s41586-019-1237-9.  Number: 7758 Publisher: Nature Publishing Group.
  * Loh and Wainwright [2013] P.-L. Loh and M. J. Wainwright.  Structure estimation for discrete graphical models: Generalized covariance matrices and their inverses.  _The Annals of Statistics_ , 41(6):3022–3049, 2013.  ISSN 0090-5364.  URL https://www.jstor.org/stable/23566757.  Publisher: Institute of Mathematical Statistics.
  * Ma and Li [2020] Z. Ma and X. Li.  Subspace perspective on canonical correlation analysis: Dimension reduction and minimax rates.  _Bernoulli_ , 26(1), Feb. 2020.  ISSN 1350-7265.  doi: 10.3150/19-BEJ1131.  URL https://projecteuclid.org/journals/bernoulli/volume-26/issue-1/Subspace-perspective-on-canonical-correlation-analysis--Dimension-reduction-and/10.3150/19-BEJ1131.full.
  * Mai and Zhang [2019] Q. Mai and X. Zhang.  An iterative penalized least squares approach to sparse canonical correlation analysis.  _Biometrics_ , 75(3):734–744, 2019.  ISSN 1541-0420.  doi: 10.1111/biom.13043.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.13043.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/biom.13043.
  * Martin et al. [2007] P. G. P. Martin, H. Guillou, F. Lasserre, S. Déjean, A. Lan, J.-M. Pascussi, M. Sancristobal, P. Legrand, P. Besse, and T. Pineau.  Novel aspects of PPARalpha-mediated regulation of lipid and xenobiotic metabolism revealed through a nutrigenomic study.  _Hepatology (Baltimore, Md.)_ , 45(3):767–777, Mar. 2007.  ISSN 0270-9139.  doi: 10.1002/hep.21510.
  * Mihalik et al. [2022] A. Mihalik, J. Chapman, R. A. Adams, N. R. Winter, F. S. Ferreira, J. Shawe-Taylor, and J. Mourão-Miranda.  Canonical Correlation Analysis and Partial Least Squares for identifying brain-behaviour associations: a tutorial and a comparative study.  _Biological Psychiatry: Cognitive Neuroscience and Neuroimaging_ , Aug. 2022.  ISSN 2451-9022.  doi: 10.1016/j.bpsc.2022.07.012.  URL https://www.sciencedirect.com/science/article/pii/S2451902222001859.
  * Neurath [2019] M. F. Neurath.  Targeting immune cell circuits and trafficking in inflammatory bowel disease.  _Nature Immunology_ , 20(8):970–979, Aug. 2019\.  ISSN 1529-2916.  doi: 10.1038/s41590-019-0415-0.
  * Parikh and Boyd [2014] N. Parikh and S. Boyd.  Proximal Algorithms.  _Foundations and Trends in Optimization_ , 1(3):127–239, 2014.  ISSN 2167-3888.  doi: 10.1561/2400000003.  URL http://dx.doi.org/10.1561/2400000003.
  * [48] A. Prokhorov.  Partial correlation coefficient - Encyclopedia of Mathematics.  URL https://encyclopediaofmath.org/wiki/Partial_correlation_coefficient#:~:text=A%20partial%20correlation%20coefficient%20is,the%20remaining%20variables%20is%20eliminated.
  * Ravikumar et al. [2011] P. Ravikumar, M. J. Wainwright, G. Raskutti, and B. Yu.  High-dimensional covariance estimation by minimizing l1-penalized log-determinant divergence.  _Electronic Journal of Statistics_ , 5(none):935–980, Jan. 2011.  ISSN 1935-7524, 1935-7524.  doi: 10.1214/11-EJS631.  URL https://projecteuclid.org/journals/electronic-journal-of-statistics/volume-5/issue-none/High-dimensional-covariance-estimation-by-minimizing-%e2%84%931-penalized-log-determinant/10.1214/11-EJS631.full.  Publisher: Institute of Mathematical Statistics and Bernoulli Society.
  * Rodosthenous et al. [2020] T. Rodosthenous, V. Shahrezaei, and M. Evangelou.  Integrating multi-OMICS data through sparse canonical correlation analysis for the prediction of complex traits: a comparison study.  _Bioinformatics_ , 36(17):4616–4625, Nov. 2020\.  ISSN 1367-4803.  doi: 10.1093/bioinformatics/btaa530.  URL https://doi.org/10.1093/bioinformatics/btaa530.
  * Rosipal and Krämer [2006] R. Rosipal and N. Krämer.  Overview and Recent Advances in Partial Least Squares.  In C. Saunders, M. Grobelnik, S. Gunn, and J. Shawe-Taylor, editors, _Subspace, Latent Structure and Feature Selection_ , Lecture Notes in Computer Science, pages 34–51, Berlin, Heidelberg, 2006. Springer.  ISBN 978-3-540-34138-3.  doi: 10.1007/11752790˙2.
  * Roweis and Ghahramani [1999] S. Roweis and Z. Ghahramani.  A Unifying Review of Linear Gaussian Models.  _Neural Computation_ , 11(2):305–345, Feb. 1999\.  ISSN 0899-7667, 1530-888X.  doi: 10.1162/089976699300016674.  URL https://direct.mit.edu/neco/article/11/2/305-345/6249.
  * Safo et al. [2018a] S. E. Safo, J. Ahn, Y. Jeon, and S. Jung.  Sparse generalized eigenvalue problem with application to canonical correlation analysis for integrative analysis of methylation and gene expression data.  _Biometrics_ , 74(4):1362–1371, 2018a.  ISSN 1541-0420.  doi: 10.1111/biom.12886.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12886.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/biom.12886.
  * Safo et al. [2018b] S. E. Safo, S. Li, and Q. Long.  Integrative analysis of transcriptomic and metabolomic data via sparse canonical correlation analysis with incorporation of biological information.  _Biometrics_ , 74(1):300–312, 2018b.  ISSN 1541-0420.  doi: 10.1111/biom.12715.  URL https://onlinelibrary.wiley.com/doi/abs/10.1111/biom.12715.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1111/biom.12715.
  * Schaipp et al. [2021] F. Schaipp, C. L. Müller, and O. Vlasovets.  GGLasso – a Python package for General Graphical Lasso computation, Oct. 2021.  URL http://arxiv.org/abs/2110.10521.  Number: arXiv:2110.10521 arXiv:2110.10521 [stat].
  * Solari et al. [2019] O. S. Solari, J. B. Brown, and P. J. Bickel.  Sparse Canonical Correlation Analysis via Concave Minimization.  Technical Report arXiv:1909.07947, arXiv, Sept. 2019.  URL http://arxiv.org/abs/1909.07947.  arXiv:1909.07947 [cs, stat] type: article.
  * Stewart and Sun [1990] G. W. Stewart and J.-G. Sun.  _Matrix Perturbation Theory_.  ACADEMIC PressINC, July 1990.  ISBN 978-1-4933-0199-7.  Google-Books-ID: bIYEogEACAAJ.
  * Suo et al. [2017] X. Suo, V. Minden, B. Nelson, R. Tibshirani, and M. Saunders.  Sparse canonical correlation analysis.  _arXiv:1705.10865 [stat]_ , June 2017.  URL http://arxiv.org/abs/1705.10865.  arXiv: 1705.10865.
  * Tan et al. [2018] K. M. Tan, Z. Wang, H. Liu, and T. Zhang.  Sparse Generalized Eigenvalue Problem: Optimal Statistical Rates via Truncated Rayleigh Flow, Aug. 2018.  URL http://arxiv.org/abs/1604.08697.  arXiv:1604.08697 [stat].
  * Tenenhaus and Tenenhaus [2011] A. Tenenhaus and M. Tenenhaus.  Regularized Generalized Canonical Correlation Analysis.  _Psychometrika_ , 76(2):257–284, Apr. 2011.  ISSN 0033-3123, 1860-0980.  doi: 10.1007/s11336-011-9206-8.  URL http://link.springer.com/10.1007/s11336-011-9206-8.
  * Tenenhaus et al. [2014] A. Tenenhaus, C. Philippe, V. Guillemot, K.-A. Le Cao, J. Grill, and V. Frouin.  Variable selection for generalized canonical correlation analysis.  _Biostatistics_ , 15(3):569–583, July 2014.  ISSN 1465-4644.  doi: 10.1093/biostatistics/kxu001.  URL https://doi.org/10.1093/biostatistics/kxu001.
  * ter Braak [1990] C. J. F. ter Braak.  Interpreting canonical correlation analysis through biplots of structure correlations and weights.  _Psychometrika_ , 55(3):519–531, Sept. 1990.  ISSN 0033-3123, 1860-0980.  doi: 10.1007/BF02294765.  URL http://link.springer.com/10.1007/BF02294765.
  * Tibshirani [2020] D. W. a. R. Tibshirani.  PMA: Penalized Multivariate Analysis, Feb. 2020.  URL https://CRAN.R-project.org/package=PMA.
  * Uurtio et al. [2017] V. Uurtio, J. M. Monteiro, J. Kandola, J. Shawe-Taylor, D. Fernandez-Reyes, and J. Rousu.  A Tutorial on Canonical Correlation Methods.  _ACM Computing Surveys_ , 50(6):95:1–95:33, Nov. 2017.  ISSN 0360-0300.  doi: 10.1145/3136624.  URL https://doi.org/10.1145/3136624.
  * Vinod [1976] H. D. Vinod.  Canonical ridge and econometrics of joint production.  _Journal of Econometrics_ , 4(2):147–166, May 1976.  ISSN 0304-4076.  doi: 10.1016/0304-4076(76)90010-5.  URL https://www.sciencedirect.com/science/article/pii/0304407676900105.
  * Wainwright [2019] M. J. Wainwright.  _High-Dimensional Statistics: A Non-Asymptotic Viewpoint_.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, Cambridge, 2019.  ISBN 978-1-108-49802-9.  doi: 10.1017/9781108627771.  URL https://www.cambridge.org/core/books/highdimensional-statistics/8A91ECEEC38F46DAB53E9FF8757C7A4E.
  * Witten et al. [2009] D. M. Witten, R. Tibshirani, and T. Hastie.  A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.  _Biostatistics_ , 10(3):515–534, July 2009.  ISSN 1465-4644, 1468-4357.  doi: 10.1093/biostatistics/kxp008.  URL https://academic.oup.com/biostatistics/article-lookup/doi/10.1093/biostatistics/kxp008.
  * Xiu et al. [2020] X. Xiu, Y. Yang, L. Kong, and W. Liu.  tSSNALM: A fast two-stage semi-smooth Newton augmented Lagrangian method for sparse CCA.  _Applied Mathematics and Computation_ , 383:125272, Oct. 2020.  ISSN 0096-3003.  doi: 10.1016/j.amc.2020.125272.  URL https://www.sciencedirect.com/science/article/pii/S0096300320302411.
  * Yu et al. [2015] Y. Yu, T. Wang, and R. J. Samworth.  A useful variant of the Davis—Kahan theorem for statisticians.  _Biometrika_ , 102(2):315–323, 2015.  ISSN 0006-3444.  URL https://www.jstor.org/stable/43908537.  Publisher: [Oxford University Press, Biometrika Trust].
  * Yuan and Lin [2007] M. Yuan and Y. Lin.  Model Selection and Estimation in the Gaussian Graphical Model.  _Biometrika_ , 94(1):19–35, 2007.  ISSN 0006-3444.  URL https://www.jstor.org/stable/20441351.  Publisher: [Oxford University Press, Biometrika Trust].
  * Zbontar et al. [2021] J. Zbontar, L. Jing, I. Misra, Y. LeCun, and S. Deny.  Barlow Twins: Self-Supervised Learning via Redundancy Reduction, June 2021.  URL http://arxiv.org/abs/2103.03230.  arXiv:2103.03230 [cs, q-bio].
  * Zhang et al. [2006] S. Z. Zhang, X. H. Zhao, and D. C. Zhang.  Cellular and molecular immunopathogenesis of ulcerative colitis.  _Cellular & Molecular Immunology_, 3(1):35–40, Feb. 2006.  ISSN 1672-7681.

##  Appendix A Matrix Analysis background

###  A.1 Warm-up: uniqueness of eigendecomposition

######  Lemma 7.

Let $V,V^{\prime}$ span the same subspace of $\mathbb{R}^{d}$, each with
orthonormal columns. Then $V^{\prime}=VO$ for some orthogonal matrix $O$.

###### Proof.

Note that $VV^{\top}$ is an orthogonal projection onto the shared subspace, so
defining $O\vcentcolon=V^{\top}V^{\prime}$ we do indeed have

$\displaystyle V^{\prime}=VV^{\top}V^{\prime}=VO\;.$

Moreover, again because $VV^{\top}$ is an orthogonal projection

$\displaystyle
O^{\top}O=V^{\prime\top}VV^{\top}V^{\prime}=V^{\prime\top}V^{\prime}=I\;,$

so $O$ is indeed orthogonal. ∎

######  Lemma 8.

Let $A$ be a symmetric matrix. Let $v,v^{\prime}$ be eigenvectors of $A$ with
eigenvalues $\lambda,\lambda^{\prime}$ respectively. If
$\lambda^{\prime}\neq\lambda$ then $\langle v,v^{\prime}\rangle=0$.

###### Proof.

Simply compute

$\displaystyle
0=v^{\top}Av^{\prime}-v^{\prime\top}Av=\lambda^{\prime}v^{\top}v^{\prime}-\lambda
v^{\prime\top}v=(\lambda^{\prime}-\lambda)\langle v,v^{\prime}\rangle$

so dividing through by $\lambda^{\prime}-\lambda$ gives the result. ∎

######  Proposition 9.

Let $A\in\mathbb{R}^{d\times d}$ have eigenvalues
$\lambda_{1}\geq\dots\geq\lambda_{d}$. Suppose we have two matrices
$V,V^{\prime}$ such that

$\displaystyle A=V\Lambda V^{\top}=V^{\prime}\Lambda V^{\prime\top}\;.$

Suppose that an eigenvalue $\lambda$ occurs with multiplicity $\delta$, and
let $V_{\lambda},V_{\lambda}^{\prime}$ denote the corresponding columns of
$V,V^{\prime}$ (i.e. for some $k$ we have
$\lambda_{k-1}>\lambda=\lambda_{k}=\dots=\lambda_{k+\delta-1}>\lambda_{k+\delta}$
and $V_{\lambda}=\begin{pmatrix}V_{k}&\dots&V_{k+\delta-1}\end{pmatrix}$).
Then there is some orthogonal matrix
$O_{\lambda}\in\mathbb{R}^{\delta\times\delta}$ such that

$\displaystyle V_{\lambda}^{\prime}=V_{\lambda}O_{\lambda}\;.$ (27)

###### Proof.

Consider $v^{\prime}_{j}$ for $j\in\\{k,\dots,k+\delta-1\\}$; then
$v^{\prime}_{j}$ is a $\lambda$-eigenvector of $A$. Now consider any
$l\notin\\{k,\dots,k+\delta-1\\}$; $v_{l}$ is a $\lambda_{l}$-eigenvector of
$A$. Therefore by Lemma 8 we have $\langle v_{l},v^{\prime}_{j}\rangle=0$.

But since $V$ spans $\mathbb{R}^{d}$ we must therefore have
$v^{\prime}_{j}\in\operatorname{span}{V_{\lambda}}$. Hence
$\operatorname{span}{V^{\prime}_{\lambda}}\subset\operatorname{span}{V_{\lambda}}$.
By symmetry we deduce that $V_{\lambda},V^{\prime}_{\lambda}$ span the same
$\delta$-dimensional subspace. We conclude by Lemma 7. ∎

###  A.2 SVD

The SVD is a widely used tool in both pure and applied mathematics, which has
many different motivations and interpretations. To make our arguments
concrete, we summarise some relevant properties of the SVD below; these are
all well known, but rarely presented precisely as we desire them, so we give
proofs. Our presentation closely follows that of Stewart and Sun [1990].
Indeed, the statement and proof of the first result below (existence of SVD)
are directly copied from Stewart and Sun [1990], but repeated here to help
with the later, less standard, results.

####  A.2.1 Definition, existence and uniqueness

######  Theorem 10 (Characterisation via $A^{\top}A$).

Let $A\in\mathbb{R}^{p\times q}$ have rank $K$. Let $U\in\mathbb{R}^{p\times
p}$ and $V\in\mathbb{R}^{q\times q}$ be orthogonal matrices partitioned as
$U=\begin{pmatrix}U_{:K}&U_{K:}\end{pmatrix},\>V=\begin{pmatrix}V_{:K}&V_{K:}\end{pmatrix}$
with $U_{:K}\in\mathbb{R}^{p\times K},V_{:K}\in\mathbb{R}^{q\times K}$.Let
$\sigma_{1}\geq\dots\geq\sigma_{K}>0$ and define
$\Sigma_{+}=\operatorname{diag}(\sigma_{1},\dots,\sigma_{K})\in\mathbb{R}^{K\times
K}$, $\Sigma_{p\times
q}=\operatorname{diag}(\sigma_{1},\dots,\sigma_{K},0,\dots,0)\in\mathbb{R}^{p\times
q}$ and $\Sigma_{q\times
q}=\operatorname{diag}(\sigma_{1},\dots,\sigma_{K},0,\dots,0)\in\mathbb{R}^{q\times
q}$.

Then the following are equivalent.

  1. (i)

$V$ defines an eigen-basis for $A^{\top}A$ and $V_{:K}$ determines $U_{:K}$
by:

$U_{:K}=AV_{:K}\Sigma_{+}^{-1}\quad\text{ and }\quad
V^{\top}A^{\top}AV=\begin{pmatrix}\Sigma_{+}^{2}&0\\\
0&0\end{pmatrix}=\Sigma_{q\times q}^{2}$ (28)

  2. (ii)

We have the singular value decomposition

$U^{\top}AV=\begin{pmatrix}\Sigma_{+}&0\\\
0&0\end{pmatrix}=\vcentcolon\Sigma_{p\times q}$ (29)

###### Proof.

(i) $\implies$ (ii):

By considering sub-matrices of Equation 28 we have

$\displaystyle V_{:K}^{\top}A^{\top}AV_{:K}=\Sigma_{+}^{2},\quad
V_{K:}^{\top}A^{\top}AV_{K:}=0\;.$ (30)

The second of these equations implies that in fact $AV_{K:}=0$.

Since $U$ is orthogonal we have $U_{:K}^{\top}U_{:K}=I$. Substituting this
into the first equation of Equation 30 gives
$I=U_{:K}^{\top}U_{:K}=U_{:K}^{\top}AV_{:K}\Sigma_{+}$ and therefore that
$U_{:K}^{\top}AV_{:K}=\Sigma_{+}$. Combining this with the fact that
$AV_{K:}=0$ then gives

$\displaystyle
U^{\top}AV=\left(\begin{array}[]{cc}U_{:K}^{\top}AV_{:K}&U_{:K}^{\top}AV_{K:}\\\
U_{K:}^{\top}AV_{:K}&U_{K:}^{\top}AV_{K:}\end{array}\right)=\left(\begin{array}[]{cc}\Sigma_{+}&0\\\
0&0\end{array}\right)$

as required.

(ii) $\implies$ (i): For the first equation in Equation 28, left multiply
Equation 29 by $U$ and right-multiply by $\begin{pmatrix}\Sigma_{+}^{-1}&0\\\
0&0\end{pmatrix}\in\mathbb{R}^{q\times q}$ to get

$\displaystyle AV\begin{pmatrix}\Sigma_{+}^{-1}&0\\\
0&0\end{pmatrix}=U\begin{pmatrix}I_{k\times k}&0\\\
0&0\end{pmatrix}=\begin{pmatrix}U_{:K}&0\end{pmatrix}$

then equate the first $K$ columns.

For the second equation in Equation 28, use Equation 29 and orthogonality of
$U$ to compute

$\displaystyle
V^{\top}A^{\top}AV=(U^{\top}AV)^{\top}(U^{\top}AV)=\Sigma_{p\times
q}^{\top}\Sigma_{p\times q}=\Sigma_{q\times q}^{2}$

as required. ∎

######  Corollary 11 (Existence and Uniqueness of SVD).

Let $A\in\mathbb{R}^{p\times q}$ have rank $K$. Then there exist orthogonal
matrices $U\in\mathbb{R}^{p\times p}$ and $V\in\mathbb{R}^{q\times q}$ such
that

$A=U\begin{pmatrix}\Sigma_{+}&0\\\ 0&0\end{pmatrix}V^{\top}=\vcentcolon
U\Sigma V^{\top}$ (31)

where $\Sigma_{+}=\operatorname{diag}(\sigma_{1},\dots,\sigma_{K})$ with
$\sigma_{1}\geq\dots\geq\sigma_{K}>0$. Here, $\sigma_{1},\dots,\sigma_{K}$ are
the square roots of the non-zero eigenvalues of $A^{\top}A$ (or equivalently
$AA^{\top}$) and so uniquely determined.

Moreover, the matrices $U,V$ are unique up to orthogonal transformations
within the spaces of repeated singular values: indeed, suppose we have another
pair of orthogonal matrices $U^{\prime}\in\mathbb{R}^{p\times p}$ and
$V^{\prime}\in\mathbb{R}^{q\times q}$ with $U^{\prime\top}AV^{\prime}=\Sigma$.
Suppose that a non-zero singular value $\sigma>0$ occurs with multiplicity
$\delta$ and let
$V_{\sigma},V_{\sigma}^{\prime},U_{\sigma},U_{\sigma}^{\prime}$ denote the
corresponding columns of $V,V^{\prime},U,U^{\prime}$ respectively (i.e. for
some $k$ we have
$\sigma_{k-1}>\sigma=\sigma_{k}=\dots=\sigma_{k+\delta-1}>\sigma_{k+\delta}$
and $V_{\sigma}=\begin{pmatrix}V_{k}\dots V_{k+\delta-1}\end{pmatrix}$ etc).
Then there is some orthogonal matrix
$O_{\sigma}\in\mathbb{R}^{\delta\times\delta}$ such that

$\displaystyle V_{\sigma}^{\prime}=V_{\sigma}O_{\sigma}\quad\text{and}\quad
U_{\sigma}^{\prime}=U_{\sigma}O_{\sigma}\;.$ (32)

By contrast, for the zero singular value, we can only say that
${\operatorname{span}\\{{U_{K:}}\\}}={\operatorname{span}\\{{U^{\prime}_{K:}}\\}}=\ker(AA^{\top})$,
and
${\operatorname{span}\\{{V_{K:}}\\}}={\operatorname{span}\\{{V^{\prime}_{K:}}\\}}=\ker(A^{\top}A)$
(i.e. we have
$V_{K:}=V^{\prime}_{K:}O_{0}^{(V)},U_{K:}=U^{\prime}_{K:}O_{0}^{(U)}$ but do
not necessarily have $O_{0}^{(V)}=O_{0}^{(U)}$).

###### Proof.

For both existence and uniqueness we use the characterisation from Theorem 10.
For existence, let
$V=\begin{pmatrix}V_{1}&V_{2}\end{pmatrix}\in\mathbb{R}^{q\times q}$ be a