matrix of successive eigenvectors for $A^{\top}A$. Let the corresponding
eigenvalues be
$\sigma_{1}^{2}\geq\dots\geq\sigma_{K}^{2}>0=\sigma_{K+1}^{2}=\dots=\sigma_{q}^{2}$,
and $\sigma_{k}=(\sigma_{k}^{2})^{\nicefrac{{1}}{{2}}}\geq 0$ be their square
roots. Then define diagonal matrices of successive $\sigma_{k}$ as in the
setup of Theorem 10. Then by definition, the second equation of Equation 28
holds.

Further define $U_{:K}=AV_{:K}$ so that the second equation of Equation 28
holds; and finally define $U_{K:}\in\mathbb{R}^{(p-K)\times p}$ be a matrix
whose columns form an orthonormal basis for
${\operatorname{span}\\{{U_{:K}}\\}}^{\perp}$. Then
$U=\begin{pmatrix}U_{:K}&U_{K:}\end{pmatrix}$ formed by concatenation is an
orthogonal matrix, so we can conclude Item (ii) by the equivalence of Theorem
10.

For uniqueness, whether $\sigma$ is non-zero or zero, we have that
$V_{\sigma},V^{\prime}_{\sigma}$ are matrices whose columns span the
$\sigma^{2}$-eigenspace of $A^{\top}A$, so by Proposition 9, there is an
orthogonal matrix $O_{\sigma}$ such that
$V^{\prime}_{\sigma}=V_{\sigma}O_{\sigma}$. If $\sigma>0$ then by the first
equation of Equation 28 we see that this same orthogonal matrix gives
$U^{\prime}_{\sigma}=U_{\sigma}O_{\sigma}$. By contrast, if $\sigma=0$, we can
apply Theorem 10 to $A^{\top}$ with $U,V$ reversed to conclude that
$U_{0},U^{\prime}_{0}$ both span the zero-eigenspace of $AA^{\top}$ and
therefore by Lemma 7 that there is an orthogonal matrix $O_{0}^{(U)}$ such
that $U^{\prime}_{0}=U_{0}O_{0}^{(U)}$ but can no longer conclude that this is
the same orthogonal matrix as transforms $V_{0}$ to $V_{0}^{\prime}$. ∎

######  Definition 5 (SVD).

We call a factorisation of the form Equation 31 from Corollary 11 a (full)
singular value decomposition of the matrix $A$. The $\sigma_{k}$ are called
the singular values of $A$. The successive columns of $U,V$ are called left
and right singular vectors of $A$ respectively; those corresponding to non-
zero singular values are naturally paired by the uniqueness statement of
Corollary 11, and so we may refer to these as pairs of singular vectors, or
simply singular pairs.

Note that we can write a more compact form of Equation 31 by cutting out the
zero components to give

$\displaystyle A$ $\displaystyle=U_{:K}\Sigma_{+}V_{:K}$ (33)

and call this a compact SVD of $A$.

######  Remark 4.

The full SVD is helpful when we need full bases for the spaces
$\mathbb{R}^{p},\mathbb{R}^{q}$ respectively, but it can be more convenient to
consider the compact SVD, which, for example inherits a cleaner notion of
uniqueness from Corollary 11. We will see that one can make a similarly
helpful distinction between full and compact CCA.

####  A.2.2 Sequential variational characterisation

######  Theorem 12 (Variational Characterisation).

Consider a sequence of vectors $(u_{k},v_{k})_{k=1}^{K}$ corresponding to
columns of the matrices $U=(u_{1},\dots,u_{K}),V=(v_{1},\dots,v_{K})$. Then
$U\Sigma V^{\top}$ is a compact SVD for $A$ if and only if
$(u_{k},v_{k})_{k=1}^{K}$ is a sequence of successive solutions to the
programs

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\>&u^{\top}Av\\\
\textrm{subject to}\>&\|u\|_{2}\leq 1,\>\|v\|_{2}\leq 1,\\\
&u^{\top}u_{j}=v^{\top}v_{j}=0\text{ for }1\leq j\leq k-1.\end{split}$ (34)

of maximal length such that each program has strictly positive objective
value.

###### Proof.

Note first that in each case there is some pair attaining the maximum because
we are optimising a continuous objective over a compact set. We will prove
inductively that sequences of maximisers to these programs are precisely
sequences of singular pairs. The $k=1$ case can be seen as a special case of
general induction step but we prove it first as a warm-up.

For $k=1$, a maximiser of Equation 34, must correspond to a maximiser of the
Lagrangian

$\displaystyle\mathcal{L}(u,v;\lambda_{1},\mu_{1})=u^{\top}Av+\tfrac{1}{2}\lambda(1-\lVert
u\rVert^{2}_{2})+\tfrac{1}{2}\mu(1-\lVert v\rVert^{2}_{2})\;.$

This leads to the first order conditions:

$\displaystyle 0=\partial_{u}\mathcal{L}$ $\displaystyle=Av-\lambda u$ (35)
$\displaystyle 0=\partial_{v}\mathcal{L}$ $\displaystyle=A^{\top}u-\mu v$ (36)

for feasible $u,v$ this gives

$\displaystyle\lambda=\lambda\lVert
u\rVert^{2}=u^{\top}Av=v^{\top}A^{\top}u=\mu\lVert v\rVert^{2}=\mu$

so in fact $\lambda=\mu$ and combining the equations we obtain

$\displaystyle A^{\top}Av=A^{\top}\lambda u=\lambda^{2}v$

whence $v$ is an eigenvector of $A^{\top}A$. Moreover, the corresponding
objective value is $u^{\top}Av=\lambda$. So any maximising $u,v$ must be
precisely some top singular pair.

For the general induction step, suppose that
$(u_{1},v_{1}),\dots,(u_{k-1},v_{k-1})$ are some top-$(k-1)$ singular pairs.
The added orthogonality constrains give the larger Lagrangian

$\displaystyle\mathcal{L}(u,v;\lambda_{1},\mu_{1})=u^{\top}Av+\tfrac{1}{2}\lambda(1-\lVert
u\rVert^{2}_{2})+\tfrac{1}{2}\mu(1-\lVert
v\rVert^{2}_{2})+\sum_{j}^{k-1}\lambda_{j}u_{j}^{\top}u+\sum_{j}^{k-1}v_{j}^{\top}v\;.$

This now gives first order conditions

$\displaystyle 0=\partial_{u}\mathcal{L}$ $\displaystyle=Av-\lambda
u+\sum_{j}^{k-1}\lambda_{j}u_{j}$ $\displaystyle 0=\partial_{v}\mathcal{L}$
$\displaystyle=A^{\top}u-\mu v+\sum_{j}^{k-1}\mu_{j}v_{j}.$

Taking the inner product with respect to some $u_{j}$ and using feasibility
shows

$\displaystyle 0=u_{j}^{\top}Av-\lambda
u_{j}^{\top}u+\lambda_{j}=\sigma_{j}v_{j}^{\top}v-0+\lambda_{j}=\lambda_{j}\;.$

Similarly $\mu_{j}=0$. Therefore we recover the previous first order
conditions Equations 35 and 36. So indeed any maximiser must be a singular
pair. Therefore any maximiser must be some choice of top remaining singular
pair, as required. ∎

####  A.2.3 SVD of a bi-linear form

Conveniently, the definition of SVD above can be extended to finite-
dimensional inner product spaces333Or indeed separable Hilbert spaces, see
Carlsson [2021] in a basis independent way.

######  Definition 6 (SVD of a bi-linear form).

Let $\mathcal{X},\mathcal{Y}$ be finite-dimensional real vector spaces of
dimensions $p,q$ respectively. Let
$\mathcal{A}:\mathcal{X},\mathcal{Y}\mapsto\mathbb{R}$ be a bi-linear form.
Then we say that the pair of orthonormal bases
$(z_{k})_{k=1}^{p},(w_{k})_{k=1}^{q}$ give a singular value decomposition of
the operator $\mathcal{A}$ if there exist constants
$\sigma_{1}\geq\dots\geq\sigma_{\min(p,q)}\geq
0=\sigma_{\min(p,q)+1}=\dots=\sigma_{\max(p,q)}$ such that

$\displaystyle\mathcal{A}(z_{k},w_{l})=\delta_{kl}\sigma_{k}\quad\forall\>k\in[p],l\in[q]\;.$

######  Remark 5 (Existence, uniqueness, variational characterisation).

Suppose we have arbitrary orthonormal bases $(x_{i})_{i=1}^{p}$ and
$(y_{i})_{i=1}^{q}$ for $\mathcal{X},\mathcal{Y}$ respectively. Write
$\operatorname{repr}_{(x_{k})_{k}}:\mathcal{X}\mapsto\mathbb{R}^{p}$,
$\operatorname{repr}_{(y_{k})_{k}}:\mathcal{Y}\to\mathbb{R}^{q}$ to be linear
maps giving the representation of vectors in the given bases, and let
$A\in\mathbb{R}^{p\times q}$ be the matrix representing $\mathcal{A}$ in the
given bases. Then for any vectors
$z,z^{\prime}\in\mathcal{X},w,w^{\prime}\in\mathcal{Y}$ we have

$\displaystyle\mathcal{A}(z,w)$
$\displaystyle=\operatorname{repr}_{(x_{k})_{k}}(z)^{\top}\>A\>\operatorname{repr}_{(y_{k})_{k}}(w),$
$\displaystyle\langle z,z^{\prime}\rangle_{\mathcal{X}}$
$\displaystyle=\operatorname{repr}_{(x_{k})_{k}}(z)^{\top}\operatorname{repr}_{(x_{k})_{k}}(z^{\prime}),$
$\displaystyle\langle w,w^{\prime}\rangle_{\mathcal{Y}}$
$\displaystyle=\operatorname{repr}_{(x_{k})_{k}}(w)^{\top}\operatorname{repr}_{(x_{k})_{k}}(w^{\prime})\;,$

by definition. Therefore, the sets of orthonormal vectors
$(z_{k})_{k=1}^{p},(w_{k})_{k=1}^{q}$ give a singular value decomposition for
$\mathcal{A}$ precisely when their vector representations give a singular
value decomposition for the matrix representation of $\mathcal{A}$. We can
therefore deduce notions of existence, uniqueness and a sequential variational
characterisation corresponding to those from the previous two subsections.

###  A.3 Canonical Bases

######  Definition 7 (Canonical bases).

Let $\mathcal{X},\mathcal{Y}$ be finite-dimensional subspaces of some inner
product space $\mathcal{H}$ of dimensions $p,q$ respectively. Consider the bi-
linear form defined by the inner product

$\displaystyle\mathcal{C}:\mathcal{X}\times\mathcal{Y}\mapsto\mathbb{R},\>(x,y)\mapsto\langle
x,y\rangle\;.$

We call a choice of pairs of singular vectors
$(z_{k})_{k=1}^{p},(w_{k})_{k=1}^{q}$ for $\mathcal{C}$ a choice of paired
canonical bases for the pair of subspaces $\mathcal{X},\mathcal{Y}$.

Firstly note that the canonical bases therefore inherit a notion of uniqueness
from the corresponding SVD. Moreover, by Cauchy–Schwarz, for any unit vectors
$z,w\in\mathcal{H}$ we have $\mathcal{C}(z,w)=\langle z,w\rangle\leq\lVert
z\rVert_{\mathcal{H}}\lVert w\rVert_{\mathcal{H}}$ and therefore the singular
values of $\mathcal{C}$ are bounded by 1. This means that they can naturally
be interpreted as cosines of angles.

######  Definition 8 (Canonical angles, $\cos^{2}\Theta$ similarity,
$\sin^{2}\Theta$-discrepancy).

We define canonical angles between the subspaces $\mathcal{X},\mathcal{Y}$ by
$\theta_{k}=\arccos{\sigma_{k}(\mathcal{C})}$ for $k=1,\dots,\min{(p,q)}$.
Define $\cos^{2}\Theta$ similarity

$\quad\cos^{2}\Theta(\mathcal{X},\mathcal{Y})=\sum_{k=1}^{\min{p,q}}\cos^{2}\theta_{k}=\lVert\mathcal{C}\rVert_{HS}^{2}\;.$
(37)

If $\mathcal{X},\mathcal{Y}\subset\mathcal{H}$ each have dimension $K$, then
we write $\Theta=\operatorname{diag}(\theta_{1},\dots,\theta_{K})$ for a
diagonal matrix of canonical angles, and $\sin\Theta,\cos\Theta$ for
corresponding diagonal matrices with $\sin,\cos$ applied element-wise
respectively.

We also use the following shorthand for the $\sin^{2}\Theta$ distance:
$\sin^{2}\Theta(\mathcal{X},\mathcal{Y})=\lVert\sin\Theta(\mathcal{X},\mathcal{Y})\rVert_{F}^{2}$.
We correspondingly recover
$\cos^{2}\Theta(\mathcal{X},\mathcal{Y})=\lVert\cos\Theta(\mathcal{X},\mathcal{Y})\rVert_{F}^{2}$
and the pleasing identity

$\displaystyle\cos^{2}\Theta(\mathcal{X},\mathcal{Y})+\sin^{2}\Theta(\mathcal{X},\mathcal{Y})=K\;.$

######  Remark 6 (Conventions and zero-padding).

There are multiple possible conventions for how many singular values to
define; $\min(p,q)$ works well for most purposes but one could also consider
there being $p,q$ or $\max(p,q)$ singular values by padding with zeros
appropriately. Conveniently, the $\cos^{2}\Theta$ definition by Equation 37
does not depend on how many zero singular values are included.

However, the $\sin^{2}\Theta$ distance does depend on the zero singular
values, because their corresponding $\sin$ is 1; the dimension of
$\Theta,\cos\Theta,\sin\Theta$ could also be defined in multiple reasonable
ways. To avoid potential ambiguity we therefore only consider these objects
when the dimensions of the subspaces are the same.

From this definition, it is easy to deduce the following lemma, which will
help us compute canonical angles in practice.

######  Lemma 13.

Let $X,Y\in\mathbb{R}^{N\times K}$ be matrices with orthonormal columns
spanning the subspaces $\mathcal{X},\mathcal{Y}\subset\mathbb{R}^{N}$. Then
the cosines of the canonical angles between these subspaces are precisely the
singular values of the matrix $X^{\top}Y\in\mathbb{R}^{K\times K}$.

###### Proof.

This follows immediately from Remark 5 upon noting that the covariance
operator is represented by the matrix $X^{\top}Y$ with respect to the
orthonormal bases for $\mathcal{X},\mathcal{Y}$ formed by the columns of $X,Y$
respectively. ∎

######  Remark 7 (Consistency with definition from Stewart and Sun [1990]).

Our definitions are consistent with standard definition of canonical angles
from the literature, and so inherit all standard results. For example, Stewart
and Sun [1990] define canonical angles via the following theorem below. It is
straightforward to prove this theorem and show that the resulting canonical
angles are consistent using the machinery we have developed; however, this is
tedious and therefore omitted.

######  Theorem 14 (Pairs of subspaces, from Stewart and Sun [1990]).

Let $X,Y\in\mathbb{R}^{N\times K}$ with $X^{\top}X=I_{K}$ and
$Y^{\top}Y=I_{K}$. If $2K\leq N$, there are unitary matrices $Q,U$ and $V$
such that

$\displaystyle QXU=\begin{pmatrix}I\\\ 0\\\
0\end{pmatrix}\begin{array}[]{l}\\}K\\\ \\}K\\\ \\}N-2K\end{array},\quad
QYV=\begin{pmatrix}\Gamma\\\ \Sigma\\\ 0\end{pmatrix}\begin{array}[]{l}\\}K\\\
\\}K\\\ \\}N-2K\end{array}$

where $\Gamma=\operatorname{diag}\left(\gamma_{1},\ldots,\gamma_{K}\right)$
and $\Sigma=\operatorname{diag}\left(\sigma_{1},\ldots,\sigma_{K}\right)$ are
such that $\sigma_{1}\geq\ldots\geq\sigma_{K}\geq 0$ and
$\Gamma^{2}+\Sigma^{2}=I_{K}$.

These canonical angles can also be written in terms of the projection matrices
on to the two respective subspaces. The proposition below is precisely theorem
5.5 from Stewart and Sun [1990].

######  Corollary 15.

Let $\mathcal{X}$ and $\mathcal{Y}$ be $K$-dimensional subspaces of
$\mathbb{R}^{N}$ with $2K\leq N$; let
$\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{K}$ be the sines of the
canonical angles between $\mathcal{X}$ and $\mathcal{Y}$. Then

  * •

The singular values of $P_{\mathcal{X}}\left(I-P_{\mathcal{Y}}\right)$ are
$\sigma_{1},\sigma_{2},\ldots,\sigma_{K},0,\ldots,0$.

  * •

The singular values of $P_{\mathcal{X}}-P_{\mathcal{Y}}$ are
$\sigma_{1},\sigma_{1},\sigma_{2},\sigma_{2},\ldots,\sigma_{K},\sigma_{K},0,\ldots,0$.

Note that in each case there are a total of $N$ singular values because the
operators send $\mathbb{R}^{N}\to\mathbb{R}^{N}$.

##  Appendix B CCA foundations: existence, uniqueness and equivalent
formulations

Standard mathematical treatments of CCA, such as Anderson [2003] typically
assume that $\Sigma_{xx},\Sigma_{yy}$ are invertible. Though this is a
reasonable assumption for a population distribution, this means that the
analysis cannot apply to sample CCA in high dimensions when
$\mathbf{C}_{xx},\mathbf{C}_{yy}$ are low rank $N$.

We give a different introduction of CCA that focuses on canonical variates,
and naturally applies both to population CCA and to sample CCA, whether or not
these within-view matrices are invertible. We also carefully specify the
appropriate notion of uniqueness for CCA, which we could not find in previous
texts; this helps provide intuition for what types of statistical estimation
are most natural, see Remark 11.

###  B.1 A general definition of CCA

######  Definition 9 (Full canonical decomposition, general definition).

Let $(x_{1},\dots,x_{p}),(y_{1},\dots,y_{q})$ be two collections of vectors
all from some common vector space $\mathcal{H}$ with inner product
$\langle\cdot,\cdot\rangle$ and corresponding norm $\lVert\cdot\rVert$. Let
the subspaces of vectors they generate be
$\mathcal{X}={\operatorname{span}\\{{x_{j}:j=1,\dots,p}\\}}$,
$\mathcal{Y}={\operatorname{span}\\{{y_{j}:j=1,\dots,q}\\}}$ and have
dimensions $R_{x},R_{y}$ respectively. For
$u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}$ introduce the shorthand

$\displaystyle\mathrm{X}^{\top}u=\sum_{j=1}^{p}u_{j}x_{j},\quad\mathrm{Y}^{\top}v=\sum_{j=1}^{q}v_{j}y_{j}\;.$

Then we say that a pair of sequences of vectors
$(u_{k})_{k=1}^{R_{x}},(v_{k})_{k=1}^{R_{y}}$ give a full basis of canonical
directions for $(x_{j})_{j},(y_{j})_{j}$ if for $k\leq\min(R_{x},R_{y})$ the
directions $u_{k},v_{k}$ solve

$\begin{split}\underset{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\operatorname{maximize}}\,&\langle\mathrm{X}^{\top}u,\mathrm{Y}^{\top}v\rangle\\\
\text{subject to }&\lVert\mathrm{X}^{\top}u\rVert^{2}\leq
1,\>\lVert\mathrm{Y}^{\top}v\rVert^{2}\leq 1,\\\
&\langle\mathrm{X}^{\top}u,\mathrm{X}^{\top}u_{l}\rangle=\langle\mathrm{Y}^{\top}v,\mathrm{Y}^{\top}v_{l}\rangle=0\text{
for }1\leq l\leq k-1.\end{split}$ (38)

and the remaining directions are such that
$(\mathrm{X}^{\top}u_{k})_{k=1}^{R_{x}},(\mathrm{Y}^{\top}v_{k})_{k=1}^{R_{y}}$
form orthonormal bases for $\mathcal{X},\mathcal{Y}$ respectively.

If so, we may concatenate the successive vectors into matrices
$U\in\mathbb{R}^{p\times R_{x}},V\in\mathbb{R}^{q\times R_{y}}$ and
equivalently write any of

$\displaystyle\left((u_{k})_{k=1}^{R_{x}},(v_{k})_{k=1}^{R_{y}}\right)\equiv(U,V)\in\operatorname{CCA}(\mathrm{X},\mathrm{Y})\equiv\operatorname{CCA}\left((x_{j})_{j=1}^{p},(y_{j})_{j=1}^{q}\right)\;.$

######  Remark 8.

Note that one could extend the $(u_{k})_{k},(v_{k})_{k}$ to full bases for
$\mathbb{R}^{p},\mathbb{R}^{q}$ respectively, but in this case it is not
obvious what the correct notion of orthonormality is. We found the formulation
above most convenient to work with, and will not dwell further on this choice.

####  B.1.1 Existence, and uniqueness in variate space

######  Proposition 16 (Equivalence to SVD of covariance operator).

Adopt all setup and notation from Definition 9. Consider the inner product
$\langle\cdot,\cdot\rangle$ as a bi-linear form
$\mathcal{H}\times\mathcal{H}\mapsto\mathbb{R}$. Then canonical correlation
analysis is equivalent to a singular value decomposition of this (inner
product) bi-linear form. Indeed, the canonical correlations are precisely the
singular values of $\langle\cdot,\cdot\rangle$ and

$\displaystyle(U,V)\in\operatorname{CCA}(\mathrm{X},\mathrm{Y})\iff(\mathrm{X}^{\top}U,\mathrm{Y}^{\top}V)\in\operatorname{SVD}(\langle\cdot,\cdot\rangle)\;.$

###### Proof.

This immediately follows by equating the sequential optimisation definition of
CCA given above Equation 38 with the equivalent sequential optimisation
formulation of SVD Theorem 12. ∎

This immediately implies that the canonical correlations are unique. The
canonical variates inherit a notion of uniqueness from the SVD formulation;
this uniqueness is up to orthogonal transformations within subspaces of the
same canonical correlation, as we make precise with the next result.

######  Corollary 17 (Existence and uniqueness).

The canonical correlations $(\rho_{k})_{k}$ are unique, while the canonical
variates $\mathrm{X}^{\top}U,\mathrm{Y}^{\top}V$ inherit the notion of
uniqueness from the corresponding SVD.

In particular, if $U^{\prime},V^{\prime}$ were alternative sets of canonical
directions satisfying Definition 9, and
$\rho=\rho_{k}=\dots=\rho_{k+\delta-1}>0$ have multiplicity $\delta$; then
there exists some orthogonal matrix $O\in\mathbb{R}^{\delta\times\delta}$ such
that

$\displaystyle(\mathrm{X}^{\top}u)^{\prime}_{(k-1)+\epsilon}=\sum_{\eta=1}^{\delta}O_{\eta,\epsilon}\mathrm{X}^{\top}u_{(k-1)+\eta},\quad\text{and}\quad\mathrm{Y}^{\top}v^{\prime}_{(k-1)+\epsilon}=\sum_{\eta=1}^{\delta}O_{\eta,\epsilon}\mathrm{Y}^{\top}v_{(k-1)+\eta}\quad\text{for}\quad\epsilon\in[\delta]\;.$

By contrast, for the zero singular value, we can only say that

$\displaystyle{\operatorname{span}\\{{\mathrm{X}^{\top}u^{\prime}_{K+1},\dots,\mathrm{X}^{\top}u^{\prime}_{p}}\\}}$
$\displaystyle={\operatorname{span}\\{{\mathrm{X}^{\top}u_{K+1},\dots,\mathrm{X}^{\top}u_{p}}\\}}=\mathcal{X}\cap\mathcal{Y}^{\perp}\text{
, and }$
$\displaystyle{\operatorname{span}\\{{\mathrm{Y}^{\top}v^{\prime}_{K+1},\dots,\mathrm{Y}^{\top}v^{\prime}_{p}}\\}}$
$\displaystyle={\operatorname{span}\\{{\mathrm{Y}^{\top}v_{K+1},\dots,\mathrm{Y}^{\top}v_{p}}\\}}=\mathcal{Y}\cap\mathcal{X}^{\perp}\;.$

###### Proof.

This follows immediately by combining the SVD formulation of CCA Proposition
16 with the uniqueness result for SVD Corollary 11. ∎

######  Remark 9 (Canonical bases).

By comparing Proposition 16 with Definition 7 we see that CCA is equivalent to
picking canonical bases in variate space. This provides another way to think
about the break down of sample CCA when $N\leq p$, and specifically when the
columns of $\mathbf{X}$ span $\mathbb{R}^{N}$; then
${\operatorname{span}\\{{\mathbf{Y}}\\}}\subset\mathbb{R}^{N}={\operatorname{span}\\{{\mathbf{X}}\\}}$
so any orthonormal basis for ${\operatorname{span}\\{{\mathbf{Y}}\\}}$ can be
extended to give a pair of canonical bases, and all the canonical angles are
zero.

######  Remark 10 (Additional non-uniqueness in direction space).

Consider the operators
$\mathcal{E}_{x}:u\mapsto\mathrm{X}^{\top}u,\mathcal{E}_{y}:v\mapsto\mathrm{Y}^{\top}v$.
Then by linearity of these operators, we have

$\displaystyle\mathrm{X}^{\top}u^{\prime}=\mathrm{X}^{\top}u\iff
u^{\prime}=u+\nu_{x},\quad\text{where}\quad\nu_{x}\in\ker(\mathcal{E}_{x}),\quad\text{i.e.,}\quad\mathrm{X}^{\top}\nu_{x}=0\;.$
(39)

We can plug this into Corollary 17 to read off the corresponding notions of
uniqueness for the canonical directions. In the subsection below we examine
what the kernel of these operators look like in more detail.

####  B.1.2 Application to population and sample CCA

The formulation above immediately applies to both population CCA and sample
CCA.

For population CCA take $\mathcal{H}\subset\mathcal{L}^{2}$ to be some space
of square integrable random variables with mean zero, and the successive
$x_{j}\vcentcolon=X_{j},y_{j}\vcentcolon=Y_{j}$ for the components of the
random vectors $X,Y$ respectively; then the inner product is precisely the
covariance operator. The additional non-uniqueness in direction space from
Equation 39 is up to vectors in the kernel of $\Sigma_{xx},\Sigma_{yy}$
respectively because

$X^{\top}u=X^{\top}\tilde{u}\text{
a.s.}\quad\iff\quad\mathbb{E}\|X^{\top}(u-\tilde{u})\|^{2}=0\quad\iff\quad(u-\tilde{u})\in\operatorname{Ker}(\Sigma_{xx})\;,$
(40)

and similarly for the $Y$-view.

For sample CCA take
$\mathcal{H}=\\{\mathbf{e}\in\mathbb{R}^{N}:\mathbf{e}^{\top}\mathrm{1}_{{n}}=0\\}$
to be the space of observed sample vectors with mean zero, and the successive
$x_{j}\vcentcolon=\mathbf{X}_{j},y_{j}\vcentcolon=\mathbf{Y}_{j}$ for the
sample vectors corresponding to successive components of the data matrices
$\mathbf{X},\mathbf{Y}$ respectively. Now the inner product is precisely the
sample covariance operator. In this case, one can interpret the non-uniqueness
from Equation 39 as being up to vectors in the kernels of the respective data
matrices, i.e. $\nu_{x}\in\ker(\mathbf{X})$ and similarly for the $Y$-view.

######  Remark 11 (Two types of uniqueness: applied relevance).

Separating out these two different forms of non-uniqueness helps explain our
empirical observations and methodological suggestions.

  * •

The additional non-uniqueness of directions due to low-rank within-view
covariance matrices helps explain why we observe better estimation in variate
space than direction space. This is one reason why we consider analysing
estimated variates and loadings rather than directions. On a more
philosophical level, one could argue from the definition in Equation 38 and
result of Proposition 16 that variates are more fundamental than the
directions.

  * •

The fundamental non-uniqueness up to orthogonal transformations within variate
subspaces helps explain why we observe better estimation with subspace losses
rather than considering separate subspaces. This motivates the need for
registration when comparing different estimates.

The statements above are deliberately blunt, but help convey our key
intuitions. When analysing estimates from some regularised CCA method in
practice, the situation is evidently more complicated. For example, when the
within-view covariances are full rank but poorly conditioned, the canonical
directions are technically unique, but we may expect estimated directions to
be less stable than estimated variates. It may also be helpful conceptually to
consider more subtle behaviour in the context of interactions between these
two different forms of non-uniqueness.

###  B.2 Further properties

####  B.2.1 Matrix SVD representation for CCA

######  Theorem 18.

Suppose $\Sigma_{xx},\Sigma_{yy}$ are both invertible. Then

$\displaystyle(U_{K},V_{K})\in\operatorname{CCA}_{K}(X,Y)\iff(\Sigma_{xx}^{\nicefrac{{1}}{{2}}}U_{K},\Sigma_{yy}^{\nicefrac{{1}}{{2}}}V_{K})\in\operatorname{SVD}(\Sigma_{xx}^{-1/2}\Sigma_{xy}\Sigma_{yy}^{-1/2})$

and the canonical correlations are precisely the singular values.

###### Proof.

To help us, recall our notation for the SVD target matrix from Equation 3 in
the main text:

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

We also now introduce notation for the transformed weights
$\alpha_{k}\vcentcolon=\Sigma_{x}^{1/2}u_{k}$ and
$\beta_{k}\vcentcolon=\Sigma_{y}^{1/2}v_{k}$.

To prove the claim, simply observe

$\displaystyle\operatorname{Cov}\left(u^{\top}X,v^{\top}Y\right)$
$\displaystyle=\>u^{\top}\Sigma_{xy}v$ $\displaystyle=\>\alpha^{\top}T\beta$
$\displaystyle\operatorname{Cov}\left(u^{\top}X,u_{j}^{\top}X\right)$
$\displaystyle=\>u^{\top}\Sigma_{xx}u_{j}^{\top}$
$\displaystyle=\>\alpha^{\top}\alpha_{j}$
$\displaystyle\operatorname{Cov}\left(v^{\top}Y,v_{j}^{\top}Y\right)$
$\displaystyle=\>v^{\top}\Sigma_{yy}v_{j}^{\top}$
$\displaystyle=\>\beta^{\top}\beta_{j}.$

So the optimisation definition of CCA Equation 2 is indeed equivalent to the
variational characterisation in Theorem 12 for finding the singular values of
$T$. ∎

Note that this theorem immediately implies

$\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}$

as we claimed in the main text.

####  B.2.2 Reconstruction and latent variables

One key motivation for using PCA is its interpretation in terms of optimal
reconstruction from low-dimensional transformed variables, or principal
component scores. The scores can therefore be thought of as latent variables.
This latent variable interpretation is particularly clean when the data is
assumed to be multi-variate Gaussian; in this setting it is well known that
PCA can be recovered from a simple Gaussian linear model [Bishop, 2006, Roweis
and Ghahramani, 1999].

CCA can be interpreted analogously in terms of reconstruction from the
canonical variates. As we now show, it is precisely the canonical loading
vectors which determine the reconstruction444Hence their name, drawing off the
factor analysis literature [Gorsuch, 2014].. So to interpret CCA via
reconstruction it is unnecessary to estimate weights, and sufficient to
estimate variates and loadings.

A formula for such reconstruction follows quickly from the orthogonality of
canonical variates. Suppose we obtain bases
$\\{u_{1},\dots,u_{R_{x}}\\},\\{v_{1},\dots,v_{R_{y}}\\}$ for
$\operatorname{range}({\Sigma_{xx}}),\operatorname{range}({\Sigma_{yy}})$
respectively from a full CCA.

Then if we express $X,Y$ in the bases
$\\{\Sigma_{xx}u_{1},\dots,\Sigma_{xx}u_{R_{x}}\\}$,
$\\{\Sigma_{yy}v_{1},\dots,\Sigma_{yy}v_{R_{y}}\\}$ the coefficients
(determined by orthogonality) are precisely the canonical variates; this gives
the expansion

$X=\sum_{k=1}^{R_{x}}\Sigma_{xx}u_{k}\>\langle u_{k},X\rangle\;,\qquad
Y=\sum_{k=1}^{R_{y}}\Sigma_{yy}v_{k}\>\langle v_{k},Y\rangle\;.$ (42)

We can now use the correlation structure of the canonical variates to come up
with more interpretable latent variable expressions; we formalise one
possibility below. Note that this result holds without any distributional
assumptions, but that the latent variables are not necessarily independent,
only uncorrelated.

######  Proposition 19.

There exist orthonormal scalar random variables
$(\xi_{k})_{k=1}^{\max(R_{x},R_{y})},(\epsilon_{k})_{k=1}^{R_{y}}$ such that

$X=\sum_{k=1}^{R_{x}}\Sigma_{xx}u_{k}\>\xi_{k},\qquad
Y=\sum_{k=1}^{R_{y}}\Sigma_{yy}v_{k}\>\left(\rho_{k}\xi_{k}+(1-\rho_{k}^{2})^{1/2}\epsilon_{k}\right)\;.$
(43)

###### Proof.

Write $\xi_{k}=u_{k}^{\top}X$ and define
$\tilde{\epsilon}_{k}=v_{k}^{\top}Y-\xi_{k}\operatorname{Cov}(\xi_{k},v_{k}^{\top}Y)=v_{k}^{\top}Y-\rho_{k}\xi_{k}$.
Then
$\operatorname{Var}(\tilde{\epsilon}_{k})=1-2\rho_{k}^{2}+\rho_{k}^{2}=1-\rho_{k}^{2}$.
So define
$\epsilon_{k}=\tilde{\epsilon}_{k}/(1-\rho_{k}^{2})^{\nicefrac{{1}}{{2}}}$.
Then the $(\xi_{k})_{k},(\epsilon_{k})_{k}$ are orthonormal by construction.
Moreover, rearranging the previous equations gives
$v_{k}^{\top}Y=\rho_{k}\xi_{k}+(1-\rho_{k}^{2})^{1/2}\epsilon_{k}$ so Equation
43 follows from Equation 42. ∎

We conclude by completing the analogy to PCA. The expressions for
reconstruction above also satisfy a notion of optimality, in terms of between-
view covariances, see Section C.2. In addition, CCA can be recovered from a
simple Gaussian linear model [Bach and Jordan, 2005], as we present in Section
F.3.1; note that, as in Equation 42, the canonical variates play a role as
latent variables, and reconstruction is via the canonical loading vectors.

##  Appendix C Registration and biplots

###  C.1 Registration

######  Lemma 20.

Let $Z_{0}\in\mathbb{R}^{n\times k_{0}},Z_{1}\in\mathbb{R}^{n\times k_{2}}$
have orthonormal columns, where $1\leq k_{0}\leq k_{1}$. Let
$Z_{1}^{\top}Z_{0}$ have singular value decomposition $UDV^{\top}$. Then
$S^{*}\vcentcolon=UV^{\top}$ satisfies

$\displaystyle S^{*}\in\underset{S\in\mathbb{R}^{k_{1}\times k_{0}}\text{ s.t.
}S^{\top}S=I_{k_{0}}}{\operatorname*{arg\,min}}\|Z_{1}S-Z_{0}\|_{F}^{2}\;.$

###### Proof.

For arbitrary $S\in\mathbb{R}^{k_{1}\times k_{0}}$ we have

$\displaystyle\|Z_{1}S-Z_{0}\|_{F}^{2}$
$\displaystyle=\|Z_{1}S\|_{F}^{2}+\|Z_{0}\|_{F}^{2}-2\langle
Z_{1}S,Z_{0}\rangle_{F}$
$\displaystyle=\operatorname{trace}S^{\top}Z_{1}^{\top}Z_{1}S+\operatorname{trace}Z_{0}^{\top}Z_{0}-2\operatorname{trace}S^{\top}Z_{1}^{\top}Z_{0}$
$\displaystyle=\operatorname{trace}S^{\top}I_{k_{1}}S+\operatorname{trace}I_{k_{0}}-2\langle
SZ_{1}^{\top}Z_{0}\rangle_{F}$ $\displaystyle=2k_{0}-2\langle
SZ_{1}^{\top}Z_{0}\rangle_{F}\;.$

Next, by von Neumann’s trace inequality [Carlsson, 2021], we have

$\displaystyle\langle SZ_{1}^{\top}Z_{0}\rangle_{F}$
$\displaystyle\leq\sum_{i=1}^{k_{0}}\sigma_{i}(S)\sigma_{i}(Z_{1}^{\top}Z_{0})$
$\displaystyle=\sum_{i=1}^{k_{0}}\sigma_{i}(Z_{1}^{\top}Z_{0})$

where we use that the singular values of a matrix with orthonormal columns are
all 1.

Finally we note that $S^{*}$ as defined in the lemma statement does satisfy
this equality; so maximises $\langle SZ_{1}^{\top}Z_{0}\rangle_{F}$ and
minimises our original objective. ∎

######  Lemma 21.

Let $Z_{0}\in\mathbb{R}^{n\times k_{0}},Z_{1}\in\mathbb{R}^{n\times k_{2}}$ be
such that $Z_{1}$ is of full rank. Then

$\displaystyle\left(Z_{1}^{\top}Z_{1}\right)^{-1}Z_{1}^{\top}Z_{0}=\underset{M\in\mathbb{R}^{k_{1}\times
k_{0}}}{\operatorname*{arg\,min}}\|Z_{1}M-Z_{0}\|_{F}^{2}\;.$

Moreover, if $Z_{0}$ has orthonormal columns then the minimum value is
precisely $\sin^{2}\Theta(Z_{0},Z_{1})$.

###### Proof.

Expanding out the Frobenius norm and completing the square gives

$\displaystyle\|Z_{1}S-Z_{0}\|_{F}^{2}$
$\displaystyle=\operatorname{trace}\left\\{M^{\top}Z_{1}^{\top}Z_{1}M-2M^{\top}Z_{1}^{\top}Z_{0}+Z_{0}^{\top}Z_{0}\right\\}$
$\displaystyle=\operatorname{trace}\left(M-\left(Z_{1}^{\top}Z_{1}\right)^{-1}Z_{1}^{\top}Z_{0}\right)\left(Z_{1}^{\top}Z_{1}\right)\left(M-\left(Z_{1}^{\top}Z_{1}\right)^{-1}Z_{1}^{\top}Z_{0}\right)$
$\displaystyle\qquad-\operatorname{trace}\left(Z_{1}^{\top}Z_{0}\right)^{\top}\left(Z_{1}^{\top}Z_{1}\right)^{-1}\left(Z_{1}^{\top}Z_{0}\right)+\operatorname{trace}Z_{0}^{\top}Z_{0}\;.$

So by full rankness of $Z_{1}$, the quadratic form is strictly convex and the
unique minimum is
$M^{*}\vcentcolon=\left(Z_{1}^{\top}Z_{1}\right)^{-1}Z_{1}^{\top}Z_{0}$ as
stated.

Moreover, the minimum value obtained is

$\displaystyle\lVert Z_{1}M^{*}-Z_{0}\rVert_{F}^{2}$
$\displaystyle=\operatorname{trace}\left\\{Z_{0}Z_{0}^{\top}(I-\mathcal{P}_{Z_{1}})\right\\}\;.$

So when $Z_{0}$ has orthonormal columns, this minimum value is precisely

$\displaystyle\operatorname{trace}\mathcal{P}_{Z_{0}}(I-\mathcal{P}_{Z_{1}})=\lVert\mathcal{P}_{Z_{0}}(I-\mathcal{P}_{Z_{1}})\rVert_{F}^{2}=\sin^{2}\Theta(Z_{0},Z_{1})$

by idempotency of projection matrices and then Corollary 15. ∎

###  C.2 Biplot background

####  C.2.1 Original definition

We now provide further background for the idea of biplots introduced in
Section 7.3. To ease notation, ter Braak [1990] assumed that all the original
variables $X_{i},Y_{i}$ have been standardised to have variance 1. We will
work in this setting for the rest of this subsection. We can replace the
correlations in Definition 3 with covariances. The population version of the
original definition of biplots in ter Braak [1990] is as follows:

######  Definition 10 (Population Covariance Biplot, adapted from ter Braak
[1990]).

Let $U_{K},V_{K}$ be matrices of canonical directions. Respectively plot the
$i^{\text{th}}$ variable in view 1 and $j^{\text{th}}$ variable in view 2 with
coordinates

$\left((\Sigma_{xx}U)_{ik}\right)_{k=1}^{K}\;,\quad\left((\Sigma_{yy}VR)_{jk}\right)_{k=1}^{K}\;.$

This is equivalent to Definition 3 under our assumption of variables with
variance 1. Indeed, Equation 4 implies

$U^{\top}\Sigma_{xy}=U^{\top}\Sigma_{xx}URV^{\top}\Sigma_{yy}=RV^{\top}\Sigma_{yy}$

so in fact

$\displaystyle\rho_{k}\operatorname{Cov}(Y_{j},v_{k}^{\top}Y)=\left(\Sigma_{yy}VR\right)_{jk}=\left(\Sigma_{yx}U\right)_{jk}=\operatorname{Cov}(Y_{j},u_{k}^{\top}X)$
(44)

as required.

######  Remark 12.

ter Braak [1990] also suggest plotting in symmetric coordinates
$\left((\Sigma_{xx}UR^{\nicefrac{{1}}{{2}}})_{ik}\right)_{k=1}^{K}$ and
$\left((\Sigma_{yy}VR^{\nicefrac{{1}}{{2}}})_{jk}\right)_{k=1}^{K}$ but this
makes the interpretation somewhat more difficult so we do not pursue the idea
further.

####  C.2.2 Optimality of Low rank approximation

Proposition 23 states a notion of optimality for biplots; this is a simple
consequence of the powerful Eckhart–Young theorem, which is standard in matrix
analysis [Stewart and Sun, 1990, Bhatia, 1997], stated below for completeness.

######  Theorem 22 (Eckhart–Young).

Let $T\in\mathbb{R}^{p\times q}$. Then $M$ minimises $\|T-M\|_{F}$ over
matrices $M$ of rank at most $K$ if and only if $M=A_{K}R_{K}B_{K}^{\top}$
where $(A_{K},R_{K},B_{K})$ is some top-$K$ SVD of the target $T$.

######  Remark 13.

In fact Eckhart–Young holds for any unitarily invariant norm.

######  Proposition 23 (Optimality of biplots).

Suppose $M^{*}$ maximises the following loss function over $p\times q$
matrices $M$ of rank at most $k$

$\|\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}(\Sigma_{xy}-M)\Sigma_{yy}^{\shortminus\nicefrac{{1}}{{2}}}\|_{F}=\|\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}\Sigma_{xy}\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}-\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}M\Sigma_{yy}^{\shortminus\nicefrac{{1}}{{2}}}\|_{F}$
(45)

then

$M^{*}=\Sigma_{xx}U_{K}R_{K}V_{K}^{\top}\Sigma_{yy}$ (46)

where $(U_{K},R_{K},V_{K})$ is some solution to the top-$K$ CCA problem.

###### Proof.

As before, write
$T=\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}\Sigma_{xy}\Sigma_{yy}^{\shortminus\nicefrac{{1}}{{2}}}$.
Then writing
$\tilde{M}=\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}M\Sigma_{yy}^{\shortminus\nicefrac{{1}}{{2}}}$
the objective is just $\|T-\tilde{M}\|$. So by Eckhart–Young (Theorem 22), we
have $\tilde{M}=A_{K}R_{K}B_{K}^{\top}$ where $(A_{K},R_{K},B_{K})$ are a
top-$K$ SVD for $T$. Now define
$U_{K}=\Sigma_{xx}^{\shortminus\nicefrac{{1}}{{2}}}A_{K},V_{K}=\Sigma_{yy}^{\shortminus\nicefrac{{1}}{{2}}}B_{K}$;
then by the arguments of Section B.2.1, this happens precisely when
$(U_{K},R_{K},V_{K})$ is some solution to the top-$K$ CCA problem. ∎

##  Appendix D Multiple correlations

###  D.1 Valid functions for aggregating correlations

We now present some further background motivating how to aggregate successive
correlations to provide a scalar objective. We show that a broad class of
aggregation methods do in fact recover CCA, and are therefore well motivated.
Oddly, we are not aware of similar analysis or discussion in the existing
literature. Lemmas 26 and 27 may appear to have long statements and proofs; we
reassure the reader that in each case the main content is given in the first
half of the statements and proofs, while the second half contains more careful
(repetitive) analysis of the equality cases.

######  Lemma 24 (Interlacing for SVD).

Let $A\in\mathbb{R}^{m\times n}$ be partitioned as

$\displaystyle A=\begin{pmatrix}A_{11}&A_{12}\\\ A_{21}&A_{22}\end{pmatrix}$

where $A_{11}\in\mathbb{R}^{m^{\prime}\times n^{\prime}}$ and
$m^{\prime}\in[m],n^{\prime}\in[n]$. Then

$\displaystyle\sigma_{i}(A_{11})\leq\sigma_{i}(A)\quad\forall i=1,\dots,n\;.$

###### Proof.

Apply Theorem 4.4 of Stewart and Sun [1990] twice. ∎

We now introduce one additional piece of notation that will be helpful for the
remainder of this section. Let $\mathcal{\rho}^{\text{CCA}}_{{k}}(X,Y)$ denote
the $k^{\text{th}}$ canonical correlation of the pair of random variables
$(X,Y)$. (Note that then $\mathcal{R}^{\text{CCA}}_{[{K}]}(X,Y)$ as defined in
Definition 1 is precisely
$\left(\mathcal{\rho}^{\text{CCA}}_{{k}}(X,Y)\right)_{k=1}^{K}$).

######  Lemma 25 (Interlacing for CCA).

Let $X,Y$ be random vectors taking values in $\mathbb{R}^{p},\mathbb{R}^{q}$
respectively. Suppose we have two arbitrary matrices
$\hat{U}\in\mathbb{R}^{p\times K},\hat{V}\in\mathbb{R}^{q\times K}$. Then we
have the element-wise inequalities

$\mathcal{\rho}^{\text{CCA}}_{{K}}(\hat{U}^{\top}X,\hat{V}^{\top}Y)\leq\mathcal{\rho}^{\text{CCA}}_{{K}}(X,Y)$
(47)

###### Proof.

Take some orthonormal bases $e_{1},\dots,e_{K}$ and $f_{1},\dots,f_{K}$ for
$\operatorname{span}(\hat{U}^{\top}X)$ and
$\operatorname{span}(\hat{V}^{\top}Y)$ respectively. Extend these to
orthonormal bases $e_{1},\dots,e_{p}$ and $f_{1},\dots,f_{q}$ for
$\operatorname{span}(X),\operatorname{span}(Y)$ respectively. Then from the
discussion in Section B, $\operatorname{CCA}(X,Y)$ is precisely the vector of
singular values of the matrix $E^{\top}F$ whose entries are
$(E^{\top}F)_{kl}=\langle e_{k},f_{l}\rangle$.

Moreover, $\mathcal{\rho}^{\text{CCA}}(\hat{U}^{\top}X,\hat{V}^{\top}Y)$ is
precisely the vector of singular values of the top-left $K\times K$ sub-matrix
of $E^{\top}F$. So by lemma 24, the interlacing theorem for singular values,
we do indeed obtain (47). ∎

######  Lemma 26 (Sums of increasing convex functions).

Let $f:\mathbb{R}\to\mathbb{R}$ be a convex function such that $f(x)\geq
f(-x)~{}\forall x>0$. Let $A\in\mathbb{R}^{K\times K}$ be a square matrix.
Then

$\displaystyle\sum_{k=1}^{K}f(A_{kk})\leq\sum_{k=1}^{K}f(\sigma_{k}(A))\;.$
(48)

Moreover, if $f$ is strictly convex then equality implies that there exists
some permutation $\pi\in S_{K}$ such that $|A_{kk}|=\sigma_{\pi(k)}(A)$.
Meanwhile, if $f$ is linear and non-constant then equality implies that there
is an orthogonal matrix $O$ such that
$O^{\top}AO=\operatorname{diag}\sigma(A)$ is a valid SVD of $A$.

###### Proof.

Fix some $f,A$. Define the function
$g:\mathbb{R}^{K}\to\mathbb{R},x\mapsto\sum_{k=1}^{K}f(x_{k})$. Then $g$ can
be viewed as the composition of convex functions and so is itself convex.

Write $a=(A_{kk})_{k=1}^{K}$ for the vector of diagonal entries of $A$. Then
the left hand side of (48) is precisely $g(a)$.

Now take some $b\in\partial g(a)$, the sub-differential of $g$ at $a$. Then
writing $m=g(a)-a^{\top}b$, we have by convexity that $\forall
x\in\mathbb{R}^{K}$

$\displaystyle m+\sum_{k}x_{k}b_{k}=m+x^{\top}b\leq g(x)=\sum_{k}f(x_{k})\;.$
(49)

Now let $B=\operatorname{diag}(b)$ be the diagonal matrix with entries
$B_{kk}=b_{k}$. Then the singular values of $B$ are precisely
$|b_{\pi(1)}|\geq\dots\geq|b_{\pi(K)}|$ where $\pi$ is some permutation of
$\\{1,\dots,K\\}$ that arranges the elements of $b$ into descending order of
modulus.

Then by von Neumann’s trace inequality [Carlsson, 2021] and re-indexing the
sum we have

$\displaystyle
a^{\top}b=\operatorname{trace}(AB)\leq\sum_{k=1}^{K}\sigma_{k}(A)\sigma_{k}(B)=\sum_{k=1}^{K}\sigma_{k}(A)|b_{\pi(k)}|=\sum_{l=1}^{K}\left(\sigma_{\pi^{-1}(l)}(A)\,\operatorname{sign}(b_{l})\right)b_{l}.$
(50)

Combining the above

$\displaystyle\sum_{k=1}^{K}f(A_{kk})$ $\displaystyle=g(a)$
$\displaystyle=m+a^{\top}b$ $\displaystyle=m+\operatorname{trace}(AB)$
$\displaystyle\leq
m+\sum_{l=1}^{K}\left(\sigma_{\pi^{-1}(l)}(A)\,\operatorname{sign}(b_{l})\right)b_{l}$
(51)
$\displaystyle\leq\sum_{l=1}^{K}f\left(\sigma_{\pi^{-1}(l)}(A)\,\operatorname{sign}(b_{l})\right)$
(52) $\displaystyle\leq\sum_{l=1}^{K}f\left(\sigma_{\pi^{-1}(l)}(A)\right)$
(53) $\displaystyle=\sum_{k=1}^{K}f\left(\sigma_{k}(A)\right)$

where the first inequality was from (49), the second inequality from (50), and
the third from the assumption that $f(x)\geq f(-x)~{}\forall x>0$.

For the moreover conclusions, note that in either case, equality in the
theorem statement implies equality in all the inequalities in the chain of
inequalities above. First consider the case where $f$ is strictly convex. Then
equality holds in Equation 49 if and only if $x=a$. This was applied at
Equation 52 to the vector $x$ with entries
$x_{l}=\sigma_{\pi^{-1}(l)}(A)\,\operatorname{sign}(b_{l})$. Equating with
$a_{l}=A_{ll}$ and taking the absolute value gives the claim.

The condition that $f$ is linear and non-constant, combined with $f(x)\geq
f(-x)$ for $x>0$ means that $f$ must have strictly positive slope. By
subtracting the constant and rescaling we may therefore assume that $f(x)=x$.
In this case, $\partial g(a)=\\{\mathrm{1}_{{K}}\\}$ so we have
$b=\mathrm{1}_{{K}}$ and $B$ is precisely the identity matrix. Singular value
decompositions of $B$ are precisely $I=OIO^{\top}$ for an orthogonal matrix
$O\in\mathbb{R}^{K\times K}$. But then the equality case of von Neumann’s
trace inequality from Carlsson [2021] states that $A$ and $B$ ‘share singular
vectors’, which in this case implies that there is some orthogonal matrix $O$
such that $O^{\top}AO=\operatorname{diag}\sigma(A)$ is a valid SVD. ∎

Finally we combine the results above to show that CCA can indeed be defined by
maximising sums of convex functions of canonical correlations.

######  Proposition 27.

Let $f:\mathbb{R}\to\mathbb{R}$ be a convex function such that $f(x)\geq
f(-x)~{}\forall x>0$. Let $X,Y$ be random vectors taking values in
$\mathbb{R}^{p},\mathbb{R}^{q}$ respectively. For
$K\in\\{1,\dots,\min(p,q)\\}$ define functions

$\displaystyle g_{K}:\mathbb{R}^{p\times K}\times\mathbb{R}^{q\times
K}\to\mathbb{R},\>(U,V)\mapsto\sum_{k=1}^{K}f\left((\rho^{\text{oracle}}(u_{k},v_{k})\right)$

where $U,V$ have columns $u_{k},v_{k}$ respectively.

Then

$\displaystyle\operatorname{CCA}(X,Y)\subset\underset{\begin{subarray}{c}U\in\mathbb{R}^{p\times
K},V\in\mathbb{R}^{q\times K}\\\
U^{\top}\Sigma_{xx}U=V^{\top}\Sigma_{yy}V=I_{K}\end{subarray}}{\operatorname*{arg\,max}}{\>\sum_{k=1}^{K}f\left((\rho^{\text{oracle}}(u_{k},v_{k})\right)}=\vcentcolon\mathcal{G}_{K}\;.$
(54)

Moreover, if $f$ is linear and non-constant then $\mathcal{G}_{K}$ consists
precisely of matrices $U,V$ whose column span forms a top-$K$ subspace for the
CCA problem. i.e.

$\displaystyle\mathcal{G}_{K}=\left\\{(UO,VO):\>O\in\mathcal{O}(K),(U,V)\in\operatorname{CCA}_{K}(X,Y)\right\\}.$

If however, $f$ is strictly convex and $f(x)=f(-x)\forall x\in\mathbb{R}$,
then $\mathcal{G}_{K}$ consists precisely of matrices $U,V$ whose columns can
be reordered and possibly negated to form a solution to the
$\operatorname{CCA}$ problem, i.e.

$\displaystyle\mathcal{G}_{K}=\left\\{(UP,VPS):\right.$
$\displaystyle\left.(U,V)\in\operatorname{CCA}_{K}(X,Y),\right.$
$\displaystyle\left.P\text{ is a permutation matrix and }S\text{ is diagonal
with entries in }\pm 1\right\\}\;.$

###### Proof.

First note that if $U^{*},V^{*}\in\operatorname{CCA}_{K}(X,Y)$ then
$g_{K}(U^{*},V^{*})=\sum_{k=1}^{K}f(\rho_{k}^{*})$ where
$\rho_{k}^{*}\vcentcolon=\mathcal{\rho}^{\text{CCA}}_{{k}}(X,Y)$. So it is
sufficient to show that this is the maximum possible value of $g_{K}$ over
feasible $(U,V)\in\mathbb{R}^{p\times K}\times\mathbb{R}^{q\times K}$, i.e.
those for which $U^{\top}\Sigma_{xx}U=I,V^{\top}\Sigma_{yy}V=I$.

So now take arbitrary feasible $U,V$. Define the matrix

$\displaystyle M=U^{\top}\Sigma_{xy}V\in\mathbb{R}^{K\times K}\;.$

Then by the orthonormality constraint, in fact we have
$M_{kk}=\operatorname{Corr}(u_{k}^{\top}X,v_{k}^{\top}Y)=\rho^{\text{oracle}}(u_{k},v_{k})$.
So the objective in (54) is precisely $\sum_{k=1}^{K}f(M_{kk})$.

The orthonormality also means that the the canonical correlations
$\mathcal{\rho}^{\text{CCA}}_{{k}}(U^{\top}X,V^{\top}Y)$ are precisely the
singular values $\sigma_{k}(M)$ of $M$. So by Lemma 25 we have that
$\sigma_{k}(M)\leq\rho_{k}^{*}$.

So by Lemma 26 we have

$\displaystyle
g_{K}(U,V)=\sum_{k=1}^{K}f(M_{kk})\leq\sum_{k=1}^{K}f(\sigma_{k}(M))\leq\sum_{k=1}^{K}f(\rho_{k}^{*})$

as required.

Suppose that $f$ is linear and non-constant, and $(U,V)$ attain the maximal
value of $\sum_{k=1}^{K}f(\rho_{k}^{*})$. Then both inequalities in the final
equation must be equalities. The latter gives $\sigma(M)=\rho$. The former
implies that we must have $O^{\top}MO=\operatorname{diag}(\sigma(M))$ for some
orthogonal matrix $O\in\mathbb{R}^{K\times K}$, by the corresponding equality
case of Lemma 26. Write $U^{*}=UO,V^{*}=V0$. Then substituting the two results
above gives
${U^{*}}^{\top}\Sigma_{xy}V^{*}=O^{\top}MO=\operatorname{diag}(\rho)$ and we
still have the required orthonormality
${U^{*}}^{\top}\Sigma_{xx}U^{*}={V^{*}}^{\top}\Sigma_{yy}V^{*}=O^{\top}O=I_{K}$.
Therefore $U^{*},V^{*}$ are indeed top-$K$ CCA matrices, as required. For the
reverse inclusion, again by removing the intercept and re-scaling we may WLOG
assume that $f(x)=x$. Take any $U=U^{*}O,V=V^{*}O$ with
$(U^{*},V^{*})\in\operatorname{CCA}_{K}(X,Y)$. Then
$\sum_{k=1}^{K}\rho^{\text{oracle}}(u_{k},v_{k})=\operatorname{trace}U^{\top}\Sigma_{xy}V=\operatorname{trace}O^{\top}{U^{*}}^{\top}\Sigma_{xy}V^{*}O=\operatorname{trace}{U^{*}}^{\top}\Sigma_{xy}V^{*}=\sum_{k=1}^{K}\rho_{k}^{*}$
does indeed attain the maximum value.

Finally, consider $f$ strictly convex with $f(x)=f(-x)\forall x\in\mathbb{R}$.
Suppose $(U,V)$ attain the maximal value of $\sum_{k=1}^{K}f(\rho_{k}^{*})$.
Then again we must have $\sigma(M)=\rho$. We must also have equality in the
application of Lemma 26 and therefore there is some permutation $\pi$ such
that $|M_{kk}|=\sigma_{\pi(k)}(M)$. So combining we get
$|M_{kk}|=\rho_{\pi(k)}$.

Now define new matrices $U^{*},V^{*}$ by their columns

$\displaystyle
U^{*}_{k}=U_{\pi^{-1}(k)},\>V^{*}_{k}=V_{\pi^{-1}(k)}\operatorname{sign}(M_{\pi^{-1}(k)\pi^{-1}(k)})\;.$

Then $U^{*},V^{*}$ are both feasible, and satisfy
${U_{k}^{*}}^{\top}\Sigma_{xy}V^{*}_{k}=|M_{\pi^{-1}(k)\pi^{-1}(k)}|=\rho_{\pi\pi^{-1}(k)}^{*}=\rho_{k}^{*}$
so are indeed top-$K$ CCA matrices. It is straightforward to rewrite
$U^{*},V^{*}$ into the matrix form in the proposition statement. For the
reverse inclusion, take any $U^{*},V^{*}$ in this matrix form, then by direct
computation one can verify that the maximal value is obtained. ∎

##  Appendix E Graphical CCA

###  E.1 Sparse precision gives sparse directions

See 4

###### Proof.

We start with the Schur complement identity that
$\Omega_{xx}^{-1}\Omega_{xy}=-\Sigma_{xy}\Sigma_{yy}^{-1}$. Substituting this
into Equation 4 gives

$\displaystyle UDV^{\top}$
$\displaystyle=\Sigma_{xx}^{-1}\Sigma_{xy}\Sigma_{yy}^{-1}$
$\displaystyle=\Sigma_{xx}^{-1}\Omega_{xx}^{-1}\Omega_{xy}$
$\displaystyle=(\Omega_{xx}-\Omega_{xy}\Omega_{yy}^{-1}\Omega_{yx})\Omega_{xx}^{-1}\Omega_{xy}$
$\displaystyle=\Omega_{xy}(I-\Omega_{yy}^{-1}\Omega_{yx}\Omega_{xx}^{-1}\Omega_{xy}).$

Hence

$\displaystyle
U=(UDV^{\top})\Sigma_{yy}VD^{-1}=\Omega_{xy}(I-\Omega_{yy}^{-1}\Omega_{yx}\Omega_{xx}^{-1}\Omega_{xy})\Sigma_{yy}VD^{-1}$

and so $U$ inherits zero-rows from $\Omega_{xy}$: for any $a\in A$,

$\displaystyle
U_{a,:}=({\Omega_{xy}})_{a,:}(I-\Omega_{yy}^{-1}\Omega_{yx}\Omega_{xx}^{-1}\Omega_{xy})\Sigma_{yy}VD^{-1}=0\;,$

as required. ∎

###  E.2 Relationship to partial correlations

Partial correlations provide an alternative lens through which to investigate
the relationship between sparse precision and sparse canonical directions.
These are often used in the statistics literature, and may help the reader
develop complementary intuition. We first give a definition and well known
result relating partial correlations to the precision matrix555Some authors
even define partial correlations this way, e.g. Lauritzen [1996]., before
applying this result to CCA.

######  Definition 11 (Partial Correlation, [Prokhorov, ]).

Let $X,Y$ be scalar random variables, and $Z$ be a random vector in
$\mathbb{R}^{d}$, each with mean zero. Let the coefficients for the optimal
linear regression fit, and the corresponding residuals be

$\displaystyle\beta_{X}$
$\displaystyle=\operatorname*{arg\,min}_{\beta\in\mathbb{R}^{d}}\mathbb{E}\lVert
X-\beta^{\top}Z\rVert^{2},$ $\displaystyle R_{X}$
$\displaystyle=X-\beta_{X}^{\top}Z,$ $\displaystyle\beta_{Y}$
$\displaystyle=\operatorname*{arg\,min}_{\beta\in\mathbb{R}^{d}}\mathbb{E}\lVert
Y-\beta^{\top}Z\rVert^{2},$ $\displaystyle R_{Y}$
$\displaystyle=Y-\beta_{Y}^{\top}Z\;.$

Then the partial correlation is defined by

$\displaystyle\operatorname{Corr}\left(X,Y|Z\right)=\frac{\operatorname{Cov}(R_{X},R_{Y})}{\operatorname{Var}(R_{X})^{\nicefrac{{1}}{{2}}}\operatorname{Var}(R_{Y})^{\nicefrac{{1}}{{2}}}}\;.$
(55)

For $X,Y,Z$ with arbitrary mean, one can define
$\operatorname{Corr}\left(X,Y|Z\right)=\operatorname{Corr}\left(X-\mathbb{E}X,Y-\mathbb{E}Y|Z-\mathbb{E}Z\right)$.

One can also extend this immediately to arbitrary random vectors
$X\in\mathbb{R}^{p},Y\in\mathbb{R}^{q}$ to give the matrix of partial
correlations

$\displaystyle\operatorname{Corr}\left(X,Y\mid
Z\right)=\left(\operatorname{Corr}\left(X_{i},Y_{j}\mid
Z\right)\right)_{i\in[p],j\in[q]}\;.$

######  Lemma 28 (Partial correlations from precision matrix).

Let $Z$ be a full-rank random vector in $\mathbb{R}^{d}$ with precision matrix
$\Omega$. Then for any $i,j\in[d]$ distinct

$\displaystyle\operatorname{Corr}\left(Z_{i},Z_{j}\mid
Z_{(-ij)}\right)=\frac{\Omega_{ij}}{\Omega_{ii}^{\nicefrac{{1}}{{2}}}\Omega_{jj}^{\nicefrac{{1}}{{2}}}}\;.$

###### Alternative proof of Proposition 4 using variational characterisation
of CCA and partial correlations.

Write $B=[p]\ A$. By Lemma 28, zeroes in the precision matrix imply zero
partial correlations. Therefore our main assumption becomes
$\operatorname{Corr}\left(X_{A};Y\mid X_{B}\right)=0$.

Now define the (multivariate) linear regression coefficients and residuals

$\displaystyle\Gamma_{A}$
$\displaystyle=\operatorname*{arg\,min}_{\Gamma\in\mathbb{R}^{|B|\times|A|}}\mathbb{E}\lVert
X_{A}-\Gamma^{\top}X_{B}\rVert^{2},$ $\displaystyle R_{A}$
$\displaystyle=X_{A}-\Gamma_{A}^{\top}X_{B},$ $\displaystyle\Gamma_{Y}$
$\displaystyle=\operatorname*{arg\,min}_{\Gamma\in\mathbb{R}^{|B|\times
q}}\mathbb{E}\lVert Y-\Gamma^{\top}X_{B}\rVert^{2},$ $\displaystyle R_{Y}$
$\displaystyle=Y-\Gamma_{Y}^{\top}X_{B}\;.$

By properties of linear regression we have
$\operatorname{Cov}(R_{A},X_{B})=\operatorname{Cov}(R_{Y},X_{B})=0$. We can
also now write the main assumption as $\operatorname{Cov}(R_{A},R_{Y})=0$. Now
take an arbitrary pair $u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}$. Then from the
definition of $R_{A}$ we have
$u_{A}^{\top}X_{A}=u_{A}^{\top}\Gamma_{A}^{\top}X_{B}+u_{A}^{\top}R_{A}$, and
therefore

$\displaystyle
u^{\top}X=u_{A}^{\top}X_{A}+u_{B}^{\top}X_{B}=(u_{B}+\Gamma_{A}u_{A})^{\top}X_{B}+u_{A}^{\top}R_{A}\;.$
(56)

So now define $\tilde{u}\in\mathbb{R}^{p}$ by
$\tilde{u}_{B}=u_{B}+\Gamma_{A}u_{A},\tilde{u}_{A}=0$. Then since
$\operatorname{Cov}(R_{A},X_{B})=0$ we have

$\displaystyle\operatorname{Var}(u^{\top}X)=\operatorname{Var}(\tilde{u}_{B}^{\top}X_{B})+\operatorname{Var}(u_{A}^{\top}R_{A})\geq\operatorname{Var}(\tilde{u}_{B}^{\top}X_{B})=\operatorname{Var}(\tilde{u}^{\top}X)$
(57)

with equality if and only if $u_{A}^{\top}R_{A}=0$, i.e. that
$u^{\top}X=\tilde{u}^{\top}X$.

We can also use the definition of residuals and orthogonality of $R_{A}$ to
both $X_{B}$ and $R_{Y}$ to compute

$\displaystyle\operatorname{Cov}(u^{\top}X,v^{\top}Y)-\operatorname{Cov}(\tilde{u}^{\top}X,v^{\top}Y)=\operatorname{Cov}(u_{A}^{\top}R_{A},v^{\top}Y)=\operatorname{Cov}(u_{A}^{\top}R_{A},v^{\top}\Gamma^{\top}X_{B}+v^{\top}R_{Y})=0.$

and conclude that the projected variables have the same covariances with
$v^{\top}Y$ as the original variables. They also have the same covariances
with $X_{B}$ as the original variables by a similar computation

$\displaystyle\operatorname{Cov}(u^{\top}X,X_{B})-\operatorname{Cov}(\tilde{u}^{\top}X,X_{B})=\operatorname{Cov}(u_{A}^{\top}R_{A},X_{B})=0\;.$
(58)

Now let $(X^{\top}u_{k},Y^{\top}v_{k})_{k=1}^{K}$ be some sequence of
canonical variates, as defined by the program of Equation 38. Inductively
apply the construction above to each pair of directions $u_{k},v_{k}$ to yield
the sequence of vectors $\tilde{u}_{k}$ in $\mathbb{R}^{p}$. By the induction
hypothesis, $u_{l}^{\top}X$ can be written only in terms of variables in $B$
for each $l<k$; so by Equation 58, $\tilde{u}_{k}^{\top}X$ satisfies the
required orthogonality to these previous canonical variates. Since the
canonical variates maximise the correlation
$\operatorname{Corr}(u_{k}^{\top}X,v_{k}^{\top}Y)$, we must have equality in
Equation 57 and therefore that $u_{k}^{\top}X=\tilde{u}_{k}^{\top}X$. Finally,
we note that $\tilde{u}_{k}$ has support in $B$, as required. ∎

###  E.3 Lipschitz Target Matrix

First we recall the assumption from the main text.

See A

We actually prove a slightly stronger version of the ‘Lipschitz plug-in’
theorem from the main text, with an additional bound that will be useful in
Section E.4. The proof spans the remainder of this section.

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

Assume Assumption A$(M,m)$. Then there exists constants
$C_{1},C_{2}\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_{1}\lVert\Omega-\hat{\Omega}\rVert_{\textrm{op}},\quad\text{and}\quad\lVert\hat{\Sigma}_{xx}^{-1}-\Sigma_{xx}^{{-1}}\rVert_{\textrm{op}}\leq
C_{2}\lVert\Omega-\hat{\Omega}\rVert_{\textrm{op}}\;.$ (59)

###### Proof.

The main idea is to show that the function $f:\Omega\mapsto T$ is Lipschitz by
decomposing it as sums, products and compositions of Lipschitz functions.

We shall write $\|\Omega-\hat{\Omega}\|=\epsilon$ for the remainder of the
proof; so we seek to bound $\|\hat{T}-T\|\leq C_{1}\epsilon$, and will track
the constants in front of epsilon over the course of the proof. Note also that
the only matrix norm needed for this proof is the operator norm, so we drop
the subscript-op from $\lVert\cdot\rVert_{\text{op}}$ for the remainder of the
proof. All vector norms refer to the Euclidean norm. In this section, we will
sometimes denote diagonal blocks with a single subscript to avoid clutter
(e.g., writing $\Sigma_{x}$ for $\Sigma_{xx}$).

####  E.3.1 Identities from Schur complements

Firstly, let’s introduce notation for the following Schur complements:

$\displaystyle
S_{x}\vcentcolon=\Omega_{x}-\Omega_{xy}\Omega_{y}^{-1}\Omega_{yx}=\Sigma_{x}^{-1}$
$\displaystyle
S_{y}\vcentcolon=\Omega_{y}-\Omega_{yx}\Omega_{x}^{-1}\Omega_{xy}=\Sigma_{y}^{-1}$
$\displaystyle\Psi_{x}\vcentcolon=\Sigma_{x}-\Sigma_{xy}\Sigma_{y}^{-1}\Sigma_{yx}=\Omega_{x}^{-1}$
$\displaystyle\Psi_{y}\vcentcolon=\Sigma_{y}-\Sigma_{yx}\Sigma_{x}^{-1}\Sigma_{xy}=\Omega_{y}^{-1}.$

Next recall the identity

$\Sigma_{xy}=\Sigma_{x}^{1/2}T\Sigma_{y}^{1/2}=\Sigma_{x}^{1/2}ADB^{\top}\Sigma_{y}^{1/2}=\Sigma_{x}UDV^{\top}\Sigma_{y}\;.$
(60)

We can now rewrite

$\Psi_{x}=\Sigma_{x}-\Sigma_{x}^{1/2}TT^{\top}\Sigma_{x}^{1/2}=\Sigma_{x}^{1/2}(I-TT^{\top})\Sigma_{x}^{1/2}=\Sigma_{x}^{1/2}(I-AD^{2}A^{\top})\Sigma_{x}^{1/2}\;.$
(61)

Starting with identity for off-diagonal inverse blocks from Schur’s formula we
also get

$\displaystyle\Omega_{xy}$
$\displaystyle=-\Psi_{x}^{-1}\Sigma_{xy}\Sigma_{y}^{-1}$
$\displaystyle=-\Sigma_{x}^{-1/2}(I-TT^{\top})^{-1}\Sigma_{x}^{-1/2}\cdot\Sigma_{x}^{1/2}T\Sigma_{y}^{1/2}\cdot\Sigma_{y}^{-1}$
$\displaystyle=-\Sigma_{x}^{-1/2}(I-TT^{\top})^{-1}T\Sigma_{y}^{-1/2}$
$\displaystyle=-\Sigma_{x}^{-1/2}A\frac{D}{(I-D^{2})}B^{\top}\Sigma_{y}^{-1/2}$
(62)

where we understand that $D$ is diagonal with entries in $(0,1)$ so the
inverse exists and we stretch notation to write
$\frac{D}{(I-D^{2})}=D(I-D^{2})^{-1}$.

####  E.3.2 Initial operator norm bounds

Firstly note that $\|S_{x}\|=\|\Sigma_{x}^{-1}\|\leq m$, and similarly
$\|S_{y}\|\leq m$. Furthermore, $0\preceq\Psi_{x}\preceq\Sigma_{x}$, therefore

$\|\Psi_{x}\|\leq\|\Sigma_{x}\|\leq M\;,$ (63)

and similarly, $\|\Psi_{y}\|\leq M$. Next by (61)

$\|\Omega_{x}\|\leq\|\Sigma_{x}^{-1/2}\|\|A(I-D^{2})^{-1}A^{\top}\|\|\Sigma_{x}^{-1/2}\|\leq\frac{1}{1-\rho_{1}^{2}}m\;,$
(64)

and similarly $\|\Omega_{y}\|\leq\frac{1}{1-\rho_{1}^{2}}m$. Next by (60) we
have

$\|\Sigma_{xy}\|\leq\|\Sigma_{x}^{1/2}\|\|ADB^{\top}\|\|\Sigma_{y}^{1/2}\|\leq\rho_{1}M\;,$
(65)

and by (62) we have

$\|\Omega_{xy}\|\leq\|\Sigma_{x}^{-1/2}\|\left\|A\frac{D}{(I-D^{2})}B^{\top}\right\|\|\Sigma_{y}^{-1/2}\|\leq\frac{\rho_{1}}{1-\rho_{1}^{2}}m=\gamma_{1}m\;,$
(66)

where we define $\gamma_{1}=\frac{\rho_{1}}{1-\rho_{1}^{2}}$ for later
convenience.

####  E.3.3 Some useful lemmas

######  Lemma 30.

$\|\Sigma\|\leq 2M$.

###### Proof.

Consider a vector $\begin{pmatrix}u\\\ v\end{pmatrix}\in\mathbb{R}^{p+q}$.
Then

$\begin{pmatrix}u\\\ v\end{pmatrix}^{\top}\Sigma\begin{pmatrix}u\\\
v\end{pmatrix}\leq\begin{pmatrix}u\\\
v\end{pmatrix}^{\top}\Sigma\begin{pmatrix}u\\\
v\end{pmatrix}+\begin{pmatrix}u\\\ \shortminus
v\end{pmatrix}^{\top}\Sigma\begin{pmatrix}u\\\ \shortminus
v\end{pmatrix}=2\left(u^{\top}\Sigma_{x}u+v^{\top}\Sigma_{y}v\right)\leq
2M\left\|\begin{pmatrix}u\\\ v\end{pmatrix}\right\|^{2}\;.\hfill\qed$

######  Lemma 31.

Restriction is 1-Lipschitz with respect to the operator norm; in particular,
$\|\hat{\Sigma}_{xy}-\Sigma_{xy}\|\leq\|\hat{\Sigma}-\Sigma\|$.

###### Proof.

Consider an arbitrary block-structured matrix $\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\in\mathbb{R}^{(p+q)\times(p+q)}$ and $u\in\mathbb{R}^{p}$. We
have

$\|Cu\|\leq\left\|\begin{pmatrix}A\\\
C\end{pmatrix}u\right\|=\left\|\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\begin{pmatrix}u\\\
0\end{pmatrix}\right\|\leq\left\|\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\right\|\|u\|\;.$

So indeed we have $\|C\|\leq\left\|\begin{pmatrix}A&B\\\
C&D\end{pmatrix}\right\|$ and the result follows by linearity of restriction.
∎

######  Lemma 32.

Let $A,B\in\mathbb{R}^{n\times n}$ with $\|B^{-1}\|\leq M$ and
$\|A-B\|\leq\epsilon$. Then
$\|A^{-1}-B^{-1}\|\leq\frac{M^{2}\epsilon}{1-M\epsilon}$. Furthermore
$\|A^{-1}\|\leq\frac{M}{1-M\epsilon}$.

###### Proof.

Write $\delta=\|A^{-1}-B^{-1}\|$. Then

$\begin{split}\delta=\|A^{-1}-B^{-1}\|=\|A^{-1}(B-A)B^{-1}\|\leq\|A^{-1}\|&\|A-B\|\|B^{-1}\|\\\
\leq(\|B^{-1}\|+\|A^{-1}-B^{-1}\|)&\|A-B\|\|B^{-1}\|\leq(M+\delta)\epsilon
M\;.\end{split}$

So rearranging gives $\delta\leq\frac{M^{2}\epsilon}{1-M\epsilon}$ and so also

$\|A^{-1}\|\leq\|B^{-1}\|+\|A^{-1}-B^{-1}\|\leq
M+\delta\leq\frac{M(1-M\epsilon)+M^{2}\epsilon}{(1-M\epsilon)}=\frac{M}{1-M\epsilon}\;.\hfill\qed$

######  Lemma 33.

Let $A,B\in\mathbb{R}^{n\times n}$ be positive definite then

$\|\sqrt{A}-\sqrt{B}\|\leq\frac{\|A-B\|}{\lambda_{\text{min}}(\sqrt{A}+\sqrt{B})}\leq\frac{\|A-B\|}{\sqrt{\lambda_{\text{min}}(A)}+\sqrt{\lambda_{\text{min}}(B)}}\;.$
(67)

###### Proof.

Let $x\in\mathbb{R}^{n}$ be a unit-norm eigenvector of $\sqrt{A}-\sqrt{B}$
with eigenvalue $\|\sqrt{A}-\sqrt{B}\|=|\mu|$; by switching $A,B$, we can
assume WLOG that $\mu\geq 0$ then

$\displaystyle\|A-B\|$ $\displaystyle\geq x^{\top}(A-B)x$
$\displaystyle=x^{\top}\\{(\sqrt{A}-\sqrt{B})\sqrt{A}+(\sqrt{A}-\sqrt{B})\sqrt{B}\\}x$
$\displaystyle=\mu x^{\top}(\sqrt{A}+\sqrt{B})x$
$\displaystyle\geq\mu\lambda_{\text{min}}(\sqrt{A}+\sqrt{B})$
$\displaystyle\geq\mu\lambda_{\text{min}}(\sqrt{A})+\lambda_{\text{min}}(\sqrt{B})\;,$

where the first equality follows from rearranging, and the final two
inequalities follow from the variational characterisation of eigenvalues. ∎

####  E.3.4 Further operator norm bounds

For clarity we shall put these in bullet form:

  * •

Recall that Equation 61 gave $\|\Omega_{x}^{-1}\|=\|\Psi_{x}\|\leq M$; and
because restriction is Lipschitz we have
$\|\hat{\Omega}_{x}-\Omega_{x}\|\leq\|\hat{\Omega}-\Omega\|\vcentcolon=\epsilon$.
Applying these to Lemma 32 gives

$\|\hat{\Omega}_{x}\|\leq\frac{M}{1-M\epsilon},\quad\text{and}\quad\|\hat{\Omega}_{x}^{-1}-\Omega_{x}^{-1}\|\leq\frac{M^{2}}{1-M\epsilon}\epsilon$
(68)

and similarly

$\|\hat{\Omega}_{y}\|\leq\frac{M}{1-M\epsilon},\quad\text{and}\quad\|\hat{\Omega}_{y}^{-1}-\Omega_{y}^{-1}\|\leq\frac{M^{2}}{1-M\epsilon}\epsilon\;.$
(69)

  * •

By combining the triangle inequality with Lemma 31 and (66) we obtain

$\|\hat{\Omega}_{xy}\|\leq\epsilon+\|\Omega_{xy}\|\leq\epsilon+\gamma_{1}m\;.$
(70)

  * •

Combining Lemmas 30 and 32 gives

$\|\hat{\Sigma}-\Sigma\|\leq\frac{(2M)^{2}}{1-2M\epsilon}\epsilon=\vcentcolon\kappa\epsilon\;.$
(71)

where we defined $\kappa=\kappa(\epsilon)=\frac{(2M)^{2}}{1-2M\epsilon}$. Then
by Lemma 31 also

$\|\hat{\Sigma}_{xy}-\Sigma_{xy}\|\leq\kappa\epsilon,\quad\text{and}\quad\lVert\hat{\Sigma}_{xx}-\Sigma_{xx}\rVert\leq\kappa\epsilon\;.$
(72)

####  E.3.5 Estimating $S_{x}^{1/2}$

Finally we can start to approach bounding our target quantity.

$\displaystyle\|\hat{S}_{x}-S_{x}\|$
$\displaystyle=\|(\hat{\Omega}_{x}-\Omega_{x})-(\hat{\Omega}_{xy}\hat{\Omega}_{y}^{-1}\hat{\Omega}_{yx}-\Omega_{xy}\Omega_{y}^{-1}\Omega_{yx})\|$
$\displaystyle\leq\|\hat{\Omega}_{x}-\Omega_{x}\|+\|(\hat{\Omega}_{xy}-\Omega_{xy})\hat{\Omega}_{y}^{-1}\hat{\Omega}_{yx}+\Omega_{xy}(\hat{\Omega}_{y}^{-1}-\Omega_{y}^{-1})\hat{\Omega}_{yx}+\Omega_{yx}\Omega_{y}^{-1}(\hat{\Omega}_{yx}-\Omega_{yx})\|$
$\displaystyle\leq\epsilon+\epsilon\cdot\frac{M}{1-M\epsilon}\cdot(\epsilon+\gamma_{1}m)+\gamma_{1}m\cdot\frac{M^{2}\epsilon}{1-M\epsilon}\cdot(\epsilon+\gamma_{1}m)+\gamma_{1}m\cdot\frac{M}{1-M\epsilon}\cdot\epsilon$
$\displaystyle=\epsilon\Bigl{\\{}1+\frac{\gamma_{1}Mm}{1-M\epsilon}(1+\gamma_{1}mM+1)\Bigr{\\}}+\epsilon^{2}\frac{M}{1-M\epsilon}(1+\gamma_{1}Mm)\;.$

Now note that $\epsilon\leq\frac{1}{4M}$ means that
$\frac{1}{1-M\epsilon}\leq\frac{4}{3}$ and so if we write
$\tau_{1}=\gamma_{1}Mm$ we can simplify this bound to

$\displaystyle\|\hat{S}_{x}-S_{x}\|$
$\displaystyle\leq\epsilon\cdot\Bigl{\\{}1+\frac{1}{1-M\epsilon}\left(2\tau_{1}+\tau_{1}^{2}+\epsilon
M(1+\tau_{1})\right)\Bigr{\\}}$
$\displaystyle\leq\epsilon\cdot\frac{4}{3}\Bigl{\\{}\frac{3}{4}+2\tau_{1}+\tau_{1}^{2}+\tfrac{1}{4}(1+\tau_{1})\Bigr{\\}}$
$\displaystyle\leq\epsilon\cdot\frac{4}{3}\left(\tau_{1}+\tfrac{9}{8}\right)^{2}=\vcentcolon
L\epsilon\;.$

where we defined
$L\vcentcolon=\frac{4}{3}\left(\tau_{1}+\tfrac{9}{8}\right)^{2}$. Next apply
Lemma 33, and note $1/\lambda_{\text{min}}(S_{x})=\lVert
S_{x}^{-1}\rVert=\lVert\Sigma_{xx}\rVert\geq M$ to get

$\|\hat{S}_{x}^{1/2}-S_{x}^{1/2}\|\leq\frac{\|\hat{S}_{x}-S_{x}\|}{\sqrt{\lambda_{\text{min}}(S_{x})}}\leq
M^{\nicefrac{{1}}{{2}}}L\epsilon\;.$ (73)

Therefore we can take $C_{2}(M,m,\rho_{1})=M^{\nicefrac{{1}}{{2}}}L$. We
similarly obtain $\|\hat{S}_{y}^{1/2}-S_{y}^{1/2}\|\leq
M^{\nicefrac{{1}}{{2}}}L\epsilon$.

####  E.3.6 Final target bound

Finally plug this in to the expression for $T$ to obtain:

$\displaystyle\|\hat{T}-T\|$
$\displaystyle=\|(\hat{S}_{x}^{1/2}-S_{x}^{1/2})\hat{\Sigma}_{xy}\hat{S}_{xy}^{1/2}+S_{x}^{1/2}(\hat{\Sigma}_{xy}-\Sigma_{xy})\hat{S}_{y}^{1/2}+S_{x}^{1/2}\Sigma_{xy}(\hat{S}_{y}^{1/2}-S_{y}^{1/2})\|$
$\displaystyle\leq
ML\epsilon\cdot(\rho_{1}M+\kappa\epsilon)\cdot(\sqrt{m}+\frac{L\epsilon}{2\sqrt{m}})+\sqrt{m}\cdot\kappa\epsilon\cdot(\sqrt{m}+\frac{L\epsilon}{2\sqrt{m}})+\sqrt{m}\cdot\rho_{1}M\cdot
ML\epsilon$
$\displaystyle=\epsilon\biggl{\\{}ML\rho_{1}M\sqrt{m}+m\kappa+\sqrt{m}\rho_{1}MML\biggl{\\}}$
$\displaystyle\qquad+\epsilon^{2}\biggl{\\{}ML(\rho_{1}M\frac{L}{2\sqrt{m}}+\sqrt{m}\kappa)+\sqrt{m}\kappa\frac{L}{2\sqrt{m}}\biggr{\\}}$
$\displaystyle\qquad+\epsilon^{3}\cdot ML\kappa\frac{L}{2\sqrt{m}}$
$\displaystyle\leq\epsilon\biggl{\\{}2\rho_{1}\sqrt{m}M^{2}L+m\kappa\biggl{\\}}$
$\displaystyle\qquad+\epsilon\biggl{\\{}\rho_{1}\frac{ML^{2}}{8\sqrt{m}}+\tfrac{1}{4}L\kappa\sqrt{m}+\tfrac{1}{2}L\cdot
2M)+\sqrt{m}\kappa\frac{L}{2\sqrt{m}}\biggr{\\}}$
$\displaystyle\qquad+\epsilon\cdot\frac{ML^{2}}{8\sqrt{m}}$
$\displaystyle\leq\epsilon\biggl{\\{}2\rho_{1}\sqrt{m}M^{2}L+8mM^{2}+\rho_{1}\frac{ML^{2}}{8\sqrt{m}}+2LM^{2}\sqrt{m}+LM+4M^{2}L\biggr{\\}}$
$\displaystyle=\vcentcolon C_{1}\epsilon\;.$

where we note $C_{1}$ does indeed only depend on $M,m,\text{ and }\rho_{1}$
and to leading order is $M^{4}m^{5/2}$. Some of the results of Gao et al.
[2016] rely on similar arguments; above we strive to make more explicit the
high-order polynomial dependence of the Lipschitz constant on condition
numbers $M$ and $m$. ∎

###  E.4 Variate space bounds

####  E.4.1 Davis-Kahan variant

The below is a special case of Theorem 4 from Yu et al. [2015]; we change to
our notation, only consider finding the top-$K$ subspace, and only consider
operator norm bounds.

######  Theorem 34.

Let $T,\hat{T}\in\mathbb{R}^{p\times q}$ have singular values
$\rho_{1}\geq\ldots\geq\rho_{\min(p,q)}$ and $\hat{\rho}_{1}\geq\ldots\geq$
$\hat{\rho}_{\min(p,q)}$ respectively. Fix $K<\operatorname{rank}(T)$ and
assume that $\rho_{K}^{2}-\rho_{K+1}^{2}>0$. Let
$V=\left(v_{1},v_{2},\ldots,v_{K}\right)\in\mathbb{R}^{q\times K}$ and
$\hat{V}=\left(\hat{v}_{1},\hat{v}_{2},\ldots,\hat{v}_{K}\right)\in\mathbb{R}^{q\times
K}$ have orthonormal columns satisfying $Tv_{j}=\rho_{j}u_{j}$ and
$\hat{T}\hat{v}_{j}=\hat{\rho}_{j}\hat{u}_{j}$ for $j=1,2,\ldots,K.$ Then

$\|\sin\Theta(\hat{V},V)\|_{\mathrm{F}}\leq\frac{2K^{1/2}\left(2\rho_{1}+\|\hat{T}-T\|\right)\|\hat{T}-T\|}{\rho_{K}^{2}-\rho_{K+1}^{2}}\;.$

To apply this result, we will take
$A=(\alpha_{1},\ldots,\alpha_{K})\in\mathbb{R}^{q\times K}$ and
$\hat{A}=(\hat{\alpha}_{1},\ldots,\hat{\alpha}_{K})\in\mathbb{R}^{q\times K}$
to be the top-$K$ left singular vectors of $T,\hat{T}$ respectively, as in
Section B.2.1. However, this $\sin\Theta$ loss between $\hat{A}$ and $A$ is
not particularly intuitive. We need a couple more steps to relate this to the
variate space loss that we used in the main text.

####  E.4.2 Converting to variate space losses

By applying Lemma 21 and noting that the true canonical variates are such that
$X^{\top}{U_{K}^{*}}$ has orthonormal components, we can rewrite the loss we
want to control as

$\displaystyle
L(\hat{U}_{K},U^{*}_{K})\vcentcolon=\sin^{2}\Theta(X^{\top}\hat{U}_{K},X^{\top}U_{K}^{*})=\inf_{W\in\mathbb{R}^{K\times
K}}\mathbb{E}\|X^{\top}\hat{U}_{K}W-X^{\top}U_{K}^{*}\|^{2}\;.$ (74)

We now compare this to the loss from Gao et al. [2016], which in our notation
is defined as

$L^{\text{Gao}}(U_{K},\hat{U}_{K})=\inf_{W\in\mathcal{O}(K)}\mathbb{E}\|X^{\top}\hat{U}_{K}W-X^{\top}U_{K}^{*}\|^{2}\;.$
(75)

We see that $L(\hat{U}_{K},U^{*}_{K})\leq
L^{\text{Gao}}(\hat{U}_{K},U^{*}_{K})$ because this second infimum is over
$\mathcal{O}(K)$, which is a subset of $\mathbb{R}^{K\times K}$. We therefore
bound this stronger loss $L^{\text{Gao}}$ and deduce the same bound must hold
for our loss $L$.

We now work on controlling the term on the right hand side of Equation 75.
From now we omit the subscript $K$ terms to declutter the notation. For
orthogonal $W$ we have

$\displaystyle\left(\mathbb{E}\|X^{\top}\hat{U}W-X^{\top}U^{*}\|^{2}\right)^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\left(\mathbb{E}\operatorname{tr}(W^{\top}\hat{U}^{\top}-U^{\top})XX^{\top}(\hat{U}W-U)\right)^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\left(\operatorname{tr}(\hat{U}W-U)^{\top}\Sigma_{xx}(\hat{U}W-U)\right)^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\|\Sigma_{xx}^{1/2}(\hat{U}W-U)\|_{F}$
$\displaystyle\leq\|\hat{\Sigma}_{xx}^{1/2}\hat{U}W-\Sigma_{xx}^{1/2}U\|_{F}+\|\hat{\Sigma}_{xx}^{1/2}\hat{U}W-\Sigma_{xx}^{1/2}\hat{U}W\|_{F}$
$\displaystyle=\|\hat{\Sigma}_{xx}^{1/2}\hat{U}W-\Sigma_{xx}^{1/2}U\|_{F}+\|(\hat{\Sigma}_{xx}^{1/2}-\Sigma_{xx}^{1/2})\hat{U}\|_{F}\;,$

where we first expand out using the trace-trick, then apply the triangle
inequality, and note that the orthogonal $W$ can be dropped from the final
Frobenius norm term. We can further bound this final term by

$\displaystyle\|(\hat{\Sigma}_{xx}^{1/2}-\Sigma_{xx}^{1/2})\hat{U}\|_{F}$
$\displaystyle=\|(I-\Sigma_{x}^{1/2}\hat{\Sigma}_{xx}^{-1/2})\hat{\Sigma}_{xx}^{1/2}\hat{U}\|_{F}$
$\displaystyle=\|(I-\Sigma_{x}^{1/2}\hat{\Sigma}_{xx}^{-1/2})\hat{A}\|_{F}$
$\displaystyle\leq
K^{\nicefrac{{1}}{{2}}}\,\|I-\Sigma_{x}^{1/2}\hat{\Sigma}_{xx}^{-1/2}\|_{\textrm{op}}$
$\displaystyle=K^{\nicefrac{{1}}{{2}}}\,\|\Sigma_{x}^{1/2}(\Sigma_{xx}^{-1/2}-\hat{\Sigma}_{xx}^{-1/2})\|_{\textrm{op}}$
$\displaystyle\leq
K^{\nicefrac{{1}}{{2}}}M^{1/2}\|\Sigma_{xx}^{-1/2}-\hat{\Sigma}_{xx}^{-1/2}\|_{\textrm{op}}\;,$

where for the first inequality we use that $\hat{A}$ has orthonormal columns
by construction.

Before we put everything together, we state the following standard lemma, e.g.
see Gao et al. [2015], Stewart and Sun [1990].

######  Lemma 35.

Suppose $A,\hat{A}\in\mathbb{R}^{p\times d}$, with $d<p$, have orthonormal
columns. Then

$\inf_{W\in\mathcal{O}(K)}\|A-\hat{A}W\|_{F}\leq\sqrt{2}\|\sin\Theta(\hat{A},A)\|_{F}\;.$

Now we can combining the previous bounds and apply Lemma 35 to give

$\displaystyle L^{\text{Gao}}(U,\hat{U})^{1/2}$
$\displaystyle\leq\inf_{W\in\mathcal{O}(K)}\left(\mathbb{E}\|X^{\top}\hat{U}W-X^{\top}U^{*}\|^{2}\right)^{\nicefrac{{1}}{{2}}}$
$\displaystyle\leq\inf_{W\in\mathcal{O}(K)}\Bigl{\\{}\|\hat{\Sigma}_{xx}^{1/2}\hat{U}W-\Sigma_{xx}^{1/2}U\|_{F}+K^{\nicefrac{{1}}{{2}}}M^{1/2}\|\Sigma_{xx}^{-1/2}-\hat{\Sigma}_{xx}^{-1/2}\|\Bigr{\\}}$
$\displaystyle\leq\sqrt{2}\|\sin\Theta(\hat{A},A)\|_{F}+K^{\nicefrac{{1}}{{2}}}M^{1/2}\|\Sigma_{xx}^{-1/2}-\hat{\Sigma}_{xx}^{-1/2}\|_{\textrm{op}}\;.$
(76)

This is now in a format that we can handle with our results from Section E.3.
Recall also that the same bound must apply to $L$, since $L^{\text{Gao}}$
bounds $L$.

####  E.4.3 Combining results

######  Proposition 36.

Assume Assumption A$(M,m)$. Then there exists a constant $C_{3}\in(0,\infty)$
only depending on $M,m,\rho_{1}$ such that whenever

$\epsilon\vcentcolon=\|\Omega-\hat{\Omega}\|\leq\min\left(\frac{1}{4M},\frac{2\rho_{1}}{C_{1}}\right)$
(77)

we have

$\lVert\sin\Theta(X^{\top}\hat{U}_{K},X^{\top}U_{K}^{*})\rVert_{F}\leq
L^{\text{Gao}}(\hat{U}_{K},U_{K}^{*})^{\nicefrac{{1}}{{2}}}\leq
C_{3}\frac{K^{\nicefrac{{1}}{{2}}}}{\rho_{K}^{2}-\rho_{K+1}^{2}}\epsilon\;.$
(78)

###### Proof.

Suppose

$\epsilon\vcentcolon=\|\Omega-\hat{\Omega}\|\leq\min\left(\frac{1}{4M},\frac{2\rho_{1}}{C_{1}}\right)$

as in the statement of this result. Then by Proposition 29, there exist
$C_{1},C_{2}$ such that

$\lVert\hat{T}-T\rVert_{\textrm{op}}\leq
C_{1}\epsilon,\quad\text{and}\quad\lVert\hat{\Sigma}_{xx}-\Sigma_{xx}\rVert_{\textrm{op}}\leq
C_{2}\epsilon\;.$

Applying Theorem 34 and using that $C_{1}\epsilon\leq 2\rho_{1}$ then yields

$\displaystyle\|\sin\Theta(\hat{A},A)\|_{\mathrm{F}}\leq\frac{2K^{1/2}\left(2\rho_{1}+\|\hat{T}-T\|\right)\|\hat{T}-T\|}{\rho_{K}^{2}-\rho_{K+1}^{2}}\leq\frac{8K^{\nicefrac{{1}}{{2}}}\rho_{1}}{\rho_{K}^{2}-\rho_{K+1}^{2}}C_{1}\epsilon\;.$

Plugging this into Equation 76 and using that
$\frac{1}{\rho_{K}^{2}-\rho_{K+1}^{2}}\geq 1$ gives

$\displaystyle L^{\text{Gao}}(U,\hat{U})^{1/2}$
$\displaystyle\leq\sqrt{2}\|\sin\Theta(\hat{A},A)\|_{F}+K^{\nicefrac{{1}}{{2}}}M^{\nicefrac{{1}}{{2}}}\|\Sigma_{xx}^{-1/2}-\hat{\Sigma}_{xx}^{-1/2}\|$
$\displaystyle\leq\frac{8\sqrt{2}K^{\nicefrac{{1}}{{2}}}\rho_{1}}{\rho_{K}^{2}-\rho_{K+1}^{2}}C_{1}\epsilon+K^{\nicefrac{{1}}{{2}}}M^{\nicefrac{{1}}{{2}}}C_{2}\epsilon$
$\displaystyle\leq(8\sqrt{2}\rho_{1}C_{1}+M^{\nicefrac{{1}}{{2}}}C_{2})\frac{K^{\nicefrac{{1}}{{2}}}}{\rho_{K}^{2}-\rho_{K+1}^{2}}\epsilon\;,$

and so the desired bound holds by taking
$C_{3}=8\sqrt{2}\rho_{1}C_{1}+M^{\nicefrac{{1}}{{2}}}C_{2}$. ∎

Note that applying instead to
$B=(\beta_{1},\ldots,\beta_{K}),\hat{B}=(\hat{\beta}_{1},\ldots,\hat{\beta}_{K})$
gives us an identical bound for the $Y$-direction estimates. These two bounds
are identical due to the symmetry of Theorem 34; we expect a stronger
asymmetric statement could be attained from the results of Cai and Zhang
[2020], but leave further exploration for future work.

###  E.5 Combining with existing bounds for graphical lasso

####  E.5.1 Structural assumptions from Ravikumar et al. [2011]

In the main text, we summarised the structural assumptions from Ravikumar et
al. [2011] in Assumption B, but did not define the terms formally. We now give
formal definitions for each of the terms. We first restate the set of
assumptions for reference.

See B

$\ell_{\infty}$-operator-norm bound:

$\displaystyle\kappa_{\Sigma}\vcentcolon=\|\Sigma\|_{\mathrm{op},\ell_{\infty}\to\ell_{\infty}}=\max_{i=1,\dots,\bar{p}}\sum_{j=1}^{\bar{p}}|\Sigma_{ij}|\;.$

Strict convexity: Let $S(\Omega)$ be the set of nonzero index pairs for
$\Omega$, augmented with the index pairs of the diagonal entries. Then we
define

$\displaystyle\Gamma_{SS}$
$\displaystyle\vcentcolon=[{\Omega}^{-1}\otimes{\Omega}^{-1}]_{SS}\in\mathbb{R}^{(s+\bar{p})\times(s+\bar{p})}\;,$
$\displaystyle\kappa_{\Gamma}$
$\displaystyle\vcentcolon=\|\left(\Gamma_{SS}\right)^{-1}\|_{\mathrm{op}}\;.$

Mutual incoherence: There exists some $\alpha\in(0,1]$ such that

$\max_{e\in
S^{c}}\|\Gamma_{eS}\left(\Gamma_{SS}\right)^{-1}\|_{1}\leq(1-\alpha)\;.$

This may also be referred to as an irrepresentability condition.

####  E.5.2 Adapting results from Ravikumar et al. [2011]

The main result in Ravikumar et al. [2011] is deliberately very general: this
covers a broad class of distributions for $Z$ with a variety of different tail
conditions, and bounds the elementwise sup-norm
$\lVert\hat{\Omega}-\Omega\rVert_{\infty}$. For our purposes, it is sufficient
to consider exponential-type tails (original random variables are sub-
Gaussian), and to obtain bounds on the operator norm
$\lVert\hat{\Omega}-\Omega\rVert_{\text{op}}$. The authors provide a number of
helpful corollaries to help the reader obtain bounds in similar, specific
situations of interest. In particular, Corollary 1 translates their general
result to the situation of exponential tails. Then Corollary 3 translates
their element-wise-sup-norm bound to an operator norm bound.

By combining the analysis from these two special cases, we can read off the
most special case of interest to us.

######  Proposition 37 (Synthesis of results from Ravikumar et al. [2011]).

Let Assumption B$\,(\kappa_{\Gamma},\kappa_{\Sigma},\alpha)$ hold. Suppose the
variables $Z_{i}/\sqrt{\Sigma_{ii}}$ are sub-Gaussian with common parameter
$\sigma$ and the samples are drawn independently. Let

$\displaystyle\gamma_{*}=8\sqrt{2}(1+4\sigma^{2})\max_{i}(\Sigma_{ii})\;.$

Let $\tau$ be a user-specified parameter, which determines the $\ell_{1}$
penalty parameter via

$\displaystyle\lambda_{N}=\frac{8\gamma_{*}}{\alpha}\sqrt{\frac{\tau\log\bar{p}+\log
4}{N}}\;.$

Suppose that the sample size $n$ satisfies the bound

$\displaystyle
N>\left\\{6\gamma_{*}\max\left\\{\kappa_{\Sigma}\kappa_{\Gamma},\kappa_{\Sigma}^{3}\kappa_{\Gamma}^{2}\right\\}\right\\}d^{2}\left(1+\frac{8}{\alpha}\right)(\tau\log\bar{p}+\log
4)\;.$

Then with probability greater than $1-\bar{p}^{2-\tau}$, the estimate
$\hat{\Omega}$ using penalty parameter $\lambda_{N}$ satisfies the bound

$\displaystyle\lVert\hat{\Omega}-\Omega\rVert_{\text{op}}\leq\left\\{2\gamma_{*}\left(1+\frac{8}{\alpha}\right)\kappa_{\Gamma}\right\\}\sqrt{\frac{\tau\log\bar{p}+\log
4}{N}}\;\min\\{\sqrt{s+\bar{p}},d\\}\;.$

######  Proposition 38 (Simplified version of Proposition 37).

Let Assumption B$\,(\kappa_{\Gamma},\kappa_{\Sigma},\alpha)$ hold. Suppose the
variables $Z_{i}/\sqrt{\Sigma_{ii}}$ are sub-Gaussian with common parameter
$\sigma$ and the samples are drawn independently. Let $\tau\geq 2$ be a user-
specified parameter, and $\lambda_{N}$ be an appropriate penalty parameter
determined by $\tau$. Suppose $p\geq 2$. Then there exist constants
$C_{4},C_{5}$ only depending on
$\left(\kappa_{\Gamma},\kappa_{\Sigma},\alpha,\sigma,\max_{i}(\Sigma_{ii})\right)$
such that when the sample size $N$ satisfies

$\displaystyle N>C_{4}\,d^{2}\tau\log\bar{p}$

we have the bound

$\displaystyle\lVert\hat{\Omega}-\Omega\rVert_{\text{op}}\leq
C_{5}\sqrt{\frac{\tau\min\\{{s+\bar{p}},d^{2}\\}\log\bar{p}}{N}}\;.$

####  E.5.3 Pulling everything together

By combining Proposition 38 with Proposition 36 we arrive at our result from
the main text. See 3

##  Appendix F Background continued

CCA has been a central concept in multi-view learning and multivariate
analysis since its proposal in Hotelling [1936]. There has been a vast amount
of work around this concept, and textbook treatments of the classical
methodology are available in Anderson [2003], for example.

Section F.1 gives a more thorough literature review on regularised CCA, which
is important to place our approach in the context of existing methods. Solari
et al. [2019] and Uurtio et al. [2017] contain solid literature reviews on
CCA, which touch on extensions and related methods. Section F.2 complements
these existing surveys by summarising selected recent developments; we
emphasise work that has been particularly helpful for developing our
intuitions for understanding CCA.

###  F.1 High-dimensional CCA: further background

####  F.1.1 Unregularised CCA: assorted results

Ma and Li [2020] provide non-asymptotic bounds for the (classical) sample CCA
estimator in the $n>p+q$ regime. Their loss functions use $\sin\Theta$
distances in variate space, and are well motivated by some insightful
discussion. We took inspiration from much of their presentation.

####  F.1.2 sPLS: selected extensions

A number of extensions and variants to the sPLS estimator we define have been
proposed; these include considering more general penalties such as fused or
group lasso penalties [Witten et al., 2009], extensions to data with more than
two views, and a computational short-cut [Solari et al., 2019]. However,
fundamentally, all these variants differ from the standard CCA problem by an
implicit assumption of identity-block covariance, as discussed in Section 3.

####  F.1.3 Sparse CCA: methodological suggestions

We now summarise a variety of recent works proposing sparse CCA methods
building upon or distinct to the sCCA [Suo et al., 2017] method we considered
in the main text. Though some of these works contain theoretical guarantees,
the main contributions appear methodological.

  * •

Earlier work: The first work we are aware of not to require the identity-
block-covariance assumption made by sPLS was that of Hardoon and Shawe-Taylor
[2011]; this uses a so-called ML-primal-dual framework, but requires a
different restrictive assumption: namely positivity for one of the canonical
directions.

  * •

Implementation of Suo et al. [2017]: For each step of the ACS procedure, Suo
et al. [2017] originally suggested to solve the resulting convex optimsation
problem with a certain linearised alternating direction method of multipliers
(LADMM) [Parikh and Boyd, 2014]. The more recent work Xiu et al. [2020]
proposes an algorithm that still uses ACS but replaces LADMM with a semi-
smooth Newton method; they claim this gives significantly increased
performance.

  * •

Iterative Penalised Least Squares: Another natural approach makes use of the
alternating-projection formulation of CCA [Anderson, 2003]; this was suggested
as an exercise in Hastie et al. [2015] and analysed in Mai and Zhang [2019],
which we shall refer to as Iterative Penalised Least Squares (IPLS). Though
relatively fast and easy to implement, this requires deflation to encourage
orthogonality of successive estimated directions; we found this deflation
caused problems with later directions when we tried to apply the method in
high-dimensional settings.

  * •

Linear Programming: Safo et al. [2018b] and its sister work Safo et al.
[2018a], both use a ridge-regularised covariance estimate666They alternatively
suggest a restrictive identity block-covariance approximation but this then
resembles a PLS algorithm, so we focus on the ridge case. to obtain initial
canonical direction estimates, then find sparse vectors which are ‘near’ these
initial estimates. Formally, if the ridge estimates are
$(\tilde{u}_{k},\tilde{v}_{k},\tilde{\rho}_{k})$ then we define

$\displaystyle\hat{u}_{k}^{\text{DS}}\in\operatorname*{arg\,min}_{u\in\mathbb{R}^{p}}\mathcal{P}(u,\tau_{x})$
$\displaystyle\quad\text{subject
to}\quad\|\mathbf{C}_{xy}\tilde{v}_{k}-\tilde{\rho}_{k}\mathbf{C}_{xx}u\|_{\infty}\leq\tau_{x}\;,$
$\displaystyle\hat{v}_{k}^{\text{DS}}\in\operatorname*{arg\,min}_{v\in\mathbb{R}^{p}}\mathcal{P}(v,\tau_{y})$
$\displaystyle\quad\text{subject
to}\quad\|\mathbf{C}_{yx}\tilde{u}_{k}-\tilde{\rho}_{k}\mathbf{C}_{yy}v\|_{\infty}\leq\tau_{y}\;,$

where $\mathcal{P}$ is some penalty function which may incorporate prior
biological knowledge; for example $\mathcal{P}(u,\tau_{x})=\|u\|_{1}$. We can
think of this as first estimating the underlying variates with a ridge penalty
and then finding sparse directions which give similar loadings.

  * •

Sparse GEP:  Tan et al. [2018] propose a Truncated Rayleigh Flow Method
(Rifle) for the more general problem of solving a sample Generalised
Eigenvalue Problem whose eigenvectors are assumed to be sparse. Their convex-
optimisation framework closely resembles that of Gao et al. [2016]. They
provide interesting theory which use matrix perturbation results which exploit
the geometry of GEPs from Stewart and Sun [1990].

####  F.1.4 Sparse CCA: theoretical contributions

We now summarise some recent works whose main contributions appear
theoretical, though they do all also provide methodological suggestions. They
each contain very impressive mathematical analysis, which require assumptions
of sparse directions. Following the discussion in the main text, we argue that
it is important in practice to be cautious about these assumptions and
contrast sCCA estimators with alternatives.

  * •

Minimax Rates:  Gao et al. [2015] and the follow-up Gao et al. [2016] contain
seminal theoretical work establishing minimax rates for sparse CCA under row-
wise sparsity conditions. The former paper proposes a combinatorial estimator
that is minimax but computationally intractable. The second paper proposes an
estimator (COLAR), which uses a convex optimisation in matrix space to obtain
a sparse-and-reduced-rank approximation to the matrix $U\Lambda V^{\top}$.
Unfortunately, their analysis in this case requires extra assumptions (of
knowing the true number of canonical directions). Though this is the most
principled sparse CCA method proposed, it is computationally expensive
(complexity $\mathcal{O}(p+q)^{3}$), and challenging to implement (2 stage
ADMM). We did not spend time implementing the method, and would welcome a
contribution to the open source python community and to cca-zoo in particular.

  * •

Support Recovery:  Laha and Mukherjee [2021] contains some sophisticated
theoretical results on information theoretic and computational hardness of
support recovery. They propose a Coordinate Thresholding algorithm for
estimation of canonical directions, and a further algorithm, RecoverSupp, for
support recovery (again via thresholding). Their method requires access to
‘good’ precision matrix estimates, among other complex assumptions, including
sparsity of the canonical directions which is essential to the premise of
support recovery.

  * •

Inference:  Laha et al. [2022] provide more sophisticated theory giving
asymptotically valid confidence intervals for sCCA. They use preliminary
estimators and then a de-biasing step which uses a smooth objective function
to characterise the first pair of canonical directions. The theory connects to
the more general area of de-biased/de-sparsified inference in high-dimensional
models. Their method only concerns the first pair of canonical directions, and
again assumes access to ‘good’ block precision estimates. The assumption of
sparse canonical directions is of course essential to the validity of
inference, and the algorithm is complicated to implement and use in practice.

###  F.2 Select extensions and alternative formulations of CCA

We briefly highlight three main avenues to reformulate or extend CCA; these
have all led to interesting empirical results and may help develop intuition
for the various different ways to look at CCA.

Firstly, CCA can be viewed through the lens of ‘multi-view’ learning; it is
therefore natural to consider generalising CCA to handle 3 or more datasets. A
large number of approaches have been proposed for multiview-CCA [Tenenhaus and
Tenenhaus, 2011, Rodosthenous et al., 2020, Fu et al., 2017] but there is no
accepted ‘canonical’ formulation and these techniques are much less widely
used than standard two-view CCA.

The second main direction is to consider learning highly correlated
representations that are nonlinear transformations of the input variables
(rather than linear transformations). Historically, this was first achieved by
kernel CCA methods, which maximise correlation between the transformed
variables over a class of nonlinear transformations in a certain Reproducing
Kernel Hilbert Space (RKHS) [Bach and Jordan, 2002]. More recently, various
notions of DeepCCA [Andrew et al., 2013, Chapman et al., 2023] instead
consider neural network transformations; these DeepCCA algorithms have had
empirical success, now enjoy fast implementations [Chapman et al., 2023], and
are closely related to recent advances in self-supervised learning (SSL)
[Zbontar et al., 2021].

Thirdly, there have been a number of probabilistic, or Bayesian approaches or
interpretations of CCA; the excellent theoretical work of Bach and Jordan
[2005] laid the foundations for later empirical work, such as Klami et al.
[2013]. These provide a natural and concrete link between CCA and the field of
latent variable models [Roweis and Ghahramani, 1999]; this is relevant to our
narrative of ‘loadings not weights’, as we describe in the following
subsection.

###  F.3 Probabilistic interpretations of CCA

####  F.3.1 Probabilistic CCA [Bach and Jordan, 2005]

We first translate the main probabilistic CCA result from Bach and Jordan
[2005] into our notation and then give further interpretation.

######  Proposition 39.

Consider the model where $X$ and $Y$ are conditionally independent given a
latent variable $Z$, with

$\displaystyle Z$ $\displaystyle\sim\mathcal{N}(0,I_{d}),$
$\displaystyle\min(p,q)\geq d\geq 1$ $\displaystyle X|Z$
$\displaystyle\sim\mathcal{N}(W_{1}Z+\mu_{1},\Psi_{1}),$ $\displaystyle
W_{1}\in\mathbb{R}^{p\times d},\Psi_{1}\succcurlyeq 0$ $\displaystyle Y|Z$
$\displaystyle\sim\mathcal{N}(W_{2}Z+\mu_{2},\Psi_{2}),$ $\displaystyle
W_{2}\in\mathbb{R}^{p\times d},\Psi_{2}\succcurlyeq 0\;.$

Then the MLEs of the parameters
$W_{1},W_{2},\Psi_{1},\Psi_{2},\mu_{1},\mu_{2}$ are

$\displaystyle\hat{W}_{1}$
$\displaystyle=\mathbf{C}_{xx}\hat{U}_{d}M_{1}\qquad$
$\displaystyle\hat{\Psi}_{1}$
$\displaystyle=\mathbf{C}_{xx}-\hat{W}_{1}\hat{W}_{1}^{\top}\qquad$
$\displaystyle\hat{\mu}_{1}$ $\displaystyle=\bar{\mathbf{X}}$
$\displaystyle\hat{W}_{2}$ $\displaystyle=\mathbf{C}_{yy}\hat{V}_{d}M_{2}$
$\displaystyle\hat{\Psi}_{2}$
$\displaystyle=\mathbf{C}_{yy}-\hat{W}_{2}\hat{W}_{2}^{\top}$
$\displaystyle\hat{\mu}_{2}$ $\displaystyle=\bar{\mathbf{Y}}$

where $\hat{U}_{d}\in\mathbb{R}^{p\times d},\hat{V}_{d}\in\mathbb{R}^{q\times
d}$ are matrices whose columns are successive sample canonical directions and
$M_{1},M_{2}\in\mathbb{R}^{d\times d}$ are arbitrary matrices with
$M_{1}M_{2}^{\top}=\hat{R},\>\|M_{1}\|_{\textrm{op}}\leq
1,\|M_{2}\|_{\textrm{op}}\leq 1$, where $\hat{R}$ is the diagonal matrix of
sample canonical correlations.

The key observation here is that the MLEs $\hat{W}_{1},\hat{W}_{2}$ are
matrices of loadings and not weights; this is exactly analogous to the results
in Section B.2.2: weights are required to obtain canonical variates from input
variables, but the loadings are the right basis with which to reconstruct the
input variables from the canonical variates, treating these canonical variates
as latent variables. This probabilistic result can be interpreted as a
rigorous statement of ‘optimal’ reconstruction, with respect to a natural
linear Gaussian latent variable model.

####  F.3.2 Mutual Information

Mutual information provides another rigorous statement of ‘optimality’ for CCA
in the Gaussian setting: CCA defines low dimensional linear projections of
maximal mutual information. Indeed there is a well known analytic formula
expressing mutual information between two Gaussian random variables as a
function of their canonical correlations [Borga, 1998]. We give a slightly
different proof here for completeness, and to aid intuition. This motivates
the alternative choice of metric for combining canonical correlations we
proposed in Section 5.1.1.

######  Proposition 40.

For $X,Y$ jointly Gaussian with canonical correlations $(\rho_{k})_{k=1}^{K}$,
the mutual information

$I(X;Y)=\frac{1}{2}\log\left(\frac{1}{\Pi_{k=1}^{K}(1-\rho_{k}^{2})}\right)=\frac{1}{2}\sum_{k=1}^{K}\shortminus\log(1-\rho_{k}^{2})\;.$
(79)

###### Proof.

We state three well known properties of (continuous) mutual information [Cover
and Thomas, 2012]:

  * •

Mutual information can be expressed in terms of differential entropies,

$\displaystyle I(X;Y)=h(X)+h(Y)-h(X,Y)\;.$ (80)

  * •

Mutual information is invariant to invertible linear transformations777In fact
also invariant to more general invertible transformations, but this is harder
to prove., to see this note that for any random vector $X$ and deterministic
matrix $A$

$\displaystyle h(AX)=h(X)+\log|A|\;,$

where $|\cdot|$ denotes the determinant, and so the $|A|$ terms cancel in
Equation 80

  * •

The differential entropy of a Gaussian random variable with covariance matrix
$\Sigma$ is

$\displaystyle\tfrac{1}{2}\log\left\\{(2\pi e)^{n}|\Sigma|\right\\}\;.$

To conclude, let $U,V$ be full bases of canonical directions for $(X,Y)$, and
let $R=U^{\top}\Sigma_{xy}V=\text{diag}(\rho_{k})_{k=1}^{K}$. Then since $U,V$
are invertible we see

$\displaystyle I(X;Y)$ $\displaystyle=I(U^{\top}X;V^{\top}Y)$
$\displaystyle=\tfrac{1}{2}\log\left\\{(2\pi
e)^{p}|I_{p}|\right\\}+\tfrac{1}{2}\log\left\\{(2\pi
e)^{q}|I_{q}|\right\\}-\tfrac{1}{2}\log\left\\{(2\pi
e)^{(p+q)}\left|\begin{pmatrix}I_{p}&R\\\
R^{\top}&I_{q}\end{pmatrix}\right|\right\\}$
$\displaystyle=-\tfrac{1}{2}\log\left\\{\left|\begin{pmatrix}I_{p}&R\\\
R^{\top}&I_{q}\end{pmatrix}\right|\right\\}$
$\displaystyle=-\tfrac{1}{2}\log\prod_{k=1}^{K}(1-\rho_{k}^{2})\;,$

as claimed. ∎

This proof, together with CCA-interlacing Lemma 25, suggest that in the joint
Gaussian setting, CCA can be characterised as finding subspaces of maximal
mutual information.

Note that we now have two straightforward ways to compute mutual information
between Gaussian distributions with known covariance. First we can compute the
canonical correlations and apply the result Proposition 40. Alternatively, we
can use the intermediate result that for that joint Gaussian $X,Y$ with
covariance $\Sigma$, the mutual information is

$I(X;Y)=I^{\text{Gauss}}(\Sigma)\vcentcolon=\frac{1}{2}\log\left(\frac{|\Sigma_{xx}||\Sigma_{yy}|}{|\Sigma|}\right)\;.$
(81)

##  Appendix G Toy synthetic data experiments

In this section, we present two further synthetic data experiments. These
consider two toy data generating mechanisms. These experiments were produced
for an earlier iteration of this manuscript, and use different error metrics
to those from the main text. These are defined by

$\displaystyle\texttt{loadu1}=\frac{\lVert\Sigma_{xx}\hat{u}_{1}-\Sigma_{xx}u_{1}\rVert}{\lVert\Sigma_{xx}u_{1}\rVert},\quad\texttt{weightu1}=\frac{\lVert\hat{u}_{1}-u_{1}\rVert}{\lVert
u_{1}\rVert},\quad\texttt{rho1}=\rho^{\text{oracle}}(\hat{u}_{1},\hat{v}_{1})\;.$

###  G.1 Data generating mechanisms

####  G.1.1 Canonical Pair model

First we consider the canonical pair model, motivated by the decomposition of
Section B.2.1, and first proposed in Chen et al. [2013]; this is a
construction for non-trivial joint covariance matrices with sparse
corresponding canonical directions; as opposed to being physically
motivated888Interestingly, one of their key observations was that high-
dimensional sparse precision matrices can be well estimated, but they did not
think to estimate the (full) joint precision, only the within-view blocks.
Their particular choices of within-view covariance blocks leads to a certain
sparsity in the joint precision matrix, which helps explain why our gCCA
approach performs well.. The construction has since been used in Gao et al.
[2016], Suo et al. [2017], Safo et al. [2018a, b], Mai and Zhang [2019], Xiu
et al. [2020], Rodosthenous et al. [2020], Laha and Mukherjee [2021].

We generate multivariate Gaussian $\mathbf{X}\in\mathbb{R}^{N\times p}$ and
$\mathbf{Y}\in\mathbb{R}^{N\times q}$ from

$\left(\begin{array}[]{l}x\\\
y\end{array}\right)\sim\mathrm{N}\left(\left(\begin{array}[]{l}0\\\
0\end{array}\right),\left(\begin{array}[]{cc}\Sigma_{xx}&\Sigma_{xy}\\\
\Sigma_{yx}&\Sigma_{yy}\end{array}\right)\right)\;,$ (82)

where we vary the form of $\Sigma_{xx},\Sigma_{yy}$ but always generate

$\Sigma_{xy}\vcentcolon=\Sigma_{xx}\left(\sum_{k=1}^{K}\rho_{k}u_{k}v_{k}^{\top}\right)\Sigma_{yy}$
(83)

with $\rho_{k}\in(0,1)$ and $u_{k}\in\mathbb{R}^{m},v_{k}\in\mathbb{R}^{p}$
such that
$u_{k}^{\top}\Sigma_{xx}u_{j}=v_{k}^{\top}\Sigma_{yy}v_{j}=\delta_{jk}$. From
the arguments of Section B.2.1, we see that $\rho_{k},u_{k},v_{k}$ are
precisely the corresponding canonical correlations and directions.

For this section, we consider constructing the covariances
$\Sigma_{xx},\Sigma_{yy}$ via
$\Sigma_{xx}=(W_{m})^{-1},\Sigma_{yy}=(W_{p})^{-1}$ where
$W_{m}\in\mathbb{R}^{m\times m}$ is defined by

$\left(W_{m}\right)_{ij}=\mathbbm{1}_{\\{i=j\\}}+0.5\times\mathbbm{1}_{\\{|i-j|=1\\}}+0.4\times\mathbbm{1}_{\\{|i-j|=2\\}}\;.$

This is referred to as the ‘sparse precision’ model in Suo et al. [2017].
There is still plenty of flexibility in the choice of $\rho_{k},u_{k},v_{k}$.
In this case, we took a single canonical pair model with $\rho_{1}=0.9$ (and
$\rho_{k}=0$ for $k>1$). We constructed $u_{1},v_{1}$ to only have non-zero
entries in the first 5 components, with each such entry generated
independently from a uniform distribution on $[-1,1]$ (a single pair
$u_{1},v_{1}$ were generated before any samples drawn, and used in all the
runs). We will therefore denote this covariance structure as ‘suo-sp-rand’.
The dimensions were taken as $p=q=150$, and we varied the number of samples
$n$.

####  G.1.2 Graphical Model

We also consider covariance structures from a sparse graphical model. We use
the ‘power law’ construction for the covariance matrix corresponding to a
sparse graphical model as implemented in the gglasso package [Schaipp et al.,
2021]. Here the precision matrix is constructed using such that the
interaction graph has nodes which satisfy a power law distribution, with some
decay constant $\gamma$. Here we took $\gamma=3$. We denote this covariance
structure as ‘powerlaw’. Again we worked in dimensions $p=q=150$, and we
varied the number of samples $n$.

###  G.2 Results

Figure 9: Evaluation of the four different algorithms for suo-sp-rand model,
relative to the known covariance structure. Dimensions $p=q=150$ are fixed and
we vary number of samples $n$. In each case we took the average over 32 random
seeds. For the solid lines we use estimates with tuning parameters selected
via CV on $\rho_{1}$ (r2s1-cv), dotted lines select tuning parameters using
$\rho^{\text{oracle}}(\hat{u}_{1},\hat{v}_{1})$, i.e. the ‘oracle’ r2s1.
Figure 10: Analogous plot to Figure 9 but now for powerlaw model; otherwise
identical setup.

####  G.2.1 Competitive performance in single canonical pair models

We summarise some observations below:

  * •

Loading recovery is easier than weight recovery; and indeed we can recover
high correlations even with poor weight estimates.

  * •

Even in this most benign case of a single strong canonical correlation, high
sparsity, and moderate dimension, it requires a few hundred samples to get
good estimates of the loadings, let alone the weights.

  * •

sPLS performs well for small $n$ but fails to converge to the true values as
$n\to\infty$ (i.e. is not consistent). This is not surprising because sPLS
should converge to PLS decomposition rather than CCA in the limit of small
regularisation, so there is no reason to expect consistency in general.

  * •

gCCA converges only somewhat slower than sCCA, and much faster than rCCA .

In other similar scenarios, we generally saw that our gCCA method performed
nearly as well as existing alternatives.

####  G.2.2 Good performance in graph-structured models

Figure 10 is the analogue of Figure 9 but now with the powerlaw model. We note
some further observations:

  * •

Generally this problem appears easier than the suo-sp-rand case.

  * •

Again loadings are easier to estimate than weights.

  * •

Now gCCA is clearly the best performing model.

  * •

Again sPLS does not appear to converge to the CCA subspace.

##  Appendix H Parametric bootstrap: further synthetic experiments

###  H.1 Alternative regularisation for oracle covariance matrix

Algorithm 2 sCCA-and-ridge-regularised joint covariance matrix

procedure sCCA-regularised-
covariance($\mathbf{X},\mathbf{Y};\lambda,\alpha,K$)

$\Sigma_{xx},\Sigma_{yy}\leftarrow\mathbf{C}_{xx}+\alpha
I_{p},\mathbf{C}_{yy}+\alpha I_{q}$

$\hat{U}_{K},\hat{V}_{K}\leftarrow sCCA{}_{K}(\mathbf{X},\mathbf{Y};\lambda)$
$\triangleright$ Columns are the successive $\hat{u}_{k},\hat{v}_{k}$ from
Equation 12

$\hat{D}_{K}=\operatorname{diag}\left(\widehat{\operatorname{Corr}}(\mathbf{X}\hat{u}_{k},\mathbf{Y}\hat{v}_{k})_{k=1}^{K}\right)$

$\hat{\Sigma}_{xy}\leftarrow\hat{\Sigma}_{xx}\hat{U}_{K}\hat{D}_{K}\hat{V}_{K}^{\top}\hat{\Sigma}_{yy}$
$\triangleright$ C.f. Equation 4

return
$\hat{\Sigma}\vcentcolon=\left(\begin{array}[]{cc}\hat{\Sigma}_{xx}&\hat{\Sigma}_{xy}\\\
\hat{\Sigma}_{xy}^{\top}&\hat{\Sigma}_{yy}\end{array}\right)$

end procedure

In Section 6 we considered a multi-variate normal model where the true
(oracle) joint covariance matrix was obtained by regularising the sample
covariance of the Microbiome dataset with the graphical lasso. Here, we
construct this oracle joint covariance matrix using sCCA to encourage
structure with sparse canonical directions.

In more detail, we construct this matrix using Algorithm 2. First we
regularise the within-view covariance matrices using a ridge penalty, choosing
small but arbitrary parameter of $\alpha=0.01$. Then we obtain $K=10$
successive canonical direction estimates using sCCA and the r2s5-cv-optimal
penalty of $\lambda=0.013$. Finally we construct the between view covariance
matrices using Equation 4.

Figure 11 is the analogous plot to Figure 1 but for this sCCA-regularised
joint covariance matrix. The primary observation is that the figure looks
extremely similar to Figure 1; now sCCA appears to perform slightly better
than gCCA in terms of correlation captured (as one would expect); otherwise
all the observations from Section 6.1 still hold.

Figure 12 is the analogous plot to Figure 2 but for this sCCA-regularised
joint covariance matrix. Again, the primary observation is that the figure
looks extremely similar to Figure 2; again sCCA appears to perform slightly
better than gCCA in terms of correlation captured; otherwise all the
observations from Section 6.2 still hold.

Figure 11: Oracle metrics: analogous plot to Figure 1 on a parametric-
bootstrapped Microbiome dataset, where now the ‘oracle’ joint covariance
matrix is constructed using sCCA. See Tables 1 and 2 for a glossary of the
legends.

Figure 12: Sample splitting metrics (with comparable oracle quantities dotted
and overlaid): analogous plot to Figure 2 on a parametric-bootstrapped
Microbiome dataset, where now the ‘oracle’ joint covariance matrix is
constructed using sCCA. See Tables 2 and 1 for a glossary of the legends.

##  Appendix I Real data supporting plots

###  I.1 Microbiome

We now provide additional figures to support the claims from Section 8.3. We
summarise the main observations below, with figures in the following
subsections.

  * •

Trajectory comparison distances larger in weights space: Figure 17 has far
higher entries than the corresponding variate version Figure 5.

  * •

Overlap matrices are versatile: the three sets of overlap matrices in Figures
14, 15 and 16 illustrate the following points:

    * –

gCCA is stable across folds: In particular, the squared overlap matrices in
Figure 14 show how estimates for the r2s5-cv-optimal gCCA estimator vary
across 3 training folds, see caption for full details. The main observation
from Figure 14 is that the Glasso estimates are reassuringly stable across the
folds: the squared overlap matrices all appear close to identity, especially
for the first 3 components; moreover, the row and column sums of the first few
components are close to one. This is perhaps surprising giving the lack of
separation of the correlation signals from Figure 4. Note that we do also
observe some mixing up of the 4/5 components in the middle plot.

    * –

visualising registration: Figure 15 uses the same estimates as Figure 6 but
now the sCCA and rCCA solutions have been registered to the gCCA solution.
Observe all the matrices look much more like identity matrices. We also plot
sPLS in this figure to show just how differently this method behaves.

    * –

mixing up over trajectories: Overlap matrices can also illustrate how
different estimators from a single algorithm change with penalty parameter;
Figure 16 shows how sCCA estimates on the full dataset change with penalty
parameter. In particular, observe that the third plot shows some significant
‘mixing-up’ of the second and third directions; but otherwise the signal is
very stable.

Figure 13 shows similar biplots but now in 3D, representing the first 3
canonical variates. Again we have one plot for each regularised CCA method,
and with variates registered to the central gCCA estimate. Our main point is
that these three plots are remarkably similar. It is expected that the three
plots should be similar given the overlap matrices. However it is remarkable
quite how similar they are; moreover, now we see that much local structure is
preserved as well as global structure. One might think of this as an illusion
of being in 3D, but we suggest it is genuinely helpful. First note that the
variables are fairly evenly spread around sphere: many of the variables have
significant correlation with all 3 of the first 3 variate estimates The result
of the extra dimension is that there are very fewer variables which are very
close together than in the 2D case. A small amount of randomness will
therefore keep the more local relationships intact.

(a)

(b)

(c)

Figure 13: 3D ‘biplots’ for microbiome dataset for rCCA, gCCA, sCCA methods
and penalty parameters 0.05, 0.01, 0.001 respectively. The rCCA and gCCA
solutions have been registered with the rCCA up to orthogonal transformations
in variate space. Variates in question are from K0s, but only the C0 variables
are plotted, with colours consistent between the three plots but otherwise
arbitrary. Figure 14: Squared overlap matrices comparing multiple gCCA
estimates for different folds of the microbiome dataset, all using penalty
parameter $\alpha=0.0059$. The y-axes correspond to weights trained on the
full dataset, while the x-axes corresponds to weights fitted only on training
data corresponding to folds $\nu=0,1,2$ respectively; but crucially, these
weights are compared ‘in variate space’ by using the full data set (i.e. take
$\textbf{Z}=\mathbf{X}\hat{U}_{K},\quad\textbf{W}=\mathbf{X}\hat{U}^{(-\nu)}_{K}$
in notation of Section 7.2).

Figure 15: Squared overlap matrices comparing sPLS, sCCA, gCCA, rCCA
estimators with maximal CV correlation (r2s5-cv), as in Figure 6 but now the
sPLS, sCCA and rCCA estimates have been registered to the gCCA estimates up to
orthogonal matrices

.

Figure 16: Overlap matrix of correlations between multiple sCCA canonical
variate estimates for multiple penalty parameters. The y-axis represents
estimates using the r2s5-cv-optimal penalty; the three different x-axes
represents that same r2s5-cv-optimal penalty and penalties a factor of 3
either side. Figure 17: Sin$\Theta$ distance for top-3 weight subspaces for
Microbiome dataset; axes are the same concatenated trajectories for sCCA,
gCCA, rCCA respectively

###  I.2 Nutrimouse

We can now justify the following observations alluded to in the main text.

  * •

correlations: see Figures 18 and 19. Observe that there are many successive
pairs with significant CV signal; moreover, genuine CCA methods capture
similar signal, while sPLS captures significantly less signal.

  * •

very small penalty parameters: appear near optimal in Figure 18 for sCCA and
gCCA.

  * •

sPLS behaves very differently to the genuine CCA methods, as illustrated by
the square overlap matrices in Figure 21.

  * •

gCCA very similar to rCCA: this is clearest to see via trajectory comparison
matrices Figure 20; indeed even the weight subspaces estimated are very
similar. But this also manifests in all the other plots: for small penalty
parameters, gCCA and rCCA have very similar correlation and stability
properties in Figure 18; there are striking similarities in the overall shape
of the correlation decay plots, and also specific arrangement of points, for
example $k=4,6$ both have extremely high test correlation on fold 2 of 5;
finally, the relevant overlap matrix in Figure 21 is very near the identity
(and this was without any registration).

  * •

sCCA is less stable: for near-optimal CV penalty parameters. This is both for
stability over folds Figure 18 and stability with respect to tuning parameter
Figure 20 (higher distances between sCCA algorithms than between the rCCA or
gCCA algorithms, especially in weight space).

(a)

(b)

Figure 18: Top: CV sums of correlations as function of regularisation path for
the four methods on the Nutrimouse dataset; error bars for the aggregated
quantities. Bottom: stability both in weight space and variate space along the
same trajectories. Figure 19: 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.

(a)

(b)

Figure 20: Top: Nutrimouse trajectory comparison top-3 variate subspace
distances.
Bottom: Nutrimouse trajectory comparison top-3 weight subspace distances.

Figure 21: Squared overlap matrices for Nutrimouse dataset comparing the sPLS,
sCCA, gCCA, rCCA estimators with maximal CV correlation (r2s5-cv); analogue of
Figure 6.

.

###  I.3 BreastData

We can now justify the following observations for the BreastData dataset
alluded to in the main text.

  * •

large correlations: see Figures 23 and 24. Again there are many successive
pairs with significant CV signal. Again, the genuine CCA methods capture
similar signal. At first glance sPLS appears to capture similar signal for
very small penalty parameters; on closer inspection this is due to the non-
orthogonality, and is significantly different when a subspace notion is used –
see Figure 27.

  * •

rCCA is almost invariant to penalty parameter: as suggested by Figure 23 and
made clear in Figure 25; it is also fairly stable to sample splitting — see
Figure 23.

  * •

sCCA and sPLS very different to rCCA: firstly there are some similarities
between the sCCA and sPLS solutions, particularly in variate space Figures 25
and 26. However, these are very different to the (more stable) rCCA estimates.
This manifests very clearly in the corresponding biplots Figure 22.

Figure 22 shows 3D biplots for rCCA and sCCA, which illustrate two main
points. Firstly, they show how rCCA and sCCA can capture very different sorts
of structure: the sCCA estimates are fairly ‘local’ with only very few
variables highly correlated with the canonical variates, while the rCCA
solutions are more ‘global’ with a much larger number of variables with
comparable large correlations with the canonical variates. Secondly, they show
that in this case the biplot structure reflects biologically meaningful
information: the biplots successfully recover the chromosomal information
(colors).

(a)

(b)

Figure 22: 3D ‘biplots’ for BreastData for rCCA (left) and sCCA (right) with
penalty parameters 0.9, 0.072 respectively. Only the DNA copy number
($X$-view) variables are plotted, coloured by their chromosome.

(a)

(b)

Figure 23: Top: CV sums of squared correlations as function of regularisation
path for sPLS, sCCA, rCCA on the BreastData dataset; error bars for the
aggregated quantities. Bottom: stability in variate space along those same
trajectories. Figure 24: CV correlations (colours) and average values (black)
for successive direction estimates using sPLS, sCCA, rCCA on microbiome
dataset; in each case r2s5-cv optimal penalty parameters were used.

(a)

(b)

Figure 25: Top: Breastdata trajectory comparison top-3 variate subspace
distances.
Bottom: Breastdata trajectory comparison top-3 weight subspace distances.

Figure 26: Squared overlap matrices for BreastData dataset comparing the sPLS,
sCCA, rCCA estimators with maximal CV correlation (r2s5-cv); analogue of
Figure 6.

.

(a) Figure 27: CV sums of squared subspace correlations as a function of
regularisation path for sPLS, sCCA, rCCA on the BreastData dataset; error bars
for standard deviation.

##  Appendix J Implementation of sCCA

Though the arXiv preprint Suo et al. [2017] is widely cited, it does not in
fact present a complete algorithm to solve Equation 12: they only present a
complete algorithm for the first pair of canonical direction estimates, and
suggest how to extend it to the general case of later canonical pairs. In
addition, we found that when implemented naively, the algorithm was very slow
to run, but could be sped up by an order of magnitude by interleaving the
update steps for the estimates from each view. We present this improved
algorithm in Section J.4. First though, we present the background for the CCA
objective from Suo et al. [2017] in Section J.1 and for the necessary
technique from convex optimisation [Parikh and Boyd, 2014] in Section J.2.

###  J.1 Objective for later canonical directions from Suo et al. [2017]

For the rest of this section we follow Suo et al. [2017] and write
$\mathbf{X},\mathbf{Y}$ for data matrices that have been downscaled by a
factor of $\sqrt{N}$; then we can drop all the factors of $\frac{1}{N}$ in
Equation 5 (for example writing
$\mathbf{C}_{xx}=\mathbf{X}^{\top}\mathbf{X}$). Their objective function
Equation 12 for the $k^{\text{th}}$ canonical pair can then be written as

$\begin{array}[]{cl}\underset{u,v}{\operatorname{minimize}}&-u^{\top}\mathbf{X}^{\top}\mathbf{Y}v+\tau_{x}|u|_{1}+\tau_{y}|v|_{1}+\mathbbm{1}\left\\{u:\|\mathbf{X}u\|_{2}\leq
1\right\\}+\mathbbm{1}\left\\{v:\|\mathbf{Y}v\|_{2}\leq 1\right\\}\\\ \text{
subject to
}&U^{\top}\mathbf{X}^{\top}\mathbf{X}u=0;V^{\top}\mathbf{Y}^{\top}\mathbf{Y}v=0\;.\end{array}$

where $U\in\mathbb{R}^{p\times(k-1)},V\in\mathbb{R}^{q\times(k-1)}$ have
columns consisting of the previous canonical direction estimates (we drop
potential subscript $(k-1)$’s to avoid clutter).

This problem is biconvex; the suggested alternating convex search (ACS)
algorithm from Suo et al. [2017] iteratively fixes one variable and solves for
the other one, alternating $u,v$. Fixing $v$, we get $\hat{u}$ by solving

$\begin{array}[]{cl}\underset{u}{\operatorname{minimize}}&-u^{\top}\mathbf{X}^{\top}\mathbf{Y}v+\tau_{x}|u|_{1}+\mathbbm{1}\left\\{\mathbf{z}:\|\mathbf{z}\|_{2}\leq
1\right\\}\\\ \text{ subject to
}&\mathbf{X}u=\mathbf{z};\;U^{\top}\mathbf{X}^{\top}\mathbf{X}u=0_{k-1}\end{array}$
(84)

and fixing $u$, we get $\hat{v}$ by solving

$\begin{array}[]{cl}\underset{v}{\operatorname{minimize}}&-u^{\top}\mathbf{X}^{\top}\mathbf{Y}v+\tau_{y}|v|_{1}+\mathbbm{1}\left\\{\mathbf{z}:\|\mathbf{z}\|_{2}\leq
1\right\\}\\\ \text{ subject to
}&\mathbf{Y}v=\mathbf{z};\;V^{\top}\mathbf{Y}^{\top}\mathbf{Y}v=0_{k-1}\;.\end{array}$

In each case, the constraints can be combined as

$\left(\begin{array}[]{c}\mathbf{X}\\\
U^{\top}\mathbf{X}^{\top}\mathbf{X}\end{array}\right)u-\left(\begin{array}[]{l}I_{N}\\\
0_{(k-1)\times
N}\end{array}\right)\mathbf{z}=0,\quad\left(\begin{array}[]{c}\mathbf{Y}\\\
V^{\top}\mathbf{Y}^{\top}\mathbf{Y}\end{array}\right)v-\left(\begin{array}[]{l}I_{N}\\\
0_{(k-1)\times N}\end{array}\right)\mathbf{z}=0\;.$

This motivates the following definitions

$\displaystyle\tilde{\mathbf{X}}=\left(\begin{array}[]{c}\mathbf{X}\\\
U^{\top}\mathbf{X}^{\top}\mathbf{X}\end{array}\right),\quad\tilde{\mathbf{Y}}=\left(\begin{array}[]{c}\mathbf{Y}\\\
V^{\top}\mathbf{Y}^{\top}\mathbf{Y}\end{array}\right),\quad\tilde{I}=\left(\begin{array}[]{l}I_{N}\\\
0_{(k-1)\times N}\end{array}\right)$ (91)

We now highlight a mistake in Suo et al. [2017]. They claim

$\displaystyle-u^{\top}\mathbf{X}^{\top}\mathbf{Y}v=-u^{\top}\tilde{\mathbf{X}}^{\top}\mathbf{Y}v$
(92)

and that one can therefore use a linearised Alternating Direction Method of
Multipliers (ADMM) routine analogous to that they propose for the first
canonical pair, but replacing $\mathbf{X},\mathbf{Y}$ by the extended matrices
$\tilde{\mathbf{X}}$ and $\tilde{\mathbf{Y}}$. However, Equation 92 is
certainly false because the dimensions for the matrix multiplication
$\tilde{\mathbf{X}}^{\top}\mathbf{Y}$ are inconsistent for $k>1$: indeed,
$\tilde{\mathbf{X}}$ has $N+(k-1)$ rows while $\mathbf{Y}$ only has $N$ rows.

Instead, we follow the spirit of Suo et al. [2017] and rewrite the $u$-update
from Equation 84 in this notation as

$\underset{v}{\operatorname{minimize}}-u^{\top}\mathbf{X}^{\top}\mathbf{Y}v+\tau_{x}|u|_{1}+\mathbbm{1}\left\\{\mathbf{z}:\|\mathbf{z}\|_{2}\leq
1\right\\}\quad\text{ subject to
}\quad\tilde{\mathbf{X}}u-\tilde{I}\mathbf{z}=0\;.$ (93)

We see that we need a version of linearised ADMM that can deal with the affine
constraint $\tilde{\mathbf{X}}u-\tilde{I}\mathbf{z}=0$. Unfortunately, the
standard reference [Parikh and Boyd, 2014] for linearised ADMM does not handle
constraints of this form out-of-the-box; we consider how to extend linearised
ADMM to this setting in the following subsection.

###  J.2 Linearised ADMM with affine constraints

Throughout this subsection we work with the notation from the monographs [Boyd
and Vandenberghe, 2004, Parikh and Boyd, 2014]. Section 4.4.2 of Parikh and
Boyd [2014] presents linearised ADMM with constraints of the form $Ax=z$. A
more general form of ADMM given in Section 3 of Boyd [2010] handles affine
constraints of the form $Ax+Bz=c$. We now work through the derivation of
linearised ADMM algorithm in this more general case of affine constraint. In
the next subsection (Section J.3), we then use this derivation to obtain an
algorithm for the $u$-update objective Equation 93.

Consider the problem

$\begin{array}[]{ll}\operatorname{minimize}&f(x)+g(z)\\\ \text{ subject to
}&Ax+Bz=c\end{array}$

with variables $x$ and $z$. The augmented Lagrangian for this problem can be
written as

$\displaystyle
L_{\rho}^{(y)}(x,z,y)=f(x)+g(z)+y^{\top}(Ax+Bz-c)+\frac{\rho}{2}\|Ax+Bz-c\|_{2}^{2}\;,$

where $y\in\mathbb{R}^{m}$ is a dual variable and $\rho=1/\lambda$. Since the
scale of $y$ is immaterial, we work with $\xi\vcentcolon=y/\rho$ for the rest
of this section999Note that this substitution is not made explicit in our main
reference of Parikh and Boyd [2014].; this transforms the objective to

$\displaystyle
L_{\rho}(x,z,\xi)=f(x)+g(z)+\rho\left(\xi^{\top}(Ax+Bz-c)+\frac{1}{2}\|Ax+Bz-c\|_{2}^{2}\right)\;,$
(94)

Usual ADMM would give updates

$\displaystyle x^{k+1}$
$\displaystyle:=\underset{x}{\operatorname{argmin}}\;L_{\rho}\left(x,z^{k},y^{k}\right)$
(95) $\displaystyle z^{k+1}$
$\displaystyle:=\underset{z}{\operatorname{argmin}}\;L_{\rho}\left(x^{k+1},z,y^{k}\right)$
$\displaystyle\xi^{k+1}$
$\displaystyle:=\xi^{k}+\left(Ax^{k+1}+Bz^{k+1}-c\right)\;.$

In linearised ADMM, we modify the usual $x$-update by replacing the
$\frac{\rho}{2}\left\|Ax+Bz^{k}-c\right\|_{2}^{2}$ term in $L_{\rho}$ by

$\rho\left(A^{\top}Ax^{k}+A^{\top}Bz^{k}\right)^{\top}x+\frac{1}{2\mu}\left\|x-x^{k}\right\|_{2}^{2}\;.$

where $\mu$ is an additional tuning parameter that satisfies
$0<\mu\leq{\lambda}/{\lVert A\rVert_{\textrm{op}}^{2}}$. Plugging this into
Equation 95 and using standard properties of proximal operators [Parikh and
Boyd, 2014] gives the alternative $x$-update defined by

$\displaystyle x^{k+1}$
$\displaystyle=\underset{x}{\operatorname{argmin}}\;f(x)+\rho(Ax)^{\top}\xi^{k}+\rho\left(A^{\top}Ax^{k}+A^{\top}Bz^{k}\right)^{\top}x+\frac{1}{2\mu}\left\|x-x^{k}\right\|_{2}^{2}$
$\displaystyle=\underset{x}{\operatorname{argmin}}\;\mu
f(x)+\frac{\mu}{\lambda}\,\langle
x,A^{\top}(Ax^{k}+Bz^{k}+\xi^{k})\rangle+\frac{1}{2}\left\|x-x^{k}\right\|_{2}^{2}\vphantom{\biggl{(}\biggr{)}}$
$\displaystyle=\operatorname{prox}_{\mu
f}\left(x^{k}-\frac{\mu}{\lambda}A^{\top}(Ax^{k}+Bz^{k}+\xi^{k})\right)$

Under the simpler setting where $Ax=z$, the $z$-update can also be simply
expressed as a proximal operator, because the augmented Lagrangian Equation 94
because the quadratic part of $z$ has isotropic Hessian. This does not hold
for general $B$, but it does hold for our special case $B=\tilde{I}$ from
Equation 91. We write out the update steps explicitly in the next subsection
for this special case that we need to implement sCCA.

###  J.3 Applying linearised ADMM to obtain updates for Equation 93

We can now translate the method from Section J.2 to our setting of Section J.1
in order to obtain update steps to solve the optimisation problem Equation 93.
We identify

$(x,z,A,B,c)\Longleftrightarrow(u,\mathbf{z},\tilde{\mathbf{X}},-\tilde{I},0)$

and

$f:u\mapsto-u^{\top}\mathbf{X}^{\top}\mathbf{Y}v+\tau_{x}\|u\|_{2},\qquad
g:\mathbf{z}\mapsto\mathbbm{1}\left\\{\mathbf{z}:\|\mathbf{z}\|_{2}\leq
1\right\\}\;.$

To show that the $\mathbf{z}$-update can be written in terms of a proximal
operator, we exploit the block structure of $\tilde{\mathbf{X}}$ and
$\tilde{I}$ to obtain

$\displaystyle\mathbf{z}^{k+1}$
$\displaystyle=\underset{\mathbf{z}}{\operatorname{argmin}}\;\mathbbm{1}\left\\{\mathbf{z}:\lVert\mathbf{z}\rVert_{2}\leq
1\right\\}+\frac{\rho}{2}\;\lVert\tilde{I}\mathbf{z}-\tilde{\mathbf{X}}u^{k+1}-\xi_{k}\rVert_{2}^{2}$
$\displaystyle=\underset{\mathbf{z}}{\operatorname{argmin}}\;\mathbbm{1}\left\\{\mathbf{z}:\lVert\mathbf{z}\rVert_{2}\leq
1\right\\}+\frac{\rho}{2}\;\lVert\mathbf{z}-\mathbf{X}u^{k+1}-\tilde{I}^{\top}\xi^{k}\rVert_{2}^{2}\;.$

Note that $\tilde{I}^{\top}\xi^{k}$ consists precisely of the first $N$
elements of $\xi^{k}$, so is very easy to compute.

Finally we obtain the updates

$\displaystyle u^{k+1}\leftarrow\operatorname{prox}_{\mu
f}\left(u^{k}-\frac{\mu}{\lambda}\tilde{\mathbf{X}}^{\top}\left(\tilde{\mathbf{X}}u^{k}-\tilde{I}\mathbf{z}^{k}+\xi^{k}\right)\right)$
(96) $\displaystyle\mathbf{z}^{k+1}\leftarrow\operatorname{prox}_{\lambda
g}\left(\mathbf{X}u^{k+1}+\tilde{I}^{\top}\xi^{k}\right)$
$\displaystyle\xi^{k+1}\leftarrow\xi^{k}+\mathbf{X}u^{k+1}-\mathbf{z}^{k+1}\;,$

where the proximal operators are defined as follows, following Suo et al.
[2017]. Write $c=\mathbf{X}^{\top}\mathbf{Y}v$ (the gradient of the objective
function with respect to $u$ for a fixed $v$). Then $f(x)$ involves (element-
wise) soft-thresholding and $g(x)$ is a projection to the unit ball defined by

$\displaystyle\operatorname{prox}_{\mu f}(x)$
$\displaystyle=\begin{cases}x+\mu c-\mu\tau_{x}&\text{ if }x+\mu
c>\mu\tau_{x}\\\ x+\mu c+\mu\tau_{x}&\text{ if }x+\mu c<-\mu\tau_{x}\\\
0&\text{ else,}\end{cases}$ $\displaystyle\operatorname{prox}_{\lambda
g}(\mathbf{z})$ $\displaystyle=\begin{cases}\mathbf{z}&\text{ if
}\|\mathbf{z}\|_{2}\leq 1\\\ \frac{\mathbf{z}}{\|\mathbf{z}\|_{2}}&\text{
else,}\end{cases}$

Note that the soft-thresholding depends on $c$ as well as $\tau_{x}$, which is
how the dependence on $v$ manifests in the update steps.

The analogous step for the $v$-optimisation for a given $u$ can be obtained by
symmetry/relabelling, as we exploit in the following subsection.

###  J.4 Complete algorithm with computational improvements

We propose Algorithm 3 to solve Equation 12. This applies the updates from
Section J.3 in a framework very similar to the original Alternating Convex
Search (ACS) procedure proposed by Suo et al. [2017], but with two small but
significant algorithmic changes:

  * •

Early stopping of the linearised ADMM sub-routines: for a fixed $v$, update
$u$ by running the linearised ADMM updates Equation 96 for some small number
of steps $n_{\text{steps\\_ADMM}}$ rather than to convergence.

  * •

Recycling optimal dual variables: after each of these sets of linearised ADMM
steps, record the values of the optimal dual variables, and use these for
initialising the subsequent set of linearised ADMM steps.

We found that these two changes reduced the number of total number of
linearised ADMM update steps required for convergence of the weight estimates
$u,v$ by an order of magnitude (on the real-world datasets of interest to us).

We now outline why these changes do not weaken the theoretical properties of
the algorithm. Suppose $u,v$ are returned by the algorithm. This requires
convergence in the outer loop, which implies convergence in the inner loop.
Since the inner loop consists of the linearised ADMM updates, this only occurs
when $u$ is the global minimum of Equation 93 for the given $v$ and vice
versa. Therefore $u,v$ must be some local minimum of the biconvex objective.
This is essentially all one could say about solutions found by the original
Alternating Convex Search (ACS) algorithm, so our algorithm has similar
theoretical justification.

A more careful convergence analysis may be possible by adapting the Lyapunov
arguments from zhang_unified_2011[Theorem 4.2]. However, the details become
involved, so we leave this for an interested reader.

We are deliberately vague about initialisation in Algorithm 3. For our
numerical experiments, we use an implementation based off ideas from Suo et
al. [2017], that starts by taking a SVD of a soft-thresholded version of the
between view covariance matrix $\mathbf{X}^{\top}\mathbf{Y}$. However, we
found that various schemes led to similar performance, and so we leave full
details to our code implementation at https://github.com/W-L-W/RegularisedCCA.

Algorithm 3 Sparse Canonical Correlation Analysis (sCCA) Next Pair Calculation

1:procedure sCCA-next-
pair($\mathbf{X},\mathbf{Y},U,V,\tau_{x},\tau_{y};\lambda,n_{\text{steps\\_ADMM}}$)

2: recall: $\tau_{x},\tau_{y}\in[0,\infty)$ are $\ell_{{1}}$ penalty
parameters, $\lambda\in[0,\infty)$ is an ADMM step size parameter,
$n_{\text{steps\\_ADMM}}$ is a small positive integer.

3: ensure: Dimensional consistency: for some positive integers $n,p,q,K$ we
have $\mathbf{X}\in\mathbb{R}^{N\times p},\mathbf{Y}\in\mathbb{R}^{N\times
q},U\in\mathbb{R}^{p\times(k-1)},V\in\mathbb{R}^{q\times(k-1)}$.

4: Construct $\tilde{\mathbf{X}},\tilde{\mathbf{Y}},\tilde{I}$ as in Equation
91

5: Initialise:

  * •

$u,v$ in some sensible manner

  * •

$\mathbf{z}_{u}\leftarrow\textsc{Normalise}(\mathbf{X}u),\quad\xi_{u}\leftarrow\tilde{\mathbf{X}}u-\tilde{I}\mathbf{z}_{u}$

  * •

$\mathbf{z}_{v}\leftarrow\textsc{Normalise}(\mathbf{Y}v),\quad\xi_{v}\leftarrow\tilde{\mathbf{Y}}v-\tilde{I}\mathbf{z}_{v}$

6: while not converged do

7: $u,\mathbf{z}_{u},\xi_{u}\leftarrow\textsc{lADMM-
steps}(u,\mathbf{z}_{u},\xi_{u},\tilde{\mathbf{X}},\mathbf{X},\mathbf{Y},v,\tau_{x};\,\lambda,n_{\text{steps\\_ADMM}})$

8: $v,\mathbf{z}_{v},\xi_{v}\leftarrow\textsc{lADMM-
steps}(v,\mathbf{z}_{v},\xi_{v},\tilde{\mathbf{Y}},\mathbf{Y},\mathbf{X},u,\tau_{y};\,\lambda,n_{\text{steps\\_ADMM}})$

9: end while

10: return $u,v$

11:end procedure

12:procedure lADMM-
steps($u,\mathbf{z},\xi,\tilde{\mathbf{X}},\mathbf{X},\mathbf{Y},v,\tau;\,\lambda,n_{\text{steps}}$)

13:
$\mu\leftarrow\lambda/(2\cdot\lVert\tilde{\mathbf{X}}\rVert_{\textrm{op}}^{2})$

14: $c\leftarrow\mathbf{X}^{\top}\mathbf{Y}v$

15: Define $\operatorname{prox}_{\mu f}:\mathbb{R}^{p}\to\mathbb{R}^{p}$ by
$x\mapsto\begin{cases}x+\mu c-\mu\tau&\text{ if }x+\mu c>\mu\tau\\\ x+\mu
c+\mu\tau&\text{ if }x+\mu c<-\mu\tau\\\ 0&\text{ else}\end{cases}$

16: for $\text{step}=1$ to $n_{\text{steps}}$ do

17: $u\leftarrow\operatorname{prox}_{\mu
f}\left(u-\frac{\mu}{\lambda}\tilde{\mathbf{X}}^{\top}\left(\tilde{\mathbf{X}}u-\tilde{I}\mathbf{z}+\xi\right)\right)$

18: $\mathbf{z}\leftarrow\textsc{Normalise}(\mathbf{X}u+\tilde{I}^{\top}\xi)$

19: $\xi\leftarrow\xi+\mathbf{X}u-\mathbf{z}$.

20: end for

21: return $u,\mathbf{z},\xi$

22:end procedure

23:procedure Normalise($\mathbf{z}$)

24: if $\lVert\mathbf{z}\rVert>1$ then

25: $\mathbf{z}\leftarrow\mathbf{z}/\lVert\mathbf{z}\rVert$

26: end if

27: return $\mathbf{z}$

28:end procedure