# On Support Recovery with Sparse CCA: Information Theoretic and Computational
Limits

Nilanjana Laha<EMAIL_ADDRESS>
Department of Biostatistics
Harvard university
Boston, MA 02115, USA Rajarshi Mukherjee<EMAIL_ADDRESS>
Department of Biostatistics
Harvard University
Boston, MA 02115, USA

###### Abstract

In this paper we consider asymptotically exact support recovery in the context
of high dimensional and sparse Canonical Correlation Analysis (CCA). Our main
results describe four regimes of interest based on information theoretic and
computational considerations. In regimes of “low” sparsity we describe a
simple, general, and computationally easy method for support recovery, whereas
in a regime of “high” sparsity, it turns out that support recovery is
information theoretically impossible. For the sake of information theoretic
lower bounds, our results also demonstrate a non-trivial requirement on the
“minimal” size of the non-zero elements of the canonical vectors that is
required for asymptotically consistent support recovery. Subsequently, the
regime of “moderate” sparsity is further divided into two sub-regimes. In the
lower of the two sparsity regimes, using a sharp analysis of a coordinate
thresholding (Deshpande and Montanari, 2014) type method, we show that
polynomial time support recovery is possible. In contrast, in the higher end
of the moderate sparsity regime, appealing to the “Low Degree Polynomial”
Conjecture (Kunisky et al., 2019), we provide evidence that polynomial time
support recovery methods are inconsistent. Finally, we carry out numerical
experiments to compare the efficacy of various methods discussed.

Keywords: Sparse Canonical Correlation Analysis, Minimax Support Recovery, Low
Degree Polynomials

##  1 Introduction

Canonical Correlation Analysis (CCA) is a highly popular technique to perform
initial dimension reduction while exploring relationships between two
multivariate objects. Due to its natural interpretability and success in
finding latent information, CCA has found enthusiasm across of vast canvass of
disciplines, which include, but are not limited to psychology and agriculture,
information retrieving (Hardoon et al., 2004; Rasiwasia et al., 2010; Gong et
al., 2014), brain-computer interface (Bin et al., 2009), neuroimaging (Avants
et al., 2010), genomics (Witten et al., 2009), organizational research
(Bagozzi, 2011), natural language processing (Dhillon et al., 2011; Faruqui
and Dyer, 2014), fMRI data analysis data analysis (Friman et al., 2003),
computer vision (Kim et al., 2007), and speech recognition (Arora and Livescu,
2013; Wang et al., 2015).

Early developments in the theory and applications of CCA have now been well
documented in statistical literature, and we refer the interested reader to
Anderson (2003) and references therein for further details. However, modern
surge in interests for CCA, often being motivated by data from high throughput
biological experiments (Lê Cao et al., 2009; Lee et al., 2011; Waaijenborg et
al., 2008), requires re-thinking several aspects of the traditional theory and
methods. A natural structural constraint that has gained popularity in this
regard, is that of sparsity, i.e. the phenomenon of an (unknown) collection of
variables being related to each other. In order to formally introduce the
framework of sparse CCA, we present our statistical set up next.  We shall
consider $n$-i.i.d. samples $(X_{i},Y_{i})\sim\mathbb{P}$ with
$X_{i}\in\mathbb{R}^{p}$ and $Y_{i}\in\mathbb{R}^{q}$ being multivariate mean
zero random variables with joint variance covariance matrix

$\Sigma=\begin{bmatrix}{\Sigma}_{x}&{\Sigma}_{yx}\\\
{\Sigma}_{yx}&{\Sigma}_{y}\\\ \end{bmatrix}.$ (1)

The first canonical correlation $\Lambda_{1}$ is then defined as the maximum
possible correlation between two linear combinations of $X$ and $Y$. This
definition interprets $\Lambda_{1}$ as the optimal value of the following
maximization problem:

$\displaystyle\mathmakebox[width("$\underset{\displaystyle
u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\mathrm{subject~{}to}}$")][l]{\underset{\displaystyle
u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}{\mathrm{maximize}}}\quad
u^{T}{\Sigma}_{xy}v\hfil\hfil\hfil\hfil$ (2)
$\displaystyle\mathmakebox[width("$\underset{\displaystyle\phantom{u\in\mathbb{R}^{p},v\in\mathbb{R}^{q}}}{\mathrm{subject~{}to}}$")][c]{{\mathrm{subject~{}to}}}\quad$
$\displaystyle u^{T}{\Sigma}_{x}u=v^{T}{\Sigma}_{y}v$ $\displaystyle=1$

The solutions to (2) are the vectors which maximize the correlation of the
projections of $X$ and $Y$ in those respective directions. Higher order
canonical correlations can thereafter be defined in a recursive fashion (cf.
Anderson, 1999). In particular, for $j\geq 1$, we define the $j^{\rm th}$
canonical correlation $\Lambda_{j}$ and the corresponding directions $u_{j}$
and $v_{j}$ by maximizing (2) with the additional constraint

$\displaystyle u^{T}{\Sigma}_{x}u_{l}=v^{T}{\Sigma}_{y}v_{l}=0,\quad 0\leq
l\leq j-1.$ (3)

As mentioned earlier, in many modern data examples, the sample size $n$ is
typically at most comparable to or much smaller than $p$ or $q$ – rendering
the classical CCA inconsistent and inadequate without further structural
assumptions (Cai et al., 2018; Ma et al., 2020; Bao et al., 2019). The
framework of Sparse Canonical Correlation Analysis (SCCA) (Witten et al.,
2009; Mai and Zhang, 2019), where the $u_{i}$’s and the $v_{i}$’s are sparse
vectors, was subsequently developed to target low dimensional structures (that
allows consistent estimation) when $p,q$ are potentially larger than $n$. The
corresponding sparse estimates of the leading canonical directions naturally
perform variable selection, thereby leading to recovery of their support
(Witten et al., 2009; Mai and Zhang, 2019; Waaijenborg et al., 2008; Solari et
al., 2019). It is unknown, however, under what settings, this naïve method of
support recovery, or any other method for the matter, is consistent. The
support recovery of the leading canonical directions serves an important
purpose of identifying groups of variables which explain the most linear
dependence among high dimensional random objects ($X$ and $Y$) under study –
and thereby renders crucial interpretability. Asymptotically optimal support
recovery is yet to be explored systematically in the context of SCCA – both
theoretically, and from the computational viewpoint. In fact, despite the
renewed enthusiasm on CCA, both the theoretical and applied communities have
mainly focused on the estimation of the leading canonical directions, and
relevant scalable algorithms – see e.g. Chen et al. (2013); Gao et al. (2015,
2017); Ma et al. (2020); Mai and Zhang (2019). This paper is motivated by
exploring the crucial question of support recovery in the context of SCCA
111In this paper, by support recovery, we refer to the exact recovery of the
combined support of the $u_{i}$’s (or the $v_{i}$’s)..

The problem of support recovery for SCCA naturally connects to a vast class of
variable selection problems (Wainwright, 2009; Amini and Wainwright, 2009;
Butucea et al., 2015; Butucea and Stepanova, 2017; Meinshausen and Bühlmann,
2010). The problem closest in terms of complexity turns out to be the sparse
PCA (SPCA) problem (Johnstone and Lu, 2009). Support recovery in the latter
problem is known to present interesting information theoretic and
computational bottlenecks (cf. Krauthgamer et al., 2015; Amini and Wainwright,
2009; Ding et al., 2019; Arous et al., 2020). Moreover, information theoretic
and computational issues also arise in context of SCCA estimation problem
(Chen et al., 2013; Gao et al., 2015, 2017; Mai and Zhang, 2019). In view of
the above, it is natural to expect that such information theoretic and
computational issues exist in context of SCCA support recovery problem as
well. However, the techniques used in the SPCA support recovery analysis is
not directly applicable to the SCCA problem, which poses additional challenges
due to the presence of high dimensional nuisance parameters ${\Sigma}_{x}$ and
${\Sigma}_{y}$. The main focus of our work is therefore retrieving the
complete picture of the information theoretic and computational limitations of
SCCA support recovery. Before going into details, we next present a brief
summary of our contributions, and defer the discussions on the main subtleties
to Section 3.

###  1.1 Summary of Main Results

We say a method successfully recovers the support if it achieves exact
recovery with probability tending to one uniformly over the sparse parameter
spaces defined in Section 2. In the sequel, we denote the cardinality of the
combined support of the $u_{i}$’s and the $v_{i}$’s by $s_{x}$ and $s_{y}$,
respectively. Thus $s_{x}$ and $s_{y}$ will be our respective sparsity
parameters. Our main contributions are listed below.

####  1.1.1 General methodology

In Section 3.1, we construct a general algorithm called RecoverSupp, which
leads to successful support recovery whenever the latter is information
theoretically tractable. This also serves as the first step in creating a
polynomial time procedure for recovering support in one of the difficult
regimes of the problem – see e,g. Corollary 17 which shows that RecoverSupp
accompanied by a coordinate thresholding type method recovers the support in
polynomial time in a regime that requires subtle analysis. Moreover, Theorem 2
shows that the minimal signal strength required by RecoverSupp is matches the
information theoretic limit whenever the nuisance precision matrices
${\Sigma}_{x}^{-1}$ and ${\Sigma}_{y}^{-1}$ are sufficiently sparse.

####  1.1.2 Information Theoretic and Computational Hardness as a Function of
Sparsity

As the sparsity level increases, we show that the CCA support recovery problem
transitions from being efficiently solvable, to NP hard (conjectured), and to
information theoretically impossible. According to this hardness pattern, the
sparsity domain can be partitioned into the following three regimes: (i)
$s_{x},s_{y}\lesssim\sqrt{n}$, (ii) $\sqrt{n}\lesssim s_{x},s_{y}\lesssim
n/\log(p+q)$, and (iii) $s_{x},s_{y}\gtrsim n/\log(p+q)$. We describe below
the distinguishing behaviours of these three regimes, which is consistent with
the sparse PCA scenario.

  * •

We show that when $s_{x},s_{y}\lesssim\sqrt{n/\log(p+q)}$ (“easy regime”),
polynomial time support recovery is possible, and well-known consistent
estimators of the canonical correlates (Mai and Zhang, 2019; Gao et al., 2017)
can be utilized to that end. When $\sqrt{n/\log(p+q)}\lesssim
s_{x},s_{y}\lesssim\sqrt{n}$ (“difficult regime”), we show that a coordinate
thresholding type algorithm (inspired by Deshpande and Montanari, 2014)
succeeds provided $p+q\asymp n$. We call the last regime “difficult” because
existing estimation methods like COLAR (Gao et al., 2017) or SCCA (Mai and
Zhang, 2019) are yet to be shown to have valid statistical guarantees in this
regime – see Section 3.1 and Section 3.4 for more details.

  * •

In Section 3.3, we show that when $\sqrt{n}\lesssim s_{x},s_{y}\lesssim
n/\log(p+q)$ (“hard regime”), support recovery is computationally hard subject
to the so called “low degree polynomial conjecture” recently popularized by
Hopkins and Steurer (2017); Hopkins (2018); Kunisky et al. (2019). Of course
this phenomenon is observable only when $p,q\gtrsim n$, because otherwise, the
problem would be solvable by the ordinary CCA analysis (Bao et al., 2019; Ma
and Yang, 2021). Our findings are consistent with the conjectured
computational barrier in context of SCCA estimation problem (Gao et al.,
2017).

  * •

When $s_{x},s_{y}\gtrsim n/\log(p+q)$, we show that support recovery is
information theoretically impossible (see Section 3.2).

####  1.1.3 Information Theoretic Hardness as a Function of Minimal Signal
Strength

In context of support recovery, the signal strength is quantified by

$\texttt{Sig}_{x}=\min_{k\in[p]}\max_{i\in[r]}|(u_{i})_{k}|\quad\text{and}\quad\texttt{Sig}_{y}=\min_{k\in[q]}\max_{i\in[r]}|(v_{i})_{k}|.$

Generally, support recovery algorithms require the signal strength to be above
some threshold. As a concrete example, the detailed analyses provided in
(Amini and Wainwright, 2009; Deshpande and Montanari, 2014; Krauthgamer et
al., 2015) all is based on the non-zero principal component elements being
$\pm 1/\sqrt{\rm sparsity}$. To the best of our knowledge, prior to our work,
there was no result in the PCA/CCA literature on the information theoretic
limit of the minimal signal strength.

  * •

In Section 3.2, we show that $\texttt{Sig}_{x}\gtrsim\sqrt{\log(p-s_{x})/n}$
and $\texttt{Sig}_{y}\gtrsim\sqrt{\log(q-s_{y})/n}$ is a necessary requirement
for successful support recovery.

###  1.2 Notation

For a vector $x\in\mathbb{R}^{p}$, we denote its support $D(x)$ by
$D(x)=\\{i:x_{i}\neq 0\\}$. We will overload notation, and for a matrix
$A\in\mathbb{R}^{p\times q}$, we will denote by $D(A)$ the indexes of the non-
zero rows of $A$. By an abuse of notation, sometimes we will refer to $D(A)$
as the support of $A$ as well. When $A\in\mathbb{R}^{p\times q}$ and
$\alpha\in\mathbb{R}^{p}$ are unknown parameters, generally the estimator of
their supports will be denoted by $\widehat{D}(A)$ and $\widehat{D}(\alpha)$,
respectively. We let $\mathbb{N}$ denote the set of all positive numbers, and
write $\mathbb{Z}$ for the set of all natural numbers $\\{0,1,2,\ldots,\\}$.
For any $n\in\mathbb{N}$, We let $[n]$ denote the set $\\{1,\ldots,n\\}$. For
any finite set $\mathcal{A}$, we denote its cardinality by $|\mathcal{A}|$.
Also, we let $1\\{\mathcal{A}\\}$ be the indicator of the event $\mathcal{A}$.

We let $\|\cdot\|_{k}$ be the usual $l_{k}$ norm in $\mathbb{R}^{k}$ for
$k\in\mathbb{Z}$. In particular, we let $\|x\|_{0}$ denote the number of non-
zero elements of a vector $x\in\mathbb{R}^{p}$. For any probability measure
$\mathbb{P}$ on the Borel sigma field of $\mathbb{R}^{p}$, we take
$L_{2}(\mathbb{P})$ be the set of all measurable functions
$f:\mathbb{R}^{p}\mapsto\mathbb{R}$ such that
$\|f\|_{L_{2}(\mathbb{P})}=\sqrt{\int f^{2}d\mathbb{P}}<\infty$. The
corresponding $L_{2}(\mathbb{P})$ inner product will be denoted by
$\langle\cdot,\cdot\rangle_{L_{2}(\mathbb{P})}$. We denote the operator norm
and the Frobenius norm of a matrix $A\in\mathbb{R}^{p\times q}$ by
$\|A\|_{op}$ and $\|A\|_{F}$, respectively. For $k\in\mathbb{N}$, we define
the norms $\|A\|_{k,\infty}=\max_{j\in[q]}\|A_{j}\|_{k}$ and
$\|A\|_{\infty,k}=\max_{i\in[q]}\|A_{i*}\|_{k}$. The maximum and minimum
eigenvalue of a square matrix $A$ will be denoted respectively by
$\Lambda_{max}(A)$ and $\Lambda_{min}(A)$. We let $A_{i*}$ and $A_{j}$ denote
the i-th row and $j$-th column of $A$, respectively. Also, we let $s(A)$
denote the maximum number of non-zero entries in any column of $A$, i.e.
$s(A)=\max_{j\in[q]}\|A_{j}\|_{0}$.

The results in this paper are mostly asymptotic (in $n$) in nature and thus
require some standard asymptotic notations. If $a_{n}$ and $b_{n}$ are two
sequences of real numbers then $a_{n}\gg b_{n}$ (and $a_{n}\ll b_{n}$) implies
that ${a_{n}}/{b_{n}}\rightarrow\infty$ (and ${a_{n}}/{b_{n}}\rightarrow 0$)
as $n\rightarrow\infty$, respectively. Similarly $a_{n}\gtrsim b_{n}$ (and
$a_{n}\lesssim b_{n}$) implies that
$\liminf_{n\rightarrow\infty}{{a_{n}}/{b_{n}}}=C$ for some $C\in(0,\infty]$
(and $\limsup_{n\rightarrow\infty}{{a_{n}}/{b_{n}}}=C$ for some
$C\in[0,\infty)$). Alternatively, $a_{n}=o(b_{n})$ will also imply $a_{n}\ll
b_{n}$ and $a_{n}=O(b_{n})$ will imply that $\limsup_{n\rightarrow\infty}\
a_{n}/b_{n}=C$ for some $C\in[0,\infty)$). We will write
$a_{n}=\tilde{\Phi}(b_{n})$ to indicate $a_{n}$ and $b_{n}$ are asymptotically
of the same order up to a poly-log term. Finally, in our mathematical
statements, $C$ and $c$ will be two different generic constants which can vary
from line to line.

##  2 Mathematical Formalism

We define the rank of ${\Sigma}_{xy}$ by $r$. It can be shown that exactly $r$
canonical correlations are positive and the rest are zero in the model (2). We
will consider the matrices $U=[u_{1},\ldots,u_{r}]$ and
$V=[v_{1},\ldots,v_{r}]$. From (2) and (3), it is not hard to see that
$U^{T}{\Sigma}_{x}U=I_{p}$ and $V^{T}{\Sigma}_{y}V=I_{q}$. The indexes of the
nonzero rows of $U$ and $V$, respectively, are the combined support of the
$u_{i}$’s and the $v_{i}$’s. Since we are interested in the recovery of the
latter, it will be useful for us to study $U$ and $V$. To that end, we often
make use of the following representation connecting ${\Sigma}_{xy}$ to $U$ and
$V$ (Anderson, 2003):

${\Sigma}_{xy}={\Sigma}_{x}U\Lambda
V^{T}{\Sigma}_{y}=\sum_{i=1}^{r}\Lambda_{i}u_{i}v_{i}^{T}.$ (4)

To keep our results straightforward, we restrict our attention to a particular
model $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ throughout, defined as follows.

###### Definition 1

Suppose $(X,Y)\sim\mathbb{P}$. Let $\mathcal{B}>1$ be a constant. We say
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ if

  * A1

(Sub-Gaussian) $X$ and $Y$ are sub-Gaussian random vectors 222See e.g.
Vershynin (2018)., with joint covariance matrix $\Sigma$ as defined in (1).
Also $\text{rank}({\Sigma}_{xy})=r$.

  * A2

Recall the definition of the canonical correlation $\Lambda_{i}$’s from (3).
Note that by definition, $0\leq\Lambda_{r}\leq\ldots\leq\Lambda_{1}$. For
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$, $\Lambda_{r}$
additionally satisfies $\Lambda_{r}\geq 1/\mathcal{B}$.

  * A3

(Sparsity) The number of nonzero rows of $U$ and $V$ are $s_{x}$ and $s_{y}$,
respectively, that is $s_{x}=|\cup_{i=1}^{r}D(u_{i})|$ and
$s_{y}=|\cup_{i=1}^{r}D(v_{i})|$. Here $U$ and $V$ are as defined in (4).

  * A4

(Bounded eigenvalue)
$1/\mathcal{B}<\Lambda_{min}({\Sigma}_{y}),\Lambda_{min}({\Sigma}_{y}),\Lambda_{max}({\Sigma}_{x}),\Lambda_{max}({\Sigma}_{y})<\mathcal{B}$.

  * A5

(Positive eigen-gap) $\Lambda_{i}-\Lambda_{i-1}\geq\mathcal{B}^{-1}$ for
$i=2,\ldots,r$.

Sometimes we will consider a sub-model of
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ where each
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ is Gaussian. This model
will be denoted by $\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$, where “$G$”
stands for the Gaussian assumption. Some remarks on the modelling assumptions
A1—A5 are in order, which we provide next.

  1. A1.

We begin by noting that we do not require $X$ and $Y$ to be jointly sub-
Gaussian. Moreover, the individual sub-Gaussian assumption itself is common in
the $s_{x},s_{y}\lesssim\sqrt{n}/\log(p+q)$ regime in the SCCA literature (Mai
and Zhang, 2019; Gao et al., 2017). For the sharper analysis in the difficult
regime ($\sqrt{n/\log{(p+q)}}\lesssim s_{x},s_{y}\lesssim\sqrt{n}$), our proof
techniques require the Gaussian model $\mathcal{P}_{G}$ – which is in parallel
with Deshpande and Montanari (2014)’s treatment of the sparse PCA in the
corresponding difficult regime. In general, the Gaussian spiked model
assumption in sparse PCA goes back to Johnstone (2001), and is common in the
PCA literature (Amini and Wainwright, 2009; Krauthgamer et al., 2015).

  2. A2-A4.

These assumptions are standard in the analysis of canonical correlations (Mai
and Zhang, 2019; Gao et al., 2017).

  3. A5.

This assumption concerns the gap between consecutive canonical correlation
strengths. However we refer to this gap as “Eigengap” because of its
similarity with the Eigengap in the sparse PCA literature (cf. Deshpande and
Montanari, 2014; Janková and van de Geer, 2018). This assumption is necessary
for the estimation of the $i$-th canonical covariates. Indeed, if
$\lambda_{i}=\lambda_{i+1}$ then there is no hope of estimating the $i$-th
canonical covariates because they are not identifiable, and so support
recovery also becomes infeasible. This assumption can be relaxed to requiring
only $k$ many $\lambda_{i}$’s to be strictly larger than $\lambda_{i-1}$’s
where $k\leq r$. In this case, we can recover the support of only the first
$k$ canonical covariates.

In the following sections, we will denote the preliminary estimators of $U$
and $V$ by $\widehat{U}$ and $\widehat{V}$, respectively. The columns of
$\widehat{U}$ and $\widehat{V}$ will be denoted by $\widehat{u}_{n,i}$ and
$\widehat{v}_{n,i}$ ($i\in[r]$), respectively. Therefore $\widehat{u}_{n,i}$
and $\widehat{v}_{n,i}$ will stand for the corresponding preliminary
estimators of $u_{i}$ and $v_{i}$. In case of CCA, the $u_{i}$’s and $v_{i}$’s
are identifiable only up to a sign flip. Hence, they are also estimable only
up to a sign flip. Finally, we denote the empirical estimates of
${\Sigma}_{x}$, ${\Sigma}_{y}$, and ${\Sigma}_{xy}$, by
$\widehat{\Sigma}_{n,x}$, $\widehat{\Sigma}_{n,y}$, and
$\widehat{\Sigma}_{n,xy}$, respectively – which will often be appended with
superscripts to denote their estimation through suitable sub-samples of the
data 333e.g. $\widehat{\Sigma}_{n,x}^{(j)}$, $\widehat{\Sigma}_{n,y}^{(j)}$,
and $\widehat{\Sigma}_{n,xy}^{(j)}$ will stand for the empirical estimators
created from the $j^{th}$-equal split of the data.. Finally, we let
$C_{\mathcal{B}}$ denote a positive constant which depends on $\mathbb{P}$
only through $\mathcal{B}$, but can vary from line to line.

##  3 Main Results

We divide our main results into the following parts based on both statistical
and computational difficulties of different regimes. First in Section 3.1 we
present a general method and associated sufficient conditions for support
recovery. This allows us to elicit a sequence of questions regarding necessity
of the conditions and remaining gaps both from statistical and computational
perspectives. Our subsequent sections are devoted to answering these very
questions. In particular, in Section 3.2 we discuss information theoretic
lower bounds followed by evidence for statistical-computational gaps in
Section 3.3. Finally, we close a final computational gap in asymptotic regime
through sharp analysis of a special coordinate-thresholding type method in
Section 3.4.

###  3.1 A Simple and General Method:

We begin with a simple method for estimating the support, which readily
establishes the result for the easy regime, and sets the directions for the
investigation into other more subtle regimes. Since the estimation of $D(U)$
and $D(V)$ are similar, we focus only on the estimation of $D(V)$ for the time
being.

Suppose $\widehat{V}$ is a row sparse estimator of $V$. The non-zero indexes
of $\widehat{V}$ is the most intuitive estimator of $D(V)$. Such an
$\widehat{V}$ is also easily attainable because most estimators of the
canonical directions in high dimension are sparse (cf. Chen et al., 2013; Gao
et al., 2017; Mai and Zhang, 2019, among others). Although we have not yet
been able to show the validity of this apparently “naïve” method, we provide
numerical results in Section 4 to explore its finite sample performance.
However, a simple method can refine these initial estimators, to often
optimally recover the support $D(V)$. We now provide the details of this
method and derive its asymptotic properties.

To that end, suppose we have at our disposal an estimating procedure for
${\Sigma}_{y}^{-1}$, which we generically denote by $\widehat{\Omega}_{n}$ and
an estimator $\widehat{U}\in\mathbb{R}^{p\times r}$ of $U$. We split the
sample in two equal parts, and compute $\widehat{U}^{(1)}$ and
$\widehat{\Omega}_{n}^{(1)}$ from the first part of the sample, and the
estimator $\widehat{\Sigma}_{n,xy}^{(2)}$ from the second part of the sample.
Define
$\widehat{V}^{clean}=\widehat{\Omega}_{n}^{(1)}\widehat{\Sigma}_{n,xy}^{(2)}\widehat{U}^{(1)}$.
Our estimator of $D(V)$ is then given by

$\widehat{D}(V):=\\{i\in[q]:|\widehat{V}^{clean}_{ij}|>\texttt{cut}\text{ for
some }j\in[r]\\},$ (5)

where cut is a pre-specified cut-off or threshold. We will discuss later how
to choose cut efficiently. The resulting algorithm, detailed as Algorithm 1
for convenience, will be referred as RecoverSupp from now on. RecoverSupp is
similar in spirit to the “cleaning” step in the sparse PCA support recovery
literature (cf. Deshpande and Montanari, 2014).

Algorithm 1 RecoverSupp
$(\widehat{U}^{(1)},\widehat{\Omega}_{n}^{(1)},\widehat{\Sigma}_{n,xy}^{(2)},\texttt{cut},r)$:
suppoet recovery of $V$

0:

  1. 1.

Preliminary estimators $\widehat{U}^{(1)}$ and $\widehat{\Omega}_{n}^{(1)}$ of
$U$ and ${\Sigma}_{y}^{-1}$, respectively, based on sample
$O_{1}=(x_{i},y_{i})_{i=1}^{[n/2]}$.

  2. 2.

Estimator $\widehat{\Sigma}_{n,xy}^{(2)}$ of ${\Sigma}_{xy}$ based on sample
$O_{2}=(x_{i},y_{i})_{i=[n/2]+1}^{n}$.

  3. 3.

Threshold level $\texttt{cut}>0$ and rank $r\in\mathbb{N}$.

0: $\widehat{D}(V)$, an estimator of ${D(V)}$.

  1. 1.

Cleaning:
$\widehat{V}^{clean}\leftarrow\widehat{\Omega}_{n}^{(1)}{{\widehat{\Sigma}_{n,yx}}}^{(2)}\widehat{U}^{(1)}$.

  2. 2.

Threshold: Compute $\widehat{D}(V)$ as in (5).

Return: $\widehat{D}(V)$.

It turns out that, albeit being so simple, RecoverSupp has desirable
statistical guarantees provided $\widehat{U}^{(1)}$ and
$\widehat{\Omega}_{n}^{(1)}$ are reasonable estimators of $U$ and
${\Sigma}_{y}^{-1}$, respectively. These theoretical properties of RecoverSupp
, and the hypotheses and queries generated thereof, lays out the roadmap for
the rest of our paper. However, before getting into the detailed theoretical
analysis of RecoverSupp , we state a $l_{2}$-consistency condition on
$\widehat{u}_{n,i}$ and $\widehat{v}_{n,i}$’s, where we remind the readers
that we let $\widehat{u}_{n,i}$ and $\widehat{v}_{n,i}$ denote the $i$-th
columns of $\widehat{V}$ and $\widehat{U}$, respectively. Recall also that the
$i$-th columns of $U$ and $V$ are denoted by $u_{i}$ and $v_{i}$,
respectively.

######  Condition 1 ($l_{2}$ consistency )

There exists a function
$\texttt{Err}\equiv\texttt{Err}:(n,p,q,s_{x},s_{y},\mathcal{B})\mapsto\mathbb{R}$
so that $|\texttt{Err}|<1/(2\mathcal{B}\sqrt{r})$ and the estimators
$\widehat{u}_{n,i}$ and $\widehat{v}_{n,i}$ of $u_{i}$ and $v_{i}$ satisfy

$\max_{i\in[r]}\min_{w\in\\{\pm
1\\}}\bigg{|}(w\widehat{u}_{n,i}-u_{i})^{T}{\Sigma}_{x}(w\widehat{u}_{n,i}-u_{i})\bigg{|}<\texttt{Err}^{2},$
$\max_{j\in[r]}\min_{w\in\\{\pm
1\\}}\bigg{|}(w\widehat{v}_{n,i}-v_{i})^{T}{\Sigma}_{y}(w\widehat{v}_{n,i}-v_{i})\bigg{|}<\texttt{Err}^{2}$

with $\mathbb{P}$ probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},B)$ .

We will discuss the estimators which satisfy Condition 1 later. Theorem 2 also
requires the signal strength $\texttt{Sig}_{y}$ to be at least of the order
$\epsilon_{n}=\xi_{n}\sqrt{\log(p+q)s({\Sigma}_{y}^{-1})/n}$, where the
parameter $\xi_{n}$ depends on the type of $\widehat{\Omega}_{n}$ as follows:

  * A.

$\widehat{\Omega}_{n}$ is of type A if there exists $C_{\text{pre}}>0$ so that
$\widehat{\Omega}_{n}$ satisfies
$\|\widehat{\Omega}_{n}-{\Sigma}_{y}^{-1}\|_{\infty,1}\leq
C_{\text{pre}}s({\Sigma}_{y}^{-1})\sqrt{(\log q)/n}$ with $\mathbb{P}$
probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$. Here we remind the
readers that
$s({\Sigma}_{y}^{-1})=\max_{j\in[q]}\|({\Sigma}_{y}^{-1})_{j}\|_{0}$. In this
case, $\xi_{n}=C_{\text{pre}}\sqrt{s({\Sigma}_{y}^{-1})}$.

  * B.

$\widehat{\Omega}_{n}$ is of type B if
$\|\widehat{\Omega}_{n}-{\Sigma}_{y}^{-1}\|_{\infty,2}\leq
C_{\text{pre}}\sqrt{s({\Sigma}_{y}^{-1})\log(q)/n}$ with $\mathbb{P}$
probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$ for some
$C_{\text{pre}}>0$. In this case, $\xi_{n}=C_{\text{pre}}\max\\{\sqrt{r(\log
q)/n},1\\}$.

  * C.

$\widehat{\Omega}_{n}$ is of type C if
$\widehat{\Omega}_{n}={\Sigma}_{y}^{-1}$. In this case, $\xi_{n}=1$.

The estimation error of $\widehat{\Omega}_{n}$ clearly decays from type A to
C, with the error being zero at type C. Because $\sqrt{r(\log q)/n}$ is
generally much smaller than $s({\Sigma}_{y}^{-1})$, $\xi_{n}$ shrinks from
Case A to Case C monotonously as well. Thus it is fair to say that $\xi_{n}$
reflects the precision of the estimator $\widehat{\Omega}_{n}$ in that
$\xi_{n}$ is smaller if $\widehat{\Omega}_{n}$ is a sharper estimator. We are
now ready to state Theorem 2. This theorem is proved in Appendix B.

###### Theorem 2

Suppose $\log(p\vee q)=o(n)$ and the estimators $\widehat{u}_{n,i}$’s satisfy
Condition 1. Further suppose $\widehat{\Omega}_{n}$ is of type A, B, or C,
which are stated above. Let
$\epsilon_{n}=\xi_{n}\sqrt{\log(p+q)s({\Sigma}_{y}^{-1})/n}$ where $\xi_{n}$
depends on the type of $\widehat{\Omega}_{n}$ as outlined above. Then there
exists a constant $C^{\prime}_{\mathcal{B}}>0$, depending only on
$\mathcal{B}>0$, so that if

$\texttt{Sig}_{y}>2C_{\mathcal{B}}^{\prime}\epsilon_{n},$ (6)

and
$\texttt{cut}\in[C_{\mathcal{B}}^{\prime}\epsilon_{n}/(2\mathcal{B}),{(\theta_{n}-1)}C_{\mathcal{B}}^{\prime}\epsilon_{n}/(2\mathcal{B})]$
with $\theta_{n}=\texttt{Sig}_{y}/(C_{\mathcal{B}}^{\prime}\epsilon_{n})$,
then the algorithm RecoverSupp fully recovers $D(V)$ with $\mathbb{P}$
probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ (for
$\widehat{\Omega}_{n}$ of type A and C), or uniformly over
$\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$ (for
$\widehat{\Omega}_{n}$ of type B).

The assumption that $\log p$ and $\log q$ are $o(n)$ appears in all
theoretical works of CCA (Gao et al., 2017; Mai and Zhang, 2019). A
requirement of this type is generally unavoidable. Note that Theorem 2 implies
a more precise estimator $\widehat{\Omega}_{n}$ requires smaller signal
strength for full support recovery. Before going into subtler implications of
Theorem 2, we make two important remarks.

###### Remark 3

Although the estimation of the high dimensional precision matrix
${\Sigma}_{y}^{-1}$ is potentially complicated, it is often unavoidable owing
to the inherent subtlety of the CCA framework due to the presence of high
dimensional nuisance parameters ${\Sigma}_{x}$ and ${\Sigma}_{y}$. Chen et al.
(2013) also used precision matrix estimator for partial recovery of the
support. In case of sparse CCA, to the best of our knowledge, there does not
exist an algorithm which can recover the support, partially or completely,
without estimating the precision matrix. However, our requirements on
$\widehat{\Omega}_{n}$ are not strict in that many common precision matrix
estimators, e.g. the nodewise Lasso (Theorem 2.4, van de Geer et al., 2014),
the thresholding estimator (cf. Theorem 1 and Section 2.3, Bickel and Levina,
2008), and the CLIME estimator (Theorem 6, Cai et al., 2011) exhibit the decay
rate of type A and B under standard sparsity assumptions on
${\Sigma}_{y}^{-1}$. We will not get into the detail of the sparsity
requirements on ${\Sigma}_{y}^{-1}$ because they are unrelated to the sparsity
of $U$ or $V$, and hence is irrelevant to the primary goal of the current
paper.

###### Remark 4

In the easy regime $s_{y}\lesssim\sqrt{n/(\log(p+q)}$, estimators satisfying
Condition 1 are already available, e.g. COLAR (cf. Theorem 4.2, Gao et al.,
2017) or SCCA (cf. Condition C4 Mai and Zhang, 2019). Thus it is easily seen
that polynomial time support recovery is possible in the easy regime provided
(6) is satisfied. That $r=O(n/(\log(p+q))$ and $s({\Sigma}_{y}^{-1})=O(1)$ are
sufficient conditions for the latter in this regime.

The implications of Theorem 2 in context of the sparsity requirements on
$D(U)$ and $D(V)$ for full support recovery are somewhat implicit through the
assumptions and conditions. However, the restriction on the sparsity is
indirectly imposed by two different sources – which we elaborate on now. To
keep the interpretations simple, throughout the following discussion, we
assume that (a) $r=O(n/\log q)$, (b) $p$ and $q$ are of the same order, and
(c) $s_{x}$ and $s_{y}$ are also of the same order. Since we separate the task
of estimating the nuisance parameter ${\Sigma}_{y}^{-1}$ from the support
recovery of $V$, we also assume that $s({\Sigma}_{y}^{-1})=O(1)$, which
reduces the minimal signal strength condition (18) to $\texttt{Sig}_{y}\geq
C_{\mathcal{B}}\sqrt{\log(p+q)/n}$.

In lieu of the discussion above, the first source of sparsity restriction is
the minimal signal strength condition on $\texttt{Sig}_{y}$ mentioned above.
It is easily seen that $\texttt{Sig}_{y}\leq s_{y}^{-1/2}$. Therefore,
implicit in Theorem 2 lies the condition

$\displaystyle s_{y}\leq\frac{C_{\mathcal{B}}^{2}n}{\log(p+q)}.$ (7)

Thus Theorem 2 does not hold for $s_{y}\gg\log(p+q)/n$ even when
$s({\Sigma}_{y}^{-1})$ and $r$ are small. This regime requires some attention
because in case of sparse PCA (Amini and Wainwright, 2009) and linear
regression (Wainwright, 2009), support recovery at $s\gg\log(p-s)/n$ 444here
and later we will use $s$ to generically denote the sparsity of relevant
parameter vectors in parallel problems like Sparse PCA or Sparse Linear
Regression. is proven to be information theoretically impossible. However,
although a parallel result can be intuited to hold for CCA, the detailed
explorations of the nuances of SCCA support recovery in this regime is yet to
be explored. Therefore, the sparsity requirement in (7) raises the question
whether support recovery for CCA is at all possible when $s_{y}\gg
n/\log(p+q)$, even if ${\Sigma}_{x}$ and ${\Sigma}_{y}$ is known.

###### Question 1

Does there exist any decoder $\widehat{D}$ such that
$\sup_{\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}(\widehat{D}(V)\neq
D(V))\to 0$ when $s_{y}\gg\log(q-s_{y})/n$?

A related question is whether the minimal signal strength requirement (6) is
necessary. To the best of our knowledge, there is no formal study on the
information theoretic limit of the minimal signal strength even in context of
the sparse PCA support recovery. Indeed, as we noted before, the detailed
analyses of support recover for SPCA provided in Amini and Wainwright (2009);
Deshpande and Montanari (2014); Krauthgamer et al. (2015) all is based on the
non-zero principal component elements being $\pm 1/\sqrt{\rm sparsity}$.
Finally, although this question is not directly related with the sparsity
conditions, it indeed probes the sharpness of the results in Theorem 2.

###### Question 2

What is the minimum signal strength required for the recovery of $D(V)$?

We will discuss Question 1 and Question 2 at greater length in Section 3.2. In
particular, Theorem 6(A) shows that there exists $C>0$ so that support
recovery at $s_{y}\geq C\mathcal{B}^{-2}n/\log(q-s_{y})$ is indeed information
theoretically intractable. On the other hand, in Theorem 6(B), we show that
the minimal signal strength has to be of the order
$\mathcal{B}\sqrt{\log(q-s_{y})/n}$ for full recovery of $D(V)$. Thus when
$p\asymp q$ , (6) is indeed necessary from information theoretic perspectives.

The second source of restriction 6on the sparsity lies in Condition 1.
Condition 1 is a $l_{2}$-consistency condition, which has sparsity requirement
itself owing the inherent hardness in the estimation of $U$. Indeed, Theorem
3.3 of Gao et al. (2017) entails that it is impossible to estimate the
canonical directions $u_{i}$’s consistently if $s_{x}>Cn/(r+\log(ep/s_{x}))$
for some large $C>0$. Hence, Condition 1 indirectly imposes the restriction
$s_{x}\lesssim n/\max\\{\log(p/s_{x}),r\\}$. However, when $s_{x}\asymp
s_{y}$, $p\asymp q$, and $r=O(1)$, the above restriction is already absorbed
into the condition $s_{y}\lesssim\mathcal{B}^{-2}n/\log(q-s_{y})$ elicited in
the last paragraph. In fact, there exist consistent estimators of $U$ whenever
$s_{x}\lesssim n/\max\\{\log(p/s_{x}),r\\}$ and $s_{y}\lesssim
n/\max\\{\log(q/s_{y}),r\\}$ (see Gao et al. (2015) or Section 3 of Gao et al.
(2017)). Therefore, in the latter regime, RecoverSupp coupled with the above-
mentioned estimators succeeds. In view of the above, it might be tempting to
think that Condition 1 does not impose significant additional restrictions.
The restriction due to Condition 1, however, is rather subtle and manifests
itself through computational challenges. Note that when support recovery is
information theoretically possible, the computational hardness of recovery by
RecoverSupp will be at least as much as that of the estimation of $U$. Indeed,
the estimators of $U$ which work in the regime $s_{x}\asymp n/\log(p/s_{x})$,
$s_{y}\asymp n/\log(q/s_{y})$ are not adaptive of the sparsity, and they
require a search over exponentially many sets of size $s_{x}$ and $s_{y}$.
Furthermore, under $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$, all polynomial
time consistent estimators of $U$ in the literature, e.g. COLAR (cf. Theorem
4.2, Gao et al., 2017) or SCCA (cf. Condition C4 Mai and Zhang, 2019), require
$s_{x}$, $s_{y}$ to be of the order $\sqrt{n/\log(p+q)}$. In fact, Gao et al.
(2017) indicates that estimation of $U$ or $V$ for much larger sparsity will
be NP hard.

The above raises the question whether RecoverSupp (or any method as such) can
succeed at polynomial time when $\sqrt{n/\log(p+q)}\ll s_{x},s_{y}\lesssim
n/\log(p+q)$. We turn to the landscape of sparse PCA for intuition. Indeed, in
case of sparse PCA, different scenarios are observed in the regime $s\lesssim
n/\log p$ depending on whether $\sqrt{n}\ll s\lesssim n/\log p$, or
$s\lesssim\sqrt{n}$ (we recall that for SPCA we denote the sparsity of the
leading principal component direction generically through $s$). We focus on
the sub-regime $\sqrt{n}\ll s\lesssim n/\log p$ first. In this case, both
estimation and support recovery for sparse PCA are conjectured to be NP hard,
which means no polynomial time method succeeds; see Section 3.3 for more
details. The above hints that the regime $s_{x},s_{y}\gg\sqrt{n}$ is NP hard
for sparse CCA as well.

###### Question 3

Is there any polynomial time method which can recover the support $D(V)$ when
$s_{x},s_{y}\gg\sqrt{n}$?

We dedicate Section 3.3 for answering this question. Subject to the recent
advances in the low degree polynomial conjecture, we establish computational
hardness of the regime $s_{x},s_{y}\gg\sqrt{n}$ (up to a logarithmic factor
gap) subject to $n\lesssim p,q$. Our results are consistent with Gao et al.
(2017)’s findings in the estimation case and covers a broader regime; see
Remark 12 for a comparison.

When the sparsity is of the order $\sqrt{n}$ and $p\asymp n$, however,
polynomial time support recovery and estimation is possible for the sparse PCA
case. Deshpande and Montanari (2014) showed that a coordinate thresholding
type spectral algorithm works in this regime. Thus the following question is
immediate.

###### Question 4

Is there any polynomial time method which can recover the support $D(V)$ when
$s_{x},s_{y}\in[\sqrt{n/\log(p+q)},\sqrt{n}]$?

We give affirmative answer to Question 4 in Section 3.4, which is parallel
with the observations for the sparse PCA. In fact, Corollary 17 shows that
when ${\Sigma}_{x}$ and ${\Sigma}_{y}$ are known, $p+q\asymp n$, and
$s_{x},s_{y}\lesssim\sqrt{n}$, estimation is possible in polynomial time.
Since estimation is possible, RecoverSupp suffices for polynomial time support
recovery in this regime, where $\sqrt{n}$ is well below the information
theoretic limit of $n/\log(p+q)$. The main tool used in Section 3.4 is
coordinate thresholding, which is originally a method for high dimensional
matrix estimation (Bickel and Levina, 2008), and apparently has nothing to do
with estimation of canonical directions. However, under our set up, if the
covariance matrix is consistently estimated in operator norm, by Wedin’s Sin
$\theta$ Theorem (Yu et al., 2015), an SVD is enough to get a consistent
estimator of $U$ and $V$ suitable for further precise analysis.

###### Remark 5

RecoverSupp uses sample splitting, which can reduce the efficiency. One can
swap between the samples and compute two estimators of the supports. One can
easily show that both the intersection and the union of the resulting supports
enjoy the asymptotic garuantees of Theorem 2.

This section can be best summarized by Figure 1, which gives the information
theoretic and computational landscape of sparse CCA analysis in terms of the
sparsity. It can be seen that our contributions (colored in red) complete the
picture, which was initiated by Gao et al. (2017).

COLAR $\sqrt{n/\log(p+q)}$$\sqrt{n}$CTNP
hardEstimation$n/\log(p+q)$Intractable(IT)RecoverSupp $+$
COLAR/SCCA$\sqrt{n/\log(p+q)}$$\sqrt{n}$$s_{x}$RecoverSupp $+$CTRecoverSupp NP
hardSupportrecovery$n/\log(p+q)$Intractable(IT) Figure 1: State of the art for
estimation informetion theoretic and computational limits in sparse CCA. We
have taken $s_{x}=s_{y}$ here. COLAR corresponds to the estimation method of
Gao et al. (2017). Our contributions are colored in red. See Gao et al. (2017)
for more details on the regions colored in blue.

###  3.2 Information Theoretic Lower Bounds: Answers to Question 1 and 2

Theorem 6 establishes the information theoretic limits on the sparsity levels
$s_{x}$, $s_{y}$, and the signal strengths $\texttt{Sig}_{x}$ and
$\texttt{Sig}_{y}$. The proof of Theorem 6 is deferred to Appendix C.

###### Theorem 6

Suppose $\widehat{D}(U)$ and $\widehat{D}(V)$ are estimators of $D(U)$ and
$D(V)$, respectively. Let $s_{x}$, $s_{y}>1$, and $p-s_{x},q-s_{y}>16$. Then
the following assertions hold:

  * A.

If $s_{x}>16n/\\{(\mathcal{B}^{2}-1)\log(p-s_{x})\\}$, then

$\inf_{\widehat{D}}\sup_{\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\bigg{(}\widehat{D}(U)\neq
D(U)\bigg{)}>1/2.$

On the other hand, if $s_{y}>16n/\\{(\mathcal{B}^{2}-1)\log(q-s_{y})\\}$, then

$\inf_{\widehat{D}}\sup_{\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\bigg{(}\widehat{D}(V)\neq
D(V)\bigg{)}>1/2.$

  * B.

Let $\mathcal{P}_{\texttt{Sig}}(r,s_{x},s_{y},\mathcal{B})$ to be the class of
distributions $\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ satisfying
$\texttt{Sig}^{2}_{x}\leq(\mathcal{B}^{2}-1)(\log(p-s_{x}))/(8n)$. Then

$\inf_{\widehat{D}}\sup_{\mathbb{P}\in\mathcal{P}_{\texttt{Sig}}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\bigg{(}\widehat{D}(U)\neq
D(U)\bigg{)}>1/2.$

On he other hand, if
$\texttt{Sig}^{2}_{y}\leq(\mathcal{B}^{2}-1)(\log(q-s_{y}))/(8n)$, then

$\inf_{\widehat{D}}\sup_{\mathbb{P}\in\mathcal{P}_{\texttt{Sig}}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\bigg{(}\widehat{D}(V)\neq
D(V)\bigg{)}>1/2.$

In both cases, the infimum is over all possible decoders $\widehat{D}(U)$ and
$\widehat{D}(V)$.

First we discuss the implications of part A of Theorem 6. This part entails
that for full support recovery of $V$, the minimum sample size requirement is
of the order $s_{y}\log(q-s_{y})$. This requirement is consistent with the
traditional lower bound on $n$ in context of support recovery for sparse PCA
(Amini and Wainwright, 2009, Theorem 3) and $L_{1}$ regression (cf. Corollary
1 of Wainwright, 2009). However, when $r=O(1)$, the sample size requirement
for estimation of $V$ is slightly relaxed, that is, $n\gg s_{y}\log(q/s_{y})$
(cf. Theorem 3.2, Gao et al., 2017). Therefore, from information theoretic
point of view, the task of full support recovery appears to be slightly harder
than the task of estimation. The scenario for partial support recovery might
be different and we do not pursue it here. Moreover, as mentioned earlier, in
the regime $s_{y}\lesssim C_{\mathcal{B}}n/\log(p+q)$, RecoverSupp works with
Gao et al. (2017)’s (see Section 3 therein) estimator of $U$. Thus part A of
Theorem 6 implies that $n/\log(p+q)$ is the information theoretic upper bound
on the sparsity for the full support recovery of sparse CCA.

Part B of Theorem 6 implies that it is not possible to push the minimum signal
strength below the level $O(\sqrt{\log(q-s_{y})/n})$. Thus the minimal signal
strength requirement (6) by Theorem 2 is indeed minimal up to a factor of
$\xi_{n}\sqrt{s({\Sigma}_{y}^{-1})}$. The last statement can be refined
further. To that end, we remind the readers that for a good estimator of
${\Sigma}_{y}^{-1}$, i.e. a type B estimator, $\xi_{n}=O(1)$ if $r=O(n/\log
q)$. However, the latter always holds if support recovery is at all possible,
because in that case $s_{y}\lesssim n/\log(p+q)$, and elementary linear
algebra gives $s_{y}\geq r$. Thus, it is fair to say that, provided a good
estimator of ${\Sigma}_{y}^{-1}$, the requirement (6) is minimal up to a
factor of $\sqrt{s({\Sigma}_{y}^{-1})}$. Indeed, this implies that for banded
inverses with finite band-width our results are rate optimal.

It is further worth comparing this part of the result to the SPCA literature.
In the SPCA support recovery literature, generally, the lower bound on the
signal strength is depicted in terms of the sparsity $s$, and usually a signal
strength of order $O(1/\sqrt{s})$ is postulated (Deshpande and Montanari,
2014; Amini and Wainwright, 2009; Krauthgamer et al., 2015). Using our proof
strategies, it can be easily shown that for SPCA, the analogous lower bound on
the signal strength would be $\sqrt{\log(p-s)/n}$. The latter is generally
much smaller than $1/\sqrt{s}$ and only when $s\asymp n/\log(p)$, the
requirement of $1/\sqrt{s}$ is close to the lower bound. Thus, in the regime
$s\lesssim\sqrt{n/\log p}$, clearly the actual lower bound should be of the
order $O(1/s)$. Therefore the signal strength requirement of $O(1/\sqrt{s})$
typically assumed in literature seems much larger than necessary even for SPCA
context.

###  3.3 Computational Limits and Low Degree Polynomials: Answer to Question 3

We have so far explored the information theoretic upper and lower bounds for
recovering the true support of leading canonical correlation directions.
However, as indicated in the discussion preceding Question 3, the
statistically optimal procedures in the regime where $\sqrt{n}\lesssim
s_{x},s_{y}\lesssim n/\log{(p+q)}$ are computationally intensive and is of
exponential complexity (as a function of $p,q$) in this regime.  In
particular, Gao et al. (2017) have already showed that when $s_{x}$ and
$s_{y}$ belong to parts of this regime, estimation of the canonical correlates
is computationally hard, subject to a computational complexity based “Planted
Clique Conjecture”. For the case of support recovery, the SPCA has been
explored in detail and the corresponding computational hardness has been
established in analogous regimes – see e.g. Amini and Wainwright (2009);
Deshpande and Montanari (2014); Krauthgamer et al. (2015) for details. A
similar phenomenon of computational hardness is observed in case of SPCA spike
detection problem (Berthet and Rigollet, 2013). In light of the above, it is
natural to believe that the SCCA support recovery is also computationally hard
in the regime $\sqrt{n}\lesssim s_{x},s_{y}\lesssim n/\log{(p+q)}$ and as a
result yields a statistical-computational gap. Although several paths exist to
provide evidence towards such gaps 555e.g. Planted Clique Conjecture (Berthet
and Rigollet, 2013; Gao et al., 2017; Brennan et al., 2018), Statistical Query
based lower bounds (Kearns, 1998; Feldman and Kanade, 2012; Brennan et al.,
2020; Dudeja and Hsu, 2021), and Overlap Gap Property based analysis (Gamarnik
and Zadik, 2017; Gamarnik et al., 2019; Arous et al., 2020) ., the recent
developments using “Predictions from Low Degree Polynomials” (Hopkins and
Steurer, 2017; Hopkins, 2018; Kunisky et al., 2019) is particularly appealing
due its simplicity in exposition. In order to show computationally hardness of
the SCCA support recovery problem in the $s\in(\sqrt{n},n/\log(p+q))$ regime,
we shall resort to this very style of ideas, which has so far been applied
successfully to explore statistical-computational gaps under sparse PCA (Ding
et al., 2019), Stochastic Block Models, and tensor PCA (Hopkins, 2018), among
others. This will allow us to explore the computational hardness of the
problem in the entire regime where

$s_{x}+s_{y}\gtrsim(\sqrt{n})(\log n)^{c},$ (8)

compared to the somewhat partial results (see Remark 12 for detailed
comparison) in earlier literature.

We divide our discussions to argue the existence of a statistical-
computational gap in this regime as follows. Starting with a brief background
on the statistical literature on such gaps we first present a natural
reduction of our problem to a suitable hypothesis testing problem in Section
3.3.1. Subsequently, in Section 3.3.2 we present the main idea of the “low
degree polynomial conjecture” by appealing to the recent developments in
Hopkins and Steurer (2017); Hopkins (2018); Kunisky et al. (2019). Finally, we
present our main result for this regime in Section 3.3.3, thereby providing
evidence of the aforementioned gap modulo the Low Degree Polynomial Conjecture
presented in Conjecture 8.

####  3.3.1 Reduction to Testing Problem:

Denote by $\mathbb{Q}$ the distribution of a $N_{p+q}(0,I_{p+q})$ random
vector. Therefore $(X,Y)\sim\mathbb{Q}$ corresponds to the case when $X$ and
$Y$ are uncorrelated. We first show that there is any scope of support
recovery in $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ only if
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ is distinguishable from $\mathbb{Q}$,
i.e. the test $H_{0}:(X,Y)\sim\mathbb{Q}$ vs
$H_{1}:(X,Y)\sim\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ has
asymptotic zero error.

To formalize the ideas, suppose we observe i.i.d random vectors
$\\{X_{i},Y_{i}\\}_{i=1}^{n}$ which are distributed either as $\mathbb{P}$ or
$\mathbb{Q}$. We denote the n-fold product measures corresponding to
$\mathbb{P}$ and $\mathbb{Q}$ by $\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$,
respectively. Note that if
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$, then
$\mathbb{P}_{n}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$. We overload
notation, and denote the combined sample $\\{X_{i}\\}_{i=1}^{n}$ and
$\\{Y_{i}\\}_{i=1}^{n}$ by $\mathbf{X}$ and $\mathbf{Y}$ respectively. In this
section, $\mathbf{X}$ and $\mathbf{Y}$ should be viewed as unordered sets. The
test $\Phi_{n}:\mathbb{R}^{pn+qn}\mapsto\\{0,1\\}$ for testing the null
$H_{0}:(\mathbf{X},\mathbf{Y})\sim\mathbb{Q}_{n}$ vs the alternative
$H_{1}:(\mathbf{X},\mathbf{Y})\sim\mathbb{P}_{n}$ is said to strongly
distinguish $\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$ if

$\lim_{n}\mathbb{Q}_{n}(\Phi_{n}(\mathbf{X},\mathbf{Y})=1)+\lim_{n}\mathbb{P}_{n}(\Phi_{n}(\mathbf{X},\mathbf{Y})=0)=0.$

The above implies that both the type I error and the type II error of
$\Phi_{n}$ converges to zero as $n\to\infty$. In case of composite alternative
$H_{1}:(\mathbf{X},\mathbf{Y})\sim\mathbb{P}_{n}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$,
the test strongly distinguishes $\mathbb{Q}_{n}$ from
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$ if

$\liminf_{n\to\infty}\bigg{\\{}\mathbb{Q}_{n}(\Phi_{n}(\mathbf{X},\mathbf{Y})=1)+\sup_{\mathbb{P}_{n}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}}\mathbb{P}_{n}(\Phi_{n}(\mathbf{X},\mathbf{Y})=0)\bigg{\\}}=0.$

Now we explain how support recovery and the testing framework are connected.
Suppose there exist decoders which exactly recover $D(U)$ and $D(V)$ under
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ for $\mathcal{B}\geq 0$. Then the
trivial test, which rejects the null if either of the estimated supports is
non-empty, strongly distinguishes $\mathbb{Q}_{n}$ from
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$. The above can be coined as the
following lemma.

###### Lemma 7

Suppose there exist polynomial time decoders $\hat{D}_{x}$ and $\hat{D}_{y}$
of $D(U)$ and $D(V)$ so that

$\liminf_{n\to\infty}\sup_{\mathbb{P}_{n}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}}\mathbb{P}_{n}\bigg{(}\hat{D}_{x}(\mathbf{X},\mathbf{Y})=D(U)\text{
and }\hat{D}_{y}(\mathbf{X},\mathbf{Y})=D(V)\bigg{)}=1$ (9)

Further assume,
$\mathbb{Q}_{n}(\hat{D}_{x}(\mathbf{X},\mathbf{Y})=\emptyset)\to 1$, and
$\mathbb{Q}_{n}(\hat{D}_{y}(\mathbf{X},\mathbf{Y})=\emptyset)\to 1$. Then
there exists a polynomial time test which strongly distinguishes
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$ and $\mathbb{Q}_{n}$.

Thus, if a regime does not allow any polynomial time test for distinguishing
$\mathbb{Q}_{n}$ from $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$, there can
be no polynomial time computable consistent decoder for $D(U)$ and $D(V)$.
Therefore, it suffices to show that there is no polynomial time test which
distinguishes $\mathbb{Q}_{n}$ from
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}$ in the regime
$s_{x},s_{y}\gg\sqrt{n}$. To be more explicit, we want to show that if
$s_{x},s_{y}\gg\sqrt{n}$, then

$\displaystyle\liminf_{n\to\infty}\bigg{\\{}\mathbb{Q}_{n}(\Phi_{n}(\mathbf{X},\mathbf{Y})=1)+\sup_{\mathbb{P}_{n}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})^{n}}\mathbb{P}_{n}(\Phi_{n}(\mathbf{X},\mathbf{Y})=0)\bigg{\\}}>0$
(10)

for any $\Phi_{n}$ that is computable in polynomial time.

The testing problem under concern is commonly known as the CCA detection
problem, owing to its alternative formulation as $H_{0}:\Lambda_{1}=0$ vs
$H_{1}:\Lambda_{1}>0$. In other words, the test tries to detect if there is
any any signal in the data. Note that, Lemma 7 also implies that detection is
an easier problem than support recovery in that the former is always possible
whenever the latter is feasible. The opposite direction may not be true,
however, since detection does not reveal much information on the support.

####  3.3.2 Background on the Low-degree Framework:

We shall provide a brief introduction to the low-degree polynomial conjecture
which forms the basis of our analyses here, and refer the interested reader to
Hopkins and Steurer (2017); Hopkins (2018); Kunisky et al. (2019) for in-depth
discussions on the topic. We will apply this method in context of the test
$H_{0}:(\mathbf{X},\mathbf{Y})\sim\mathbb{Q}_{n}$ vs
$H_{1}:(\mathbf{X},\mathbf{Y})\sim\mathbb{P}_{n}$. The low-degree method
centres around the likelihood ratio $\mathbb{L}_{n}$ which takes the form
$\frac{d\mathbb{P}_{n}}{d\mathbb{Q}_{n}}$ in the above framework. Our key tool
here will be the Hermite polynomials, which form a basis system of
$L_{2}(\mathbb{Q}_{n})$ (Szegö, 1939). Central to the low-degree approach lies
the projection of $\mathbb{L}_{n}$ onto the subspace (of
$L_{2}(\mathbb{Q}_{n})$) formed by the Hermite polynomials of degree at most
$D_{n}\in\mathbb{N}$. The latter projection, to be denoted by
$\mathbb{L}_{n}^{\leq D_{n}}$ from now on, is important because it measures
how well polynomials of degree $\leq D_{n}$ can distinguish $\mathbb{P}_{n}$
from $\mathbb{Q}_{n}$. In particular,

$\|\mathbb{L}_{n}^{\leq D_{n}}\|_{L_{2}(\mathbb{Q}_{n})}:=\max_{f\text{ deg
}{\leq
D_{n}}}\frac{\mathbb{E}_{\mathbb{P}_{n}}[f(\mathbf{X},\mathbf{Y})]}{\sqrt{\mathbb{E}_{\mathbb{Q}_{n}}[f(\mathbf{X},\mathbf{Y})^{2}]}},$
(11)

where the maximization is over polynomials
$f:\mathbb{R}^{n(p+q)}\mapsto\mathbb{R}$ of degree at most $D_{n}$ (Ding et
al., 2019).

The $L_{2}(\mathbb{Q}_{n})$ norm of the un-truncated likelihood ration
$\mathbb{L}_{n}$ has long held an important place in the theory hypothesis
testing since $\|\mathbb{L}_{n}\|_{L_{2}(\mathbb{Q}_{n})}=O(1)$ implies
$\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$ are asymptotically indistinguishable.
While the un-truncated likelihood ratio $\mathbb{L}_{n}$ is connected to the
existence of _any_ distinguishing test, degree $D_{n}$ projections of
$\mathbb{L}_{n}$ are connected to the existence of polynomial time
distinguishing tests. The implications of the above heuristics are made
precise by the following conjecture (cf. Hypothesis 2.1.5 of Hopkins, 2018).

######  Conjecture 8 (Informal)

Suppose $t:\mathbb{N}\mapsto\mathbb{N}$. For “nice” sequences of distributions
$\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$, if $\|\mathbb{L}_{n}^{\leq
D_{n}}\|_{L_{2}(\mathbb{Q}_{n})}=O(1)$ as $n\to\infty$ whenever $D_{n}\leq
t(n)\text{polylog}(n)$, then there is no time-$n^{t(n)}$ test
$\Phi_{n}:\mathbb{R}^{n(p+q)}\mapsto\\{0,1\\}$ that strongly distinguishes
$\mathbb{P}_{n}$ and $\mathbb{Q}_{n}$.

Thus Conjecture 8 implies that the degree-$D_{n}$ polynomial
$\mathbb{L}_{n}^{\leq D_{n}}$ is a proxy for time-$n^{t(n)}$ algorithms
(Kunisky et al., 2019). If we can show that $\|\mathbb{L}_{n}^{\leq
D_{n}}\|_{L_{2}(\mathbb{Q}_{n})}=O(1)$ for a $D_{n}$ of the order $(\log
n)^{1+\epsilon}$ for some $\epsilon>0$, then the low degree Conjecture says
that no polynomial time test can strongly distinguish $\mathbb{P}_{n}$ and
$\mathbb{Q}_{n}$ (Conjecture 1.16 of Kunisky et al., 2019).

Conjecture 8 is informal in the sense that we do not specify the “nice”
distributions, which are defined in Section 4.2.4 of Kunisky et al. (2019)
(see also Conjecture 2.2.4 of Hopkins, 2018). Niceness requires
$\mathbb{P}_{n}$ to be sufficiently symmetric, which is generally guaranteed
by naturally occuring high dimensional problems like ours. The condition of
“niceness” is attributed to eliminate pathological cases where the testing can
be made easier by methods like Gaussian elimination. See Hopkins (2018) for
more details.

####  3.3.3 Main Result

Similar to Ding et al. (2019), we will consider a Bayesian framework. It might
not be immediately clear how a Bayesian formulation will fit into the low-
degree framework, and lead to (10). However, the connection will be clear
soon. We put independent Rademacher priors $\pi_{x}$ and $\pi_{y}$ on $\alpha$
and $\beta$. We say $\alpha\sim\pi_{x}$ if $\alpha_{1},\ldots,\alpha_{p}$ are
i.i.d., and for each $i\in[p]$,

$\displaystyle\alpha_{i}=\begin{cases}1/\sqrt{s_{x}}&w.p.\quad s_{x}/(2p)\\\
-1/\sqrt{s_{x}}&w.p.\quad s_{x}/(2p)\\\ 0&w.p.\quad 1-s_{x}/p.\end{cases}$
(12)

The Rademacher prior $\pi_{y}$ can be defined similarly. We will denote the
product measure $\pi_{x}\times\pi_{y}$ by $\pi$. Let us define

$\displaystyle\Sigma(\alpha,\beta,\rho)=\begin{bmatrix}I_{p}&\rho\alpha\beta^{T}\\\
\rho\beta\alpha^{T}&I_{q}\end{bmatrix},\quad\alpha\in\mathbb{R}^{p},\
\beta\in\mathbb{R}^{q},\ \rho>0.$ (13)

When $\rho\|\alpha\|_{2}\|\beta\|_{2}<1$, $\Sigma(\alpha,\beta,\rho)$ is the
covariance matrix corresponding to $X\sim N_{p}(0,I_{p})$ and $Y\sim
N_{q}(0,I_{q})$ with covariance $\text{cov}(X,Y)=\rho\alpha\beta^{T}$. Hence,
for $\Sigma(\alpha,\beta,\rho)$ to be positive definite,
$\|\alpha\|_{2}\|\beta\|_{2}<1/\rho$ is a sufficient condition. The priors
$\pi_{x}$ and $\pi_{y}$ put positive weight on $\alpha$ and $\beta$ that do
not lead to a positive definite $\Sigma(\alpha,\beta,\rho)$, and hence calls
for extra care during the low-degree analysis. This subtlety is absent in the
sparse PCA analogue (Ding et al., 2019).

Let us define

$\mathbb{P}_{\alpha,\beta}=\begin{cases}N(0,\Sigma(\alpha,\beta,1/\mathcal{B}))&\text{when
}\|\alpha\|_{2}\|\beta\|_{2}<\mathcal{B}\\\ \mathbb{Q}&\text{o.w.}\end{cases}$
(14)

We denote the $n$-fold product measure corresponding to
$\mathbb{P}_{\alpha,\beta}$ by $\mathbb{P}_{n,\alpha,\beta}$. If
$(\mathbf{X},\mathbf{Y})\mid\alpha,\beta\sim\mathbb{P}_{n,\alpha,\beta}$, then
the marginal density of $(\mathbf{X},\mathbf{Y})$ is
$\mathbb{E}_{\alpha\sim\pi_{x},\beta\sim\pi_{y}}d\mathbb{P}_{n,\alpha,\beta}$.
The following lemma, which is proved in Appendix G.3, explains how the
Bayesian framework is connected to (10).

###### Lemma 9

Suppose $\mathcal{B}>2$ and $s_{x},s_{y}\to\infty$. Then

$\liminf_{n}\sup_{\mathbb{P}_{n}\in\mathcal{P}_{G}(r,2s_{x},2s_{y},\mathcal{B})^{n}}\mathbb{P}_{n}\Big{(}\Phi_{n}(\mathbf{X},\mathbf{Y})=0\Big{)}\geq\liminf_{n}\mathbb{E}_{\pi}\mathbb{P}_{n,\alpha,\beta}\Big{(}\Phi_{n}(\mathbf{X},\mathbf{Y})=0\Big{)},$

where $\mathbb{E}_{\pi}$ is the shorthand for
$\mathbb{E}_{\alpha\sim\pi_{x},\beta\sim\pi_{y}}$.

Note that a similar result holds for $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$
as well because
$\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})\subset\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$.
Lemma 9 implies that to show (10), it suffices to show that a polynomial time
computable $\Phi_{n}$ fails to strongly distinguish the marginal distribution
of $\mathbf{X}$ and $\mathbf{Y}$ from $\mathbb{Q}_{n}$. However, the latter
falls within the realms of the low degree framework. To see this, note that
the corresponding likelihood ratio takes the form

$\mathbb{L}_{n}=\frac{\mathbb{E}_{\alpha\sim\pi_{x},\beta\sim\pi_{y}}d\mathbb{P}_{n,\alpha,\beta,\mathcal{B}}}{d\mathbb{Q}_{n}(\mathbf{X},\mathbf{Y})}.$
(15)

If we can show that $\|\mathbb{L}_{n}^{\leq
D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}=O(1)$ for some $D_{n}=O(\log n)$, then
Conjecture 8 would indicate that a $n^{\tilde{\Theta}(D_{n})}$-time computable
$\Phi_{n}$ fails to distinguish the distribution of
$\mathbb{E}_{\alpha\sim\pi_{x},\beta\sim\pi_{y}}d\mathbb{P}_{n,\alpha,\beta,\mathcal{B}}$
from $\mathbb{Q}_{n}$. Theorem 10 accomplishes the above under some additional
conditions on $p$, $q$, and $n$, which we will discuss shortly. Theorem 10 is
proved in Appendix D.

###### Theorem 10

Suppose $D_{n}\leq\min(\sqrt{p},\sqrt{q},n)$,

$s_{x},s_{y}\geq\sqrt{enD_{n}}/\mathcal{B}\quad\text{and}\quad p,q\geq
3en/\mathcal{B}^{2}.$ (16)

Then $\|\mathbb{L}_{n}^{\leq D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}$ is $O(1)$
where $\mathbb{L}_{n}$ is as defined in (15).

The following Corollary results from combining Lemma 9 with Theorem 10.

###### Corollary 11

Suppose

$s_{x},s_{y}\geq 2\sqrt{enD_{n}}/\mathcal{B}\quad\text{and}\quad p,q\geq
3en/\mathcal{B}^{2}.$ (17)

If Conjecture 8 is true, then for $D_{n}\leq\min(\sqrt{p},\sqrt{q},n)$, there
is no time-$n^{\tilde{\Theta}(D_{n})}$ test that strongly distinguishes
$\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$ and $\mathbb{Q}_{n}$.

Corollary 11 conjectures that polynomial time algorithms can not strongly
distinguish $\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})^{n}$ and
$\mathbb{Q}_{n}$ provided $s_{x},s_{y}$, $p$, and $q$ satisfy (17). Therefore
under (17), Lemma 7 conjectures support recovery to be NP hard.

Now we discuss a bit on condition (17). The first constraint in (17) is
expected because it ensures $s_{x},s_{y}\gg\sqrt{n}$, which indicates that the
sparsity is in the hard regime. We need to explain a bit on why the other
constraint $p,q>3en/\mathcal{B}^{2}$ is needed. If $n\gg p,q$, the sample
canonical correlations are consistent, and therefore strong separation is
possible in polynomial time without any restriction on the sparsity (Bao et
al., 2019; Ma and Yang, 2021). Even if $p/n\to c_{1}\in(0,1)$ and $q/n\to
c_{2}\in(0,1)$, then also strong separation is possible in model 13 provided
the canonical correlation $\rho$ is larger than some threshold depending on
$c_{1}$ and $c_{2}$ (Bao et al., 2019). The restriction
$p,q>3en/\mathcal{B}^{2}$ ensures that the problem is hard enough so that the
vanilla CCA does not lead to successful detection. The constant $3e$ is not
sharp and possibly can be improved. The necessity of the condition $p,q\gtrsim
n/\mathcal{B}^{2}$ is unknown for support recovery, however. Since support
recovery is a harder problem than detection, in the hard regime, polynomial
time support recovery algorithms may fail at a weaker condition on $n$, $p$,
and $q$.

###### Remark 12

[Comparison with previous work:] As mentioned earlier, Gao et al. (2017) was
the first to discover the existence of computational gap in context of sparse
CCA. In their seminal work, Gao et al. (2017) established the computational
hardness of CCA estimation problem at a particular subregime of
$s_{x},s_{y}\gg\sqrt{n}/(\mathcal{B}\sqrt{\log(p+q)})$ provided
$\mathcal{B}\to\infty$ is allowed. In view of the above, it was hinted that
sparse CCA becomes computationally hard when
$s_{x},s_{y}\gg\sqrt{n}/(\mathcal{B}\sqrt{\log(p+q)})$. However, when
$\mathcal{B}$ is bounded, the entire regime
$s_{x},s_{y}\gg\sqrt{n}/(\mathcal{B}\sqrt{\log(p+q)})$ is probably not
computationally hard. In Section 3.4, we show that if $p+q\asymp n$, then both
polynomial time estimation and support recovery are possible even if
$s_{x}+s_{y}\lesssim\sqrt{n}$, at least in the known ${\Sigma}_{x}$ and
${\Sigma}_{y}$ case. The latter sparsity regime can be considerably larger
than $s_{x},s_{y}\lesssim\sqrt{n/\log(p+q)}$. Together, Section 3.4 and the
current section indicate that in the bounded $\mathcal{B}$ case, the
transition of computational hardness for sparse CCA probably happens at the
sparsity level $\sqrt{n}$, not $\sqrt{n/\log(p+q)}$, which is consistent with
sparse PCA. Also, the low-degree polynomial conjecture allowed us to explore
almost the entire targeted regime $s_{x},s_{y}\gg\sqrt{n}$, where Gao et al.
(2017), who used the planted clique conjecture, considers only a subregime of
$s_{x},s_{y}\gg\sqrt{n}/(\mathcal{B}\sqrt{\log(p+q)})$.

###  3.4 A Polynomial Time Algorithm for $\sqrt{n}/\log{(p+q)}\ll
s_{x},s_{y}\ll\sqrt{n}$ Regime : Answer to Question 4

In this subsection, we show that in the difficult regime
$s_{x}+s_{y}\in[\sqrt{n/\log(p+q)},\sqrt{n}]$, using a soft coordinate
thresholding (CT) type algorithm, we can estimate the canonical directions
consistently for our purpose when $p+q\asymp n$. CT was introduced by the the
seminal work of Bickel and Levina (2008) for the purpose of estimating high
dimensional covariance matrices. For SPCA, Deshpande and Montanari (2014)’s CT
is the only algorithm which provably recovers the full support in the
difficult regime (Krauthgamer et al., 2015). In context of CCA, Chen et al.
(2013) uses CT for partial support recovery in the rank one model under what
we referred to as the easy regime. However, Chen et al. (2013)’s main goal was
the estimation of the leading canonical vectors, not support recovery. As a
result, Chen et al. (2013) detects the support of the relatively large
elements of the leading canonical directions, which are subsequently used to
obtain consistent preliminary estimators of the leading canonical directions.
Our thresholding level and theoretical analysis is different from that of Chen
et al. (2013) because the analytical tools used in the easy regime does not
work in the difficult regime.

####  3.4.1 Methodology: Estimation via CT

By thresholding a matrix $A$ coordinate-wise, we will roughly mean the process
of assigning the value zero to any element of $A$ which is below a certain
threshold in absolute value. Similar to Deshpande and Montanari (2014), we
will consider the soft thresholding operator, which, at threshold level $t$,
takes the form

$\eta(x,t)=\begin{cases}x-t&x>t\\\ 0&|x|<t\\\ x+t&x<-t.\end{cases}$

It will be worth noting that the soft thresholding operator
$x\mapsto\eta(x,t)$ is continuous.

Algorithm 2 Coordinate Thresholding (CT) for CCA

0:

  1. 1.

Sample covariance matrices $\widehat{\Sigma}_{n,xy}^{(1)}$ and
$\widehat{\Sigma}_{n,xy}^{(2)}$ based on samples
$O_{1}=(x_{i},y_{i})_{i=1}^{[n/2]}$ and $O_{2}=(x_{i},y_{i})_{i=[n/2]+1}^{n}$,
respectively.

  2. 2.

Variances ${\Sigma}_{x}$ and ${\Sigma}_{y}$.

  3. 3.

Parameters Thr and cut.

  4. 4.

$r$, i.e. rank of ${\Sigma}_{xy}$

0: $\widehat{D}(V)$.

  1. 1.

Peeling: calculate
$\tilde{\Sigma}_{xy}={\Sigma}_{x}^{-1}\widehat{\Sigma}_{n,xy}^{(1)}{\Sigma}_{y}^{-1}$.

  2. 2.

Threshold: Letting $N=m+n$, perform soft thresholding
$x\mapsto\eta(x;\texttt{Thr}/\sqrt{N})$ entrywise on $\tilde{\Sigma}_{xy}$ to
obtain thresholded $\eta(\tilde{\Sigma}_{xy})$.

  3. 3.

Sandwitch:
$\eta(\tilde{\Sigma}_{xy})\mapsto{\Sigma}_{x}^{1/2}\eta(\tilde{\Sigma}_{xy}){\Sigma}_{y}^{1/2}$.

  4. 4.

SVD:  Find $\widehat{U}_{pre}$, the matrix of the leading $r$ singular vector
of ${\Sigma}_{x}^{1/2}\eta(\tilde{\Sigma}_{xy}){\Sigma}_{y}^{1/2}$.

  5. 5.

Premultiply: Set $\widehat{U}^{(1)}={\Sigma}_{x}^{-1/2}\widehat{U}_{pre}$.

Return: RecoverSupp
$(\widehat{U}^{(1)},\texttt{cut},{\Sigma}_{y}^{-1},\widehat{\Sigma}_{n,xy}^{(2)},r)$
where RecoverSupp is given by Algorithm 1.

We will also assume that the covariance matrices ${\Sigma}_{x}$ and
${\Sigma}_{y}$ are known. To understand the difficulty of unknown
${\Sigma}_{x}$ and ${\Sigma}_{y}$, we remind the readers that
${\Sigma}_{xy}={\Sigma}_{x}U\Lambda V^{T}{\Sigma}_{y}$. Because the matrices
$U$ and $V$ are sandwiched by the matrices ${\Sigma}_{x}$ and ${\Sigma}_{y}$,
their sparsity pattern does not get reflected in the sparsity pattern of
${\Sigma}_{xy}$. Therefore, if one blindly applies CT to
$\widehat{\Sigma}_{n,xy}$, they can at best hope to recover the sparsity
pattern of the outer matrices ${\Sigma}_{x}$ and ${\Sigma}_{y}$. If the
supports of the matrices $U$ and $V$ are of main concern, CT should rather be
applied on the matrix
$\tilde{\Sigma}_{xy}={\Sigma}_{x}^{-1}\widehat{\Sigma}_{n,xy}{\Sigma}_{y}^{-1}$.
If ${\Sigma}_{x}$ and ${\Sigma}_{y}$ are unknown, one needs to efficiently
estimate $\tilde{\Sigma}_{xy}$ before the application of CT. Although under
certain structural conditions, it is possible to find rate optimal estimators
$\widehat{\Sigma}_{n,x}^{-1}$ and $\widehat{\Sigma}_{n,y}^{-1}$ of
${\Sigma}_{x}^{-1}$ and ${\Sigma}_{y}^{-1}$ at least in theory, the errors
$\|(\widehat{\Sigma}_{n,x}^{-1}-{\Sigma}_{x}^{-1})\widehat{\Sigma}_{n,xy}{\Sigma}_{y}^{-1}\|_{op}$
and
$\|{\Sigma}_{x}^{-1}\widehat{\Sigma}_{n,xy}(\widehat{\Sigma}_{n,y}^{-1}-{\Sigma}_{y}^{-1})\|_{op}$
may still blow up due to the presence of the high dimensional matrix
$\widehat{\Sigma}_{n,xy}$, which can be as big as $O(\sqrt{(p+q)/n})$ in
operator norm. One may be tempted to replace $\widehat{\Sigma}_{n,xy}$ by a
sparse estimator of ${\Sigma}_{xy}$ to facilitate faster estimation, but that
does not work because we explicitly require the formulation of
$\widehat{\Sigma}_{n,xy}$ as the sum of Wishart matrices (see equation 28 in
the proof). The latter representation, which is critical for the sharp
analysis, may not be preserved by a CLIME (Cai et al., 2011) or nodewise Lasso
estimator (van de Geer et al., 2014) of ${\Sigma}_{xy}$.

We remark in passing that it is possible to obtain an estimator $\widehat{A}$
so that $|\widehat{A}-\tilde{\Sigma}_{xy}|_{\infty}=o_{p}(1)$. Although the
latter does not provide much control over the operator norm of
$\widehat{A}-\tilde{\Sigma}_{xy}$, it is sufficient for partial support
recovery, e.g. the recovery of the rows of $U$ or $V$ with strongest signals.
(See Appendix B of Chen et al., 2013, for example, for some results in this
direction under the easy regime when $r=1$.)

As indicated by the previous paragraph, we apply coordinate thresholding to
the matrix
$\tilde{\Sigma}_{xy}={\Sigma}_{x}^{-1}\widehat{\Sigma}_{n,xy}{\Sigma}_{y}^{-1}$,
which directly targets the matrix
${\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}=U\Lambda V^{T}$. We call this
step the peeling step because it extracts the matrix $\tilde{\Sigma}_{xy}$
from the sandwiched matrix
$\widehat{\Sigma}_{n,xy}={\Sigma}_{x}\tilde{\Sigma}_{xy}{\Sigma}_{y}$. We then
perform the entry-wise coordinate thresholding algorithm on the peeled form
$\tilde{\Sigma}_{xy}$ with threshold Thr so as to obtain
$\eta(\tilde{\Sigma}_{xy};\texttt{Thr}/\sqrt{n})$. We postpone the discussion
on Thr to Section 3.4.2. The thresholded matrix is a good estimator of
${\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}$, but we need an estimator of
${\Sigma}_{x}^{-1/2}{\Sigma}_{xy}{\Sigma}_{y}^{-1/2}$. Therefore, we again
sandwich $\tilde{\Sigma}_{xy}$ between ${\Sigma}_{x}^{1/2}$ and
${\Sigma}_{y}^{1/2}$. The motivation behind this sandwiching is that if
$\|\tilde{\Sigma}_{xy}-{\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}\|_{op}=\epsilon_{n}$,
then ${\Sigma}_{x}^{1/2}\tilde{\Sigma}_{xy}{\Sigma}_{y}^{1/2}$ is a good
estimator of ${\Sigma}_{x}^{-1/2}{\Sigma}_{xy}{\Sigma}_{y}^{-1/2}$ in that

$\|{\Sigma}_{x}^{1/2}\tilde{\Sigma}_{xy}{\Sigma}_{y}^{1/2}-{\Sigma}_{x}^{-1/2}{\Sigma}_{xy}{\Sigma}_{y}^{-1/2}\|_{op}\leq\sqrt{\|{\Sigma}_{x}\|_{op}\|{\Sigma}_{y}\|_{op}}\epsilon_{n}\leq\mathcal{B}\epsilon_{n}.$

However, ${\Sigma}_{x}^{1/2}U\Lambda V^{T}{\Sigma}_{y}^{1/2}$ is an SVD of
${\Sigma}_{x}^{-1/2}{\Sigma}_{xy}{\Sigma}_{y}^{-1/2}$. Therefore, using Davis-
Kahan sin theta theorem (Yu et al., 2015), one can easily see that the SVD of
${\Sigma}_{x}^{1/2}\tilde{\Sigma}_{xy}{\Sigma}_{y}^{1/2}$ produces estimators
$\widehat{U}^{\prime}$ and $\widehat{V}^{\prime}$ whose columns are
$\epsilon_{n}$-consistent in $l_{2}$ norm for the columns of
${\Sigma}_{x}^{1/2}U$ and ${\Sigma}_{y}^{1/2}V$ up to a sign flip (cf. Theorem
2, Yu et al., 2015). Pre-multiplying the resulting $U^{\prime}$ by
${\Sigma}_{x}^{-1/2}$ yields an estimator $\widehat{U}$ of $U$ up to a sign
flip of the columns. We do not worry about the sign flip because Condition 1
allows for the sign flips of the columns. Therefore, we feed this
$\widehat{U}$ into RecoverSupp as our final step. See Algorithm 2 for more
details.

###### Remark 13

In case of electronics health records data, it is possible to obtain large
surrogate data on $X$ and $Y$ separately and thus might allow relaxing the
known precision matrices assumption above. We do not pursue such semi-
supervised setups here.

####  3.4.2 Analysis of the CT Algorithm

For the asymptotic analysis of the CT algorithm, we will assume the underlying
distribution to be Gaussian, i.e.
$\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$. This Gaussian
assumption will be used to perform a crucial decomposition of sample
covariance matrix which typically holds for Gaussian random vectors (see
equation 28). Deshpande and Montanari (2014), who used similar devices for
obtaining the sharp rate results in SPCA, also required a similar Gaussian
assumption. We do not yet know how to extend these results to sub-Gaussian
random vectors.

Let us consider the threshold $\texttt{Thr}/\sqrt{n}$, where Thr is explicitly
given in Theorem 14. Unfortunately, tuning of Thr requires the knowledge of
the underlying sparsity $s_{x}$ and $s_{y}$. Similar to Deshpande and
Montanari (2014), our thresholding level is different than the traditional
choice of order $O(\sqrt{\log(p+q)/n})$ in the easy regime (Bickel and Levina,
2008; Cai et al., 2012; Chen et al., 2013). The latter level is too large to
successfully recover all the non-zero elements in the difficult regime. We
threshold $\tilde{\Sigma}_{xy}$ at a lower level at higher sparsity, which in
its turn, complicates the analysis to a greater degree. Our main result in
this direction, stated in Theorem 14, is proved in Appendix E.

###### Theorem 14

Suppose
$(X_{i},Y_{i})\sim\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$.
Further suppose $s_{x}+s_{y}<\sqrt{n}$, $p\vee q=o(\log n)$, and $\log
n=o(\sqrt{p}\vee\sqrt{q})$. Let $K$ and $C_{1}$ be constants so that $K\geq
1288\mathcal{B}^{4}$ and $C_{1}\geq C\mathcal{B}^{4}$, where $C>0$ is an
absolute constant. Suppose the threshold level Thr is defined by

$\displaystyle\texttt{Thr}=\begin{cases}\sqrt{C_{1}\log(p+q)}&\text{ if
}\quad(s_{x}+s_{y})^{2}<2^{1/4}(p+q)^{3/4}\text{ (case i)}\\\
\Big{(}K\log(\frac{p+q}{(s_{x}+s_{y})^{2}})\Big{)}^{1/2}&\text{ if }\quad
2^{1/4}(p+q)^{3/4}\leq(s_{x}+s_{y})^{2}\leq(p+q)/e\text{ (case ii)}\\\
0&\text{ o.w. (case iii).}\end{cases}$

Suppose $c_{\mathcal{B}}$ is a constant that takes the value $K$, $C_{1}$, or
one in case (i), (ii), and (iii), respectively. Then there exists an absolute
constant $C>0$ so that the following holds with probability $1-o(1)$ for
$\tilde{\Sigma}_{xy}={\Sigma}_{x}^{-1}\widehat{\Sigma}_{n,xy}{\Sigma}_{y}^{-1}:$

$\|\eta(\tilde{\Sigma}_{xy};\eta)-{\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}\|_{op}\leq
C\mathcal{B}^{2}\frac{(s_{x}+s_{y})}{\sqrt{n}}\max\bigg{\\{}\bigg{(}c_{\mathcal{B}}\log(\frac{p+q}{(s_{x}+s_{y})^{2}})\bigg{)}^{1/2},1\bigg{\\}}.$

To disentangle the implications of the theorem above, let us assume $p+q\asymp
n$ for the time being. Then case (ii) in the theorem corresponds to
$n^{3/4}\lesssim(s_{x}+s_{y})^{2}\leq n$. Thus, CT works in the difficult
regime provided $p+q\asymp n$. It should be noted that the threshold for this
case is almost of the order $O(1/\sqrt{n})$, which is much smaller than
$O(\sqrt{\log(p+q)/n})$, the traditional threshold for the easy regime. Next,
observe that case (i) is an easy case because $s_{x}+s_{y}$ is much smaller
than $\sqrt{n}$. Therefore, in this case, the traditional threshold of the
easy regime works. Case (iii) includes the hard regime, where polynomial time
support recovery is probably impossible. Because it is unlikely that CT can
improve over the vanilla estimator $\tilde{\Sigma}_{xy}$ in this regime, a
threshold of zero is set.

###### Remark 15

Theorem 14 requires $\log n=o(\sqrt{p}\vee\sqrt{q})$ because one of our
concentration inequalities in the analysis of case (ii) needs this technical
condition (see Lemma 27). The omitted regime $\log n>C(\sqrt{p}\vee\sqrt{q})$
is indeed an easier one, where special methods like CT is not even required.
In fact, it is well known that subgaussian $X$ and $Y$ satisfy (cf. Theorem
4.7.1 of of Vershynin, 2018)

$\|\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy}\|_{op}\leq
C\bigg{(}\bigg{(}\frac{p+q}{n}\bigg{)}^{1/2}+\frac{p+q}{n}\bigg{)},$

which is $O(\log n/\sqrt{n})$ in the regime under concern. We decided include
this result in the statement of Theorem 14 since it would unnecessarily
lengthen the exposition. Therefore, in this section, we exclude this regime
from our consideration to focus more on the $s_{x}+s_{y}\approx\sqrt{p+q}$
regime.

###### Remark 16

The statement of Theorem 14 is not explicit on the lower bound of the constant
$C_{1}$. However, our simulation shows setting $C_{1}\geq 50\mathcal{B}^{4}$
works. Both threshold parameters $C_{1}$ and $K$ in Theorem 14 depend on the
unknown $\mathcal{B}>0$. The proof actually shows that $\mathcal{B}$ can be
replaced by
$\max\\{\Lambda_{max}({\Sigma}_{x}),\Lambda_{max}({\Sigma}_{y}),\Lambda_{max}({\Sigma}_{x}^{-1}),\Lambda_{max}({\Sigma}_{y}^{-1})\\}$.
In practice, estimating the latter also can be difficult. Therefore our
suggestion is to set $K$ and $C_{1}$ to be some large numbers, or to use cross
validation to choose them.

Finally, Theorem 14 leads to the following corollary, which establishes that
in the difficult regime, there exist estimators which satisfy Condition 1, and
Algorithm 2 succeeds with probability one provided $p+q\asymp n$. This answers
Question 4 in the affirmative for Gaussian distributions.

###### Corollary 17

Instate the conditions of Theorem 14. Then there exists $C_{\mathcal{B}}>0$ so
that if

$n\geq
C_{\mathcal{B}}r(s_{x}+s_{y})^{2}\max\bigg{\\{}\log\bigg{(}\frac{p+q}{(s_{x}+s_{y})^{2}}\bigg{)},1\bigg{\\}},$
(18)

then the $\widehat{U}^{(1)}$ defined in Algorithm 2 satisfies Condition 1, and
$\inf_{\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}($Algorithm
2 correctly recovers $D(V)$ $)\to_{n}1$.

We defer the proof of Corollary 17 to Appendix F.

##  4 Numerical Experiments

This section illustrates the performance of different polynomial time CCA
support recovery methods when the sparsity transitions from the easy to
difficult regime. We base our demonstration on a Gaussian rank one model, i.e.
$(X,Y)$ are jointly Gaussian with covariance matrix
${\Sigma}_{xy}=\rho{\Sigma}_{x}\alpha\beta^{T}{\Sigma}_{y}$. For simplicity,
we take $p=q$ and $s_{x}=s_{y}=s$. In all our simulations, $\rho$ is set to be
$0.5$, and $\alpha=\alpha^{*}/\sqrt{(\alpha^{*})^{T}{\Sigma}_{x}\alpha^{*}}$,
$\beta=\beta^{*}/\sqrt{(\beta^{*})^{T}{\Sigma}_{y}\beta^{*}}$ where

$\displaystyle\alpha^{*}=$ $\displaystyle\
(1/\sqrt{s},\ldots,1/\sqrt{s},0,\ldots,0),$ $\displaystyle\beta^{*}=$
$\displaystyle\
\Big{(}\sqrt{1-(s-1)s^{-4/3}},s^{-2/3},\dots,s^{-2/3},0,\ldots,0\Big{)}$

are unit norm vectors. Note that the order of most elements of $\beta$ is
$O(s^{-2/3})$, where a typical element of $\alpha$ is $O(s^{-1/2})$.
Therefore, we will refer to $\alpha$ and $\beta$ as the moderate and the small
signal case, respectively. For the population covariance matrices
${\Sigma}_{x}$ and ${\Sigma}_{y}$ of $X$ and $Y$, we consider the following
two scenarios:

A (Identity):

${\Sigma}_{x}=I_{p}$ and ${\Sigma}_{y}=I_{q}$. Since $p=q$, they are
essentially the same.

B (Sparse inverse):

This example is taken from Gao et al. (2017). In this case,
${\Sigma}_{x}^{-1}={\Sigma}_{y}^{-1}$ are banded matrices, whose entries are
given by

$({\Sigma}_{x}^{-1})_{i,j}=1\\{i=j\\}+0.65\times 1\\{|i-j|=1\\}+0.4\times
1\\{|i-j|=2\\}.$

Now we explain our common simulation scheme. We take the sample size $n$ to be
$1000$, and consider three values for $p$: $100$, $200$, and $300$. The
highest value of $p+q$ is thus $600$, which is smaller than but in proportion
to $n$ regime. Our simulations indicate that all of the methods considered
here requires $n$ to be quite larger than $p+q$ for the asymptotics to kick in
at $\rho=0.5$ and we comment on it as required below. We further let
$s/\sqrt{n}$ vary in the set $[0.01,2]$. To be more specific, we consider $16$
equidistant points in the set $[0.01,2]$ for the ratio $s/\sqrt{n}$.

Now we discuss the error metric used here to compare the performance of
different support recovery methods. Type I and type II errors are commonly
used tools to measure the performance of support recovery (Deshpande and
Montanari, 2014). In case of support recovery of $\alpha$, we define the type
I error to be the proportion of zero elements in $\alpha$ that appear in the
estimated support $\widehat{D}(\alpha)$. Thus, we quantify the type I error of
$\alpha$ by $|\widehat{D}(\alpha)\setminus D(\alpha)|/(p-s)$. On the other
hand, the type II error for $\alpha$ is the proportion of elements in
$D(\alpha)$ which are absent in $\widehat{D}(\alpha)$, i.e. the type II error
is quantified by $|D(\alpha)\setminus\widehat{D}(\alpha)|/s$. One can define
the type I and type II errors corresponding to $\beta$ similarly. Our
simulations demonstrate that often the methods with low type I error exhibit
high type II error, and vice versa. In such situations, comparison between the
corresponding methods become difficult if one uses the type I and type II
errors separately. Therefore, we consider a scaled Hamming loss type metric
which suitably combines the type I and type II error. The symmetric Hamming
error of estimating $D(\alpha)$ by $\widehat{D}(\alpha)$ is (cf. Section 2.1
of Wang et al., 2021)

$\bigg{(}1-\frac{|D(\alpha)\cap\widehat{D}(\alpha)|}{\sqrt{|D(\alpha)||\widehat{D}(\alpha)|}}\bigg{)}.$

Note that the above quantity is always bounded above by one. We can similarly
define the symmetric Hamming distance between $D(\beta)$ and
$\widehat{D}(\beta)$. Finally, the estimates of these three errors (Type I,
Type, and scaled Hamming Loss) are obtained based on $1000$ Monte Carlo
replications.

Now we discuss the support recovery methods we compare here.

Naïve SCCA

We estimate $\alpha$ and $\beta$ using the SCCA method of Mai and Zhang
(2019), and set $\widehat{D}(\alpha)=\\{i\in[p]:\widehat{\alpha}_{i}\neq 0\\}$
and $\widehat{D}(\beta)=\\{i\in[q]:\widehat{\beta}_{i}\neq 0\\}$, where
$\widehat{\alpha}$ and $\widehat{\beta}$ are the corresponding SCCA
estimators. To implement the SCCA method of Mai and Zhang (2019), we use the R
code referred therein with default tuning parameters.

Cleaned SCCA

This method implements RecoverSupp with the above mentioned SCCA estimators of
$\alpha$ and $\beta$ as the preliminary estimators.

CT

This is the method outlined in Algorithm 2, which is RecoverSupp coupled with
the CT estimators of $\alpha$ and $\beta$.

Our CT method requires the knowledge of the population covariance matrices
${\Sigma}_{x}$ and ${\Sigma}_{y}$. Therefore, to keep the comparison fair, in
case of the cleaned SCCA method as well, we implement RecoverSupp with the
popular covariance matrices. Because of their reliance on RecoverSupp, both
cleaned SCCA and CT depend on the threshold cut, tuning which seems to be a
non-trivial task. We set
$\texttt{cut}=C\sqrt{\log(p+q)s({\Sigma}_{x}^{-1})/n}$, where $C$ is the
thresholding constant. Our simulations show that a large $C$ results in high
type II error, where insufficient thresholding inflates the type I error.
Taking the hamming loss into account, we observe that $C\approx 1$ leads to a
better performance in case A in an overall sense. On the other hand, case B
requires a smaller value of thresholding parameter. In particular, we let $C$
to be one in case A, and set $C=0.05$ and $0.2$, respectively, for the support
recovery of $\alpha$ and $\beta$ in case B. The CT algorithm requires an extra
threshold parameter, namely the parameter Thr in Algorithm 2, which
corresponds to the coordinate thresholding step. We set Thr in accordance with
Theorem 14 and Remark 16, with $K$ being $1288\mathcal{B}^{4}$ and $C_{1}$
being $50\mathcal{B}^{4}$. We set $\mathcal{B}$ as in Remark 16, that is

$\mathcal{B}=\max\\{\Lambda_{max}({\Sigma}_{x}),\Lambda_{max}({\Sigma}_{y}),\Lambda_{max}({\Sigma}_{x}^{-1}),\Lambda_{max}({\Sigma}_{y}^{-1})\\}.$

The errors incurred by our methods in case A are displayed in Figure 2 (for
$\alpha$) and Figure 3 (for $\beta$). Figures 4 and 5, on the other hand,
display the errors in the recovery of $\alpha$ and $\beta$, respectively, in
case B.

Now we discuss the main observations from the above plots. When the sparsity
parameter $s$ is considerably low (less than ten in the current settings), the
naïve SCCA method is sufficient in the sense that the specialized methods do
not perform any better. Moreover, the naïve method is the most conservative
one among all three methods. As a consequence, the associated type I error is
always small, although the type II error of the naïve method grows faster than
any other method. The specialized methods are able to improve the type II
error at the cost of higher type I error. At a higher sparsity level, the
specialized methods can outperform the naïve method in terms of the Hamming
error, however. This is most evident when the setting is also complex, i.e.
the signal is small, or the underlying covariance matrices are not identity.
In particular, Figure 2 and 4 entail that when the signal strength is moderate
and the sparsity is high, the cleaned SCCA has the lowest hamming error. In
the small signal case, however, CT exhibits the best hamming error as
$s/\sqrt{n}$ increases; cf. Figure 3 and 5.

The Type I error of CT can be slightly improved if the sparsity information
can be incorporated during the thresholding step. We simply replace cut by the
maximum of cut and the $s$-th largest element of $\widehat{V}^{clean}$, where
the latter is as in Algorithm RecoverSupp. See for example Figure 6, which
entails that this modification reduces the Hamming error of the CT algorithm
in case A. Our empirical analysis hints that the CT algorithm has potential
for improvement from the implementation perspective. In particular, it may be
desirable to obtain a more efficient procedure for choosing cut in a
systematic way. However, such a detailed numerical analysis is beyond the
scope of the current paper and will require further modifications of the
initial methods for estimation of $\alpha,\beta$ both for scalability and
finite sample performance reasons. We keep these explorations as important
future directions.

It is natural to wonder what is the effect of cleaning via RecoverSupp on
SCCA. As mentioned earlier, during our simulations we observed that a cleaning
step generally improves the type II error of the naïve SCCA, but it also
increases the type I error. In terms of the combined measure, i.e. the Hamming
error, it turns out that cleaning does have an edge at higher sparsity levels
in case B; cf. Figure 4 and Figure 5. However, the scenario is different in
case A. Although Figures 2 and 3 indicate that almost no cleaning occurs at
the set threshold level of one, we saw that cleaning happens at lower
threshold levels. However, the latter does not improve the overall Hamming
error of naïve SCCA. The consequence of cleaning may be different for other
SCCA methods.

To summarize, when the sparsity is low, support recovery using the naïve SCCA
is probably as good as the specialized methods. However, at higher sparsity
level, specialized support recovery methods may be preferable. Consequently,
the precise analysis of the apparently naïve SCCA will indeed be an
interesting future direction.

(a) Type I error for support recovery of $\alpha$

(b) Type II error for support recovery of $\alpha$

(c) Symmetrized Hamming error for support recovery of $\alpha$

Figure 2: Support recovery for $\alpha$ when ${\Sigma}_{x}=I_{p}$ and
${\Sigma}_{y}=I_{q}$.

(a) Type I error for support recovery of $\beta$

(b) Type II error for support recovery of $\beta$

(c) Symmetrized Hamming error for support recovery of $\beta$

Figure 3: Support recovery for $\beta$ when ${\Sigma}_{x}=I_{p}$ and
${\Sigma}_{y}=I_{q}$.

(a) Type I error for support recovery of $\alpha$

(b) Type II error for support recovery of $\alpha$

(c) Symmetrized Hamming error for support recovery of $\alpha$

Figure 4: Support recovery for $\alpha$ when ${\Sigma}_{x}$ and ${\Sigma}_{y}$
are the sparse covariance matrices.

(a) Type I error for support recovery of $\beta$

(b) Type II error for the support recovery of $\beta$

(c) Symmetrized Hamming error for support recovery of $\beta$

Figure 5: Support recovery for $\beta$ when ${\Sigma}_{x}$ and ${\Sigma}_{y}$
are the sparse covariance matrices.

(a) Errors for support recovery of $\alpha$

(b) Errors for support recovery of $\beta$

Figure 6: Support recovery by the CT algorithm when we use the information on
sparsity to improve the type I error. Here ${\Sigma}_{x}$ and ${\Sigma}_{y}$
are $I_{p}$ and $I_{q}$, respectively. To see the decrease in type I error,
compare the errors with that of Figure 2 and Figure 3.

##  5 Discussion

In this paper we have discussed rate optimal behavior of information theoretic
and computational limits of the joint support recovery problem for the sparse
canonical correlation analysis problem. Inspired by recent results in
estimation theory of sparse CCA, flurry of results in sparse PCA, and related
developments based on low-degree polynomial conjecture – we are able to paint
a complete picture of the landscape of support recovery for SCCA. For future
directions, it is worth noting that our results are so far not designed to
recover $D(\beta_{i})$ for individual $i\in[r]$ separately (and hence the term
joint recovery). Although this is also the case for most state of the art in
the case of the sparse PCA problem (results often exist only for the combined
support (Deshpande and Montanari, 2014) or the single spike model where $r=1$
(Wainwright, 2009).), we believe that this is an interesting question for
deeper explorations in the future. Moreover, moving beyond asymptotically
exact recovery of support to more nuanced metrics (e.g. Hamming Loss) will
also require new ideas worth studying. Finally, it remains an interesting
question to pursue whether polynomial time support recovery is possible in the
$\sqrt{n}/\log{(p+q)}\ll s_{x},s_{y}\ll\sqrt{n}$ regime using a CT type idea –
but for unknown yet structured high dimensional nuisance parameters
$\Sigma_{x},\Sigma_{y}$.

Acknowledgments

This work was supported by National Institutes of Health grant P42ES030990.

##  A Proof preliminaries

The Appendix collects the proof of all our theorems and lemmas. This section
introduces some new notations and collects some facts, which are used
repeatedly in our proofs.

###  A.1 New Notations

Since the columns of ${\Sigma}_{x}^{1/2}U$, i.e.
$[{\Sigma}_{x}^{1/2}U_{1},\ldots,{\Sigma}_{x}^{1/2}U_{r}]$ are orthogonal, we
can extend it to an orthogonal basis of $\mathbb{R}^{p}$, which can also be
expressed in the form
$[{\Sigma}_{x}^{1/2}u_{1},\ldots,{\Sigma}_{x}^{1/2}u_{p}]$ since
${\Sigma}_{x}$ is non-singular. Let us denote the matrix
$[u_{1},\ldots,u_{p}]$ by $\tilde{U}$, whose first $r$ columns form the matrix
$U$. Along the same line, we can define $\tilde{V}$, whose first $q$ columns
constitute the matrix $V$.

Suppose $A\in\mathbb{R}^{p\times q}$ is a matrix. We define the projection of
$A$ onto $D\subset[p]\times[q]$ by

$\bigg{(}\mathcal{P}_{D}\\{A\\}\bigg{)}_{i,j}=\begin{cases}A_{i,j}&\text{ if
}(i,j)\in D,\\\ 0&\text{otherwise.}\end{cases}$

Also, for any $S\subset[p]$, we let $A_{S*}$ denote the matrix
$\mathcal{P}_{S\times[q]}\\{A\\}$. Similarly, for $F\subset[q]$, we let
$A_{F}$ be the matrix $\mathcal{P}_{[p]\times F}\\{A\\}$. For
$k\in\mathbb{N}$, we define the norms
$\|A\|_{k,\infty}=\max_{j\in[q]}\|A_{j}\|_{k}$ and
$\|A\|_{\infty,k}=\max_{i\in[q]}\|A_{i}\|_{k}$. We will use the notation
$|A|_{\infty}$ to denote the quantity $\sup_{1\in[p],j\in[q]}|A_{i,j}|$.

The Kullback Leibler (KL) divergence between two probability distributions
$P_{1}$ and $P_{2}$ will be denoted by $KL(P_{1}\mid P_{2})$. For
$x\in\mathbb{R}$, we let $\left\lfloor x\right\rfloor$ denote greatest integer
less than or equal to $x\in\mathbb{R}$.

###  A.2 Facts on $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$

First note that since $v_{i}^{T}{\Sigma}_{y}v_{i}=1$ by (2) for all $i\in[q]$,
we have $\|v_{i}\|_{2}\leq\sqrt{\mathcal{B}}$. Similarly, we can also show
that $\|u_{i}\|_{2}\leq\sqrt{\mathcal{B}}$. Second, we note that
$\|{\Sigma}_{x}^{1/2}U\|_{op}=\|{\Sigma}_{y}^{1/2}V\|_{op}=1$, and

$\displaystyle|{\Sigma}_{yx}|_{\infty}\leq\|{\Sigma}_{yx}\|_{op}=$
$\displaystyle\ \|{\Sigma}_{y}V\Lambda
U^{T}{\Sigma}_{x}\|_{op}\leq\|{\Sigma}_{y}^{1/2}\|_{op}\|{\Sigma}_{y}^{1/2}V\|_{op}\|\Lambda\|_{op}\|{\Sigma}_{x}^{1/2}U\|_{op}\|{\Sigma}_{x}^{1/2}\|_{op}\leq\mathcal{B}$
(19)

because the largest element of $\Lambda$ is not larger than one. Since
$X_{i}$’s and $Y_{i}$’s are Subgaussian, for any random vector $v$ independent
of $\mathbf{X}$ and $\mathbf{Y}$, it follows that (cf. Lemma 7 of Janková and
van de Geer, 2018)

$|(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})v|_{\infty}\leq
C_{\mathcal{B}}\|v\|_{2}\sqrt{\frac{\log(p+q)}{n}}$ (20)

with $\mathbb{P}$ probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$. Also, we can show that
$\Phi_{0}={\Sigma}_{y}^{-1}$ satisfies

$\|(\Phi_{0})_{k}\|_{1,\infty}\leq\sqrt{s({{\Sigma}_{x}})}\|(\Phi_{0})_{k}\|_{2,\infty}\leq\sqrt{s({{\Sigma}_{x}})}\|\Phi_{0}\|_{op}\leq\sqrt{s({{\Sigma}_{x}})}\mathcal{B},$

where Cauchy-Schwarz inequality was used in the first step.

###  A.3 General Technical Facts

###### Fact 18

For two matrices $A\in\mathbb{R}^{m\times n}$ and $B\in\mathbb{R}^{n\times
q}$, we have

$\|AB\|_{F}^{2}\leq\|A\|_{op}^{2}\|B\|_{F}^{2},\quad\|AB\|_{F}^{2}\leq\|A\|_{F}^{2}\|B\|_{op}^{2}$

######  Fact 19 (Lemma 11 of Deshpande and Montanari (2014))

Let $\mathbf{Z}\in\mathbb{R}^{n\times p}$ be a matrix with i.i.d. standard
normal entries, i.e. $Z_{i,j}\sim N(0,1)$. Then for every $t>0$,

$\mathbb{P}(\|\mathbf{Z}\|_{op}\geq\sqrt{p}+\sqrt{n}+t)\leq\exp(-t^{2}/2).$

As a consequence, there exists an absolute constant $C>0$ such that

$\mathbb{P}\Big{(}\|\mathbf{Z}\|_{op}\geq\sqrt{2}(\sqrt{p}+\sqrt{n})\Big{)}\leq\exp(-C(p+n)).$

Recall that for $A\in\mathbb{R}^{p\times q}$, in Appendix A.1, we defined
$\|A\|_{1,\infty}$ and $\|A\|_{\infty,1}$ to be the matrix norms
$\max_{j\in[q]}\|A_{j}\|_{1}$ and $\max_{i\in[p]}\|A_{i*}\|_{1}$,
respectively.

The following fact is a Corollary to (20).

###### Fact 20

Suppose $X$ and $Y$ are jointly subgaussian. Then
$|\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy}|_{\infty}=O_{p}(\sqrt{\log(p+q)/n})$.

######  Fact 21 (Chi-square tail bound)

Suppose $\mathbb{Z}_{1},\ldots,\mathbb{Z}_{k}\stackrel{{\scriptstyle
iid}}{{\sim}}N(0,1)$. Then for any $y>5$, we have

$\mathbb{P}\Big{(}\sum_{l=1}^{k}\mathbb{Z}_{l}^{2}\geq
yk\Big{)}\leq\exp(-yk/5).$

Proof [Proof of Fact 21] Since $Z_{l}$’s are independent standard Gaussian
random variables, by tail bounds on Chi-squared random variables (The form
below is from Lemma 12 of Deshpande and Montanari, 2014),

$\mathbb{P}\Big{(}\sum_{l=1}^{k}\mathbb{Z}_{l}^{2}\geq
k+2\sqrt{kx}+2x\Big{)}\leq\exp(-x).$

Plugging in $x=yk$, we obtain that

$\mathbb{P}\Big{(}\sum_{l=1}^{k}\mathbb{Z}_{l}^{2}\geq(1+2\sqrt{y}+2y)k\Big{)}\leq\exp(-yk),$

which implies for $y>1$,

$\mathbb{P}\Big{(}\sum_{l=1}^{k}\mathbb{Z}_{l}^{2}\geq
5yk\Big{)}\leq\exp(-yk),$

which can be rewritten as

$\mathbb{P}\Big{(}\sum_{l=1}^{k}\mathbb{Z}_{l}^{2}\geq
yk\Big{)}\leq\exp(-yk/5)$

as long as $y>5$.

##  B Proof of Theorem 2

For the sake of simplicity, we denote $\widehat{U}^{(1)}$,
$\widehat{\Sigma}_{n,xy}^{(2)}$, and $\widehat{\Omega}_{n}^{(1)}$ by
$\widehat{U}$, $\widehat{\Sigma}_{n,xy}$, and $\widehat{\Omega}_{n}$,
respectively. The reader should keep in mind that $\widehat{U}$ is independent
of $\widehat{\Sigma}_{n,xy}$ and $\widehat{\Omega}_{n}$ because it is
constructed from a different sample. Next, using Condition 1, we can show that
there exists $(w_{i},\ldots,w_{p})\in\\{\pm 1\\}^{p}$ so that

$\inf_{\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}\max_{i\in[r]}\Big{|}(w_{i}\widehat{u}_{n,i}-u_{i})^{T}{\Sigma}_{x}(w_{i}\widehat{u}_{n,i}-u_{i})\Big{|}<\texttt{Err}^{2}\Big{)}\to
1$

as $n\to\infty$. Without loss of generality, we assume $w_{i}=1$ for all
$i\in[r]$. The proof will be similar for general $w_{i}$’s. Thus

$\displaystyle\inf_{\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}\max_{i\in[r]}\Big{|}(\widehat{u}_{n,i}-u_{i})^{T}{\Sigma}_{x}(\widehat{u}_{n,i}-u_{i})\Big{|}<\texttt{Err}^{2}\Big{)}\to
1$ (21)

Therefore $\|\widehat{u}_{n,i}-u_{i}\|_{2}\leq\texttt{Err}\sqrt{\mathcal{B}}$
for all $i\in[r]$ with $\mathbb{P}$ probability tending to one.

Now we will collect some facts which will be used during the proof. Because
$\widehat{u}_{n,i}$ and $\widehat{\Sigma}_{n,yx}$ are independent, (20)
implies that

$|(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}|_{\infty}\leq
C_{\mathcal{B}}\|\widehat{u}_{n,i}\|_{2}\sqrt{\frac{\log(p+q)}{n}}.$

Using (21), we obtain that
$\|\widehat{u}_{n,i}\|_{2}\leq\|\widehat{u}_{n,i}-u_{i}\|_{2}+\|u_{i}\|_{2}\leq\sqrt{B}(\texttt{Err}+1)$.
Because $\texttt{Err}<\mathcal{B}^{-1}\leq 1$, we have

$\inf_{\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\bigg{(}\max_{i\in[r]}|(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}|_{\infty}\leq
C_{\mathcal{B}}\sqrt{\frac{\log(p+q)}{n}}\bigg{)}=1-o(1).$ (22)

Noting (19) implies
$|{\Sigma}_{yx}\widehat{u}_{n,i}|_{\infty}\leq\|{\Sigma}_{yx}\|_{op}\|\widehat{u}_{n,i}\|_{2}\leq
2\mathcal{B}^{3/2}$, and that $\log(p+q)=o(n)$, using (22), we obtain that

$\displaystyle\max_{i\in[r]}|\widehat{\Sigma}_{n,yx}\widehat{u}_{n,i}|_{\infty}\leq|(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}|_{\infty}+|{\Sigma}_{yx}\widehat{u}_{n,i}|_{\infty}\leq
3\mathcal{B}^{3/2}$ (23)

with $\mathbb{P}$ probability $1-o(1)$.

Now we are ready to prove Theorem 2. Because
$\Lambda_{i}(v_{i})_{k}=e_{k}^{T}{\Sigma}_{y}^{-1}{\Sigma}_{yx}u_{i}$, it
holds that

$\displaystyle(\widehat{v}_{n,i}^{clean})_{k}-\Lambda_{i}(v_{i})_{k}=$
$\displaystyle\
e_{k}^{T}(\widehat{\Omega}_{n}-\Phi_{0})\widehat{\Sigma}_{n,yx}\widehat{u}_{n,i}+e_{k}^{T}\Phi_{0}(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}+e_{k}^{T}\Phi_{0}{\Sigma}_{yx}(\widehat{u}_{n,i}-u_{i})$

leading to

$\displaystyle|(\widehat{v}_{n,i}^{clean})_{k}-\Lambda_{i}(v_{i})_{k}|\leq$
$\displaystyle\
\underbrace{|e_{k}^{T}(\widehat{\Omega}_{n}-\Phi_{0})\widehat{\Sigma}_{n,yx}\widehat{u}_{n,i}|}_{T_{1}(i,k)}+\underbrace{|e_{k}^{T}\Phi_{0}(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}|}_{T_{2}(i,k)}+\underbrace{|e_{k}^{T}\Phi_{0}{\Sigma}_{yx}(\widehat{u}_{n,i}-u_{i})|}_{T_{3}(i,k)}.$

Handling the term $T_{2}$ is the easiest because

$\displaystyle\max_{i\in[r],k\in[q]}T_{2}(i,k)\leq\|\Phi_{0}\|_{1,\infty}|(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}|_{\infty}\leq
C_{\mathcal{B}}\sqrt{\frac{s({\Sigma}_{y}^{-1})\log(p+q)}{n}}$

with $\mathbb{P}$ probability $1-o(1)$ uniformly over
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$, where we used (22) and the fact that
$\|\Phi_{0}\|_{1,\infty}\leq\sqrt{s({{\Sigma}_{x}})}\mathcal{B}$. The
difference in cases (A), (B), (C) arises only due to different bounds on
$T_{1}(i,k)$ in these cases. We demonstrate the whole proof only for case (A).
For the other two cases, we only discuss the analysis of $T_{1}(i,k)$ because
the rest of the proof remains identical in these cases.

####  B.0.1 Case (A)

Since we have shown in (23) that
$|\widehat{\Sigma}_{n,yx}\widehat{u}_{n,i}|_{\infty}\leq 3\mathcal{B}^{3/2}$,
we calculate

$\displaystyle\max_{i\in[r],k\in[q]}T_{1}(i,k)\leq\|\widehat{\Omega}_{n}-\Phi_{0}\|_{1,\infty}\max_{i\in[r]}|\widehat{\Sigma}_{n,yx}\widehat{u}_{n,i}|_{\infty}\leq
3\mathcal{B}^{3/2}C_{\text{pre}}s({\Sigma}_{y}^{-1})\sqrt{\frac{\log q}{n}}$

with $\mathbb{P}$ probability tending to one, uniformly over
$\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$, where to get the last inequality, we
also used the bound on $\|\widehat{\Omega}_{n}-\Phi_{0}\|_{\infty,1}$ in case
(A).

Finally, for $T_{3}$, we notice that

$T_{3}(i,k)=\big{|}e_{k}^{T}\Phi_{0}{\Sigma}_{yx}(\widehat{u}_{n,i}-u_{i})\big{|}=\bigg{|}e_{k}^{T}\sum_{j=1}^{r}\Lambda_{j}v_{j}u_{j}^{T}{\Sigma}_{x}(\widehat{u}_{n,i}-u_{i})\bigg{|}\leq\max_{j\in[r]}\big{|}(v_{j})_{k}\big{|}\bigg{|}\sum_{j=1}^{r}u_{j}^{T}{\Sigma}_{x}(\widehat{u}_{n,i}-u_{i})\bigg{|}$

since $\Lambda_{1}\leq 1$. Since $(v_{j})_{k}=V_{kj}$, it is clear that
$T_{3}(i,k)$ is identically zero if $k\notin D(V)$. Otherwise, Cauchy Schwarz
inequality implies,

$\bigg{|}\sum_{j=1}^{r}u_{j}^{T}{\Sigma}_{x}(\widehat{u}_{n,i}-u_{i})\bigg{|}\leq\sqrt{r}\bigg{(}\sum_{j=1}^{r}(u_{j}^{T}{\Sigma}_{x}(\widehat{u}_{n,i}-u_{i}))^{2}\bigg{)}^{1/2}\leq\sqrt{r}\|{\Sigma}_{x}^{1/2}(\widehat{u}_{n,i}-u_{i})\|_{2}$

because ${\Sigma}_{x}^{1/2}u_{j}$’s are orthogonal. Thus

$\displaystyle\max_{i\in[r],k\in
D(V)}|T_{3}(i,k)|\leq\sqrt{r}\max_{j\in[r]}\big{|}(v_{j})_{k}\big{|}\texttt{Err}.$

Now we will combine the above pieces together. Note that

$\displaystyle\max_{i\in[q]}\max_{k\in[r]}(|T_{1}(i,k)|+|T_{2}(i,k)|)\leq
C_{\mathcal{B}}\underbrace{C_{\text{pre}}s({\Sigma}_{y}^{-1})\sqrt{\frac{\log(p+q)}{n}}}_{\epsilon_{n}}.$
(24)

For $k\notin D(V)$, denoting the $i$-th column of $\widehat{V}^{clean}$ by
$\widehat{v}_{n,i}^{clean}$ we observe that,

$\displaystyle\max_{k\notin
D(V)}\max_{i\in[r]}|\widehat{V}_{ki}^{clean}|=\max_{k\notin
D(V)}\max_{i\in[r]}|(\widehat{v}_{n,i}^{clean})_{k}|\leq\max_{i\in[q]}\max_{k\in[r]}(|T_{1}(i,k)|+|T_{2}(i,k)|)\leq
C_{\mathcal{B}}{\epsilon_{n}}$ (25)

with $\mathbb{P}$ probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$. On the other hand, if
$k\in D(v_{i})$, then we have for all $i\in[r]$,

$|(\widehat{v}_{n,i}^{clean})_{k}|>\Lambda_{i}|(v_{i})_{k}|-\sqrt{r}\max_{j\in[r]}\big{|}(v_{j})_{k}\big{|}\texttt{Err}-\max_{i\in[q]}\max_{k\in[r]}(|T_{1}(i,k)|+|T_{2}(i,k)|),$

which implies

$\max_{i\in[r]}|\widehat{V}_{ki}^{clean}|>\max_{i\in[r]}\Lambda_{i}|(v_{i})_{k}|-\sqrt{r}\max_{i\in[r]}\big{|}(v_{i})_{k}\big{|}\texttt{Err}-C_{\mathcal{B}}\epsilon_{n}.$

Since $\texttt{Err}<\mathcal{B}^{-1}/(2\sqrt{r})$ and
$\mathcal{B}^{-1}<\min_{i\in[r]}\Lambda_{i}$, we have

$\max_{i\in[r]}\Lambda_{i}|(v_{i})_{k}|-\sqrt{r}\max_{i\in[r]}\big{|}(v_{i})_{k}\big{|}\texttt{Err}>(\mathcal{B}^{-1}-\sqrt{r}\texttt{Err})\max_{i\in[r]}\big{|}(v_{i})_{k}\big{|}>\mathcal{B}^{-1}\max_{i\in[r]}\big{|}(v_{i})_{k}\big{|}/2.$

Thus, noting $V_{ki}=(v_{i})_{k}$, we obtain that

$\min_{k\in D(V)}\max_{i\in[r]}|(\widehat{v}_{n,i}^{clean})_{k}|=\min_{k\in
D(V)}\max_{i\in[r]}|\widehat{V}_{ki}^{clean}|>\min_{k\in
D(V)}\max_{i\in[r]}\big{|}V_{ki}\big{|}/(2\mathcal{B})-C_{\mathcal{B}}\epsilon_{n}$

with $\mathbb{P}$ probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$. Suppose
$C_{\mathcal{B}}^{\prime}=2\mathcal{B}C_{\mathcal{B}}$. Note that
$\min_{k\in[p]}\max_{i\in[r]}\big{|}(v_{i})_{k}\big{|}=\theta_{n}C_{\mathcal{B}}^{\prime}\epsilon_{n}$
where $\theta_{n}>2$. Then

$\min_{k\in
D(V)}\max_{i\in[r]}\widehat{V}_{ki}^{clean}>(\theta_{n}-1)C_{\mathcal{B}}^{\prime}\epsilon_{n}/(2\mathcal{B}).$

with $\mathbb{P}$ probability $1-o(1)$ uniformly over
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$. This, combined with
(25) implies setting
$\texttt{cut}\in[C_{\mathcal{B}}^{\prime}\epsilon_{n}/(2\mathcal{B}),(\theta_{n}-1)C_{\mathcal{B}}^{\prime}\epsilon_{n}/(2\mathcal{B})]$
leads to full support recovery with $\mathbb{P}$ probability $1-o(1)$. The
proof of the first part follows.

####  B.0.2 Case (B)

In the Gaussian case, we resort to the hidden variable representation of $X$
and $Y$ due to Bach and Jordan (2005), which enables sharper bound on the term
$T_{1}(i,k)$. Suppose $\mathbf{Z}\sim N_{r}(0,I_{r})$ where $r$ is the rank of
${\Sigma}_{xy}$. Consider $Z_{1}\sim N_{p}(0,I_{p})$ and $Z_{2}\sim
N_{q}(0,I_{q})$ independent of $\mathbf{Z}$. Then $X$ and $Y$ can be
represented as

$X=\mathcal{W}_{1}Z+\ \mathcal{H}_{1}Z_{1}\quad\text{ and }\quad
Y=\mathcal{W}_{2}Z+\mathcal{H}_{2}Z_{2},$ (26)

where

$\mathcal{W}_{1}={\Sigma}_{x}U\Lambda^{1/2},\
\mathcal{W}_{2}={\Sigma}_{y}V\Lambda^{1/2},\
\mathcal{H}_{1}=({\Sigma}_{x}-\mathcal{W}_{1}\mathcal{W}_{1}^{T})^{1/2},\ \
\text{and}\
\mathcal{H}_{2}=({\Sigma}_{y}-\mathcal{W}_{2}\mathcal{W}_{2}^{T})^{1/2}.$

Here $({\Sigma}_{x}-\mathcal{W}_{1}\mathcal{W}_{1}^{T})^{1/2}$ is well defined
because
${\Sigma}_{x}-\mathcal{W}_{1}\mathcal{W}_{1}^{T}={\Sigma}_{x}\tilde{U}(I_{p}-\Lambda_{x})\tilde{U}^{T}{\Sigma}_{x}$,
where $\Lambda_{x}$ is a $p\times p$ diagonal matrix whose first $p$ elements
are $\Lambda_{1},\ldots,\Lambda_{r}$, and they rest are zero. Because
$\Lambda_{1}\leq 1$, we have
$({\Sigma}_{x}-\mathcal{W}_{1}\mathcal{W}_{1}^{T})^{1/2}={\Sigma}_{x}\tilde{U}(I_{p}-\Lambda_{x})^{1/2}\tilde{U}^{T}{\Sigma}_{x}$.
Similarly, we can show that
$({\Sigma}_{y}-\mathcal{W}_{2}\mathcal{W}_{2}^{T})^{1/2}={\Sigma}_{y}\tilde{V}(I_{q}-\Lambda_{y})^{1/2}\tilde{V}^{T}{\Sigma}_{y}$
where $\Lambda_{y}$ is the diagonal matrix whose first $r$ elements are
$\Lambda_{1},\ldots,\Lambda_{r}$, and the rest are zero. It can be easily
verified that

$Var(X)=\mathcal{W}_{1}\mathcal{W}_{1}^{T}+\mathcal{H}_{1}={\Sigma}_{x},\quad
Var(Y)=\mathcal{W}_{2}\mathcal{W}_{2}^{T}+\mathcal{H}_{2}={\Sigma}_{y},\quad{and}\quad{\Sigma}_{xy}=\mathcal{W}_{1}\mathcal{W}_{2}^{T}={\Sigma}_{x}U\Lambda
V^{T}{\Sigma}_{y},$

which ensures that the joint variance of $(X,Y)$ is still $\Sigma$. Also, some
linear algebra leads to

$\displaystyle\max\bigg{\\{}\|\mathcal{H}_{1}\|^{2}_{op},\|\mathcal{H}_{2}\|^{2}_{op},\|\mathcal{W}_{1}\|_{op},\|\mathcal{W}_{2}\|_{op}\bigg{\\}}<\mathcal{B}.$
(27)

Suppose we have $n$ independent realizations of the pseudo-observations
$Z_{1}$, $Z_{2}$, and $Z$. Denote by $\mathbf{Z}_{1}$, $\mathbf{Z}_{2}$, and
$\mathbf{Z}$, the stacked data matrices with the i-th row as $(Z_{1})_{i}$,
$(Z_{2})_{i}$, and $Z_{i}$, respectively, where $i\in[n]$. Here we used the
term data-matrix although we do not observe $\mathbf{Z}$, $\mathbf{Z}_{1}$ and
$\mathbf{Z}_{2}$ directly. Due to the representation in (26), the data
matrices $\mathbf{X}$ and $\mathbf{Y}$ have the form

$\mathbf{X}=\mathbf{Z}\mathcal{W}_{1}^{T}+\mathbf{Z}_{1}\mathcal{H}_{1},\quad\mathbf{Y}=\mathbf{Z}\mathcal{W}_{2}^{T}+\mathbf{Z}_{2}\mathcal{H}_{2}.$

We can write the covariance matrix
$\widehat{\Sigma}_{n,xy}=\mathbf{X}^{T}\mathbf{Y}/n$ as

$\widehat{\Sigma}_{n,xy}=\frac{1}{n}\bigg{\\{}\mathcal{W}_{1}\mathbf{Z}^{T}\mathbf{Z}\mathcal{W}_{2}^{T}+\mathcal{W}_{1}\mathbf{Z}^{T}\mathbf{Z}_{2}\mathcal{H}_{2}+\mathcal{H}_{1}^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}\mathcal{W}_{2}^{T}+\mathcal{H}_{1}^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}_{2}\mathcal{H}_{2}\bigg{\\}}.$
(28)

Therefore, for any vector $\theta_{1}\in\mathbb{R}^{p}$ and
$\theta_{2}\in\mathbb{R}^{q}$, we have

$\theta_{1}^{T}(\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy})\theta_{2}=\theta_{1}^{T}\mathcal{W}_{1}^{T}\Big{(}\frac{\mathbf{Z}^{T}\mathbf{Z}}{n}-I_{r}\Big{)}\mathcal{W}_{2}\theta_{2}+\frac{1}{n}\theta_{1}^{T}\Big{(}\mathcal{W}_{1}\mathbf{Z}^{T}\mathbf{Z}_{2}\mathcal{H}_{2}+\mathcal{H}_{1}^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}\mathcal{W}_{2}^{T}+\mathcal{H}_{1}^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}_{2}\mathcal{H}_{2}\Big{)}\theta_{2}.$
(29)

By Bai-Yin law on eigenvalues of Wishart matrices (Bai and Yin, 1993), there
exists abolute constant $C>0$ so that for any $t>1$,

$P\bigg{(}\norm{\frac{\mathbf{Z}^{T}\mathbf{Z}}{n}-I_{r}}_{op}<t\sqrt{r/n}\bigg{)}\geq
1-2\exp(-Ct^{2}r),$

which, combined with (27), implies

$\inf_{\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}\Big{|}\theta_{1}^{T}\mathcal{W}_{1}^{T}({\mathbf{Z}^{T}\mathbf{Z}}/{n}-I_{r})\mathcal{W}_{2}\theta_{2}\Big{|}\leq
t\mathcal{B}^{2}\|\theta_{1}\|_{2}\|\theta_{2}\|_{2}\sqrt{r/n}\Big{)}\geq
1-2\exp(-Ct^{2}r).$

Now we will state a lemma which will be required to control the other terms on
the right hand side of (29).

###### Lemma 22

Suppose $\mathbf{Z}_{1}\in\mathbb{R}^{n\times p}$ and
$\mathbf{Z}_{2}\in\mathbb{R}^{n\times q}$ are independent Gaussian data
matrices. Further suppose $x\in\mathbb{R}^{p}$ and $y\in\mathbb{R}^{q}$ are
either deterministic or independent of both $\mathbf{Z}_{1}$ and
$\mathbf{Z}_{2}$. Then there exists a constant $C>0$ so that for any $t>1$,

$P\Big{(}\absolutevalue{x^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}_{2}y}>t\|x\|_{2}\|y\|_{2}\sqrt{n}\Big{)}\leq\exp(-Cn)-\exp(t^{2}/2).$

The proof of Lemma 22 follows directly setting $b=1$ in the following Lemma,
which is proved in Appendix G.4.

###### Lemma 23

Suppose $\mathbf{Z}_{1}\in\mathbb{R}^{n\times p}$ and
$\mathbf{Z}_{2}\in\mathbb{R}^{n\times q}$ are independent standard Gaussian
data matrices, and $D\in\mathbb{R}^{n\times k_{1}}$ and
$B\in\mathbb{R}^{n\times k_{2}}$ are deterministic matrices with rank $a$ and
$b$, respectively. Let $a\leq b\leq n$. Then there exists an absolute constant
$C>0$ so that for any $t\geq 0$, the following holds with probability at least
$1-\exp(-Cn)-\exp(-t^{2}/2)$:

$\|D^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}_{2}B\|_{op}\leq
C\|D\|_{op}\|B\|_{op}\sqrt{n}\max\\{\sqrt{b},t\\}.$

Lemma 22, in conjunction with (27), implies that there exists an absolute
constant $C>0$ so that

$\frac{1}{n}\Big{|}\theta_{1}^{T}\Big{(}\mathcal{W}_{1}\mathbf{Z}^{T}\mathbf{Z}_{2}\mathcal{H}_{2}+\mathcal{H}_{1}^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}\mathcal{W}_{2}^{T}+\mathcal{H}_{1}^{T}\mathbf{Z}_{1}^{T}\mathbf{Z}_{2}\mathcal{H}_{2}\Big{)}\theta_{2}\Big{|}\leq
t\mathcal{B}^{2}\|\theta_{1}\|_{2}\|\theta_{2}\|_{2}n^{-1/2}$

with $\mathbb{P}$ probability at least $1-\exp(-Cn)-\exp(t^{2}/2)$ for all
$\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$. Therefore, there
exists $C>0$ so that

$\inf_{\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}|\theta_{1}^{T}(\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy})\theta_{2}|\leq
t\sqrt{r}\mathcal{B}^{2}\|\theta_{1}\|_{2}\|\theta_{2}\|_{2}n^{-1/2}\Big{)}\geq
1-\exp(-Cn)-\exp(-Ct^{2}).$ (30)

. Note that

$T_{1}(i,k)\leq\underbrace{\Big{|}\Big{(}(\widehat{\Omega}_{n})_{k*}-({\Sigma}_{y}^{-1})_{k*}\Big{)}^{T}(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})\widehat{u}_{n,i}\Big{|}}_{T_{11}(i,k)}+\underbrace{\Big{|}\Big{(}(\widehat{\Omega}_{n})_{k*}-({\Sigma}_{y}^{-1})_{k*}\Big{)}^{T}{\Sigma}_{yx}\widehat{u}_{n,i}\Big{|}}_{T_{12}(i,k)}.$

Now suppose $\theta_{1}=(\widehat{\Omega}_{n})_{k*}-({\Sigma}_{y}^{-1})_{k*}$
and $\theta_{2}=\widehat{u}_{n,i}$. By our assumption, $\|\theta_{1}\|_{2}\leq
C_{\text{pre}}\sqrt{s({\Sigma}_{y}^{-1})(\log q)/n}$ with $\mathbb{P}$
probability $1-o(1)$ uniformly across
$\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$. We also showed that
$\|\widehat{u}_{n,i}\|_{2}\leq 2\sqrt{\mathcal{B}}$. It is not had to see that

$\displaystyle\sup_{i\in[q],k\in[r]}T_{12}(i,k)\leq
2\mathcal{B}^{3/2}C_{\text{pre}}\sqrt{s({\Sigma}_{y}^{-1})(\log q)/n}$ (31)

with $\mathbb{P}$ probability $1-o(1)$ uniformly across
$\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})$. For $T_{11}$,
observe that (30) applies because
$\theta_{i}=(\widehat{\Omega}_{n})_{k*}-({\Sigma}_{y}^{-1})_{k*}$ and
$\theta_{2}=\widehat{u}_{n,i}$ are independent of $\widehat{\Sigma}_{n,xy}$.
Thus we can write that for any $t>1$, there exists $C_{\mathcal{B}}>1$ such
that

$\sup_{\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}|T_{11}(i,k)|>tC_{\mathcal{B}}C_{\text{pre}}\sqrt{rs({\Sigma}_{y}^{-1})\log
q}/n\Big{)}\leq\exp(-Cn)+\exp(-Ct^{2}).$

Applying union bound, we obtain that

$\displaystyle\sup_{\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}\max_{i\in[q]}\max_{k\in[r]}|T_{11}(i,k)|>tC_{\mathcal{B}}C_{\text{pre}}\sqrt{rs({\Sigma}_{y}^{-1})\log
q}/n\Big{)}$ $\displaystyle\leq\exp(-Cn+\log(qr))+\exp(-Ct^{2}+\log(qr)).$

Since $r<q$ and $\log q=o(n)$, setting $t=2\sqrt{\log q}/C$, we obtain that

$\sup_{\mathbb{P}\in\mathcal{P}_{G}(r,s_{x},s_{y},\mathcal{B})}\mathbb{P}\Big{(}\max_{i\in[q]}\max_{k\in[r]}|T_{11}(i,k)|>C_{\mathcal{B}}C_{\text{pre}}\sqrt{rs({\Sigma}_{y}^{-1})}\log
q/n\Big{)}=o(1).$

Using (24) and (31), one can show that
$\epsilon_{n}=C_{\text{pre}}\sqrt{s({\Sigma}_{y}^{-1})(\log(p+q))/n}\max\\{\sqrt{r(\log
q)/n},1\\}$ in this case.

####  B.0.3 Case (C)

Note that when $\widehat{\Omega}_{n}={\Sigma}_{y}^{-1}$, $T_{1}(i,k)=0$.
Therefore, (24) implies
$\epsilon_{n}=\sqrt{{s({\Sigma}_{y}^{-1})\log(p+q)}/{n}}$ in this case.

##  C Proof of Theorem 6

Since the proof for $U$ and $V$ follows in a similar way, we will only
consider the support recovery of $U$. The proof for both cases follows a
common structure. Therefore, we will elaborate the common structure first.
Since the model $\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ is fairly large, we
will work with a smaller submodel. Specifically, we will consider a subclass
of the single spike models, i.e. $r=1$. Because we are concerned with only the
support recovery of the left singular vectors, we fix $\beta_{0}$ in
$\mathbb{R}^{q}$ so that $\|\beta_{0}\|_{2}=1$. We also fix $\rho\in(0,1)$ and
consider the subset
$\mathcal{E}\subset\\{\alpha\in\mathbb{R}^{p}:\|\alpha\|_{2}=1\\}$. Both
$\rho$ and $\mathcal{E}$ will be chosen later. We restrict our attention to
the submodel $\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})$ given by

$\bigg{\\{}\mathbb{P}\in\mathcal{P}(1,s_{x},s_{y},\mathcal{B})\ :\
\mathbb{P}\equiv N_{p+q}(0,\Sigma)\text{ where }\Sigma\text{ is of the form
}\eqref{def: sigma in Theorem 3}\text{ with
}\alpha\in\mathcal{E},\beta=\beta_{0}\bigg{\\}},$

where (32) is as follows:

$\Sigma=\begin{bmatrix}I_{p}&\rho\alpha\beta^{T}\\\
\rho\beta\alpha^{T}&I_{q}\end{bmatrix}.$ (32)

That $\Sigma$ is positive definite for $\rho\in(0,1)$ can be shown either
using elementary linear algebra or the the hidden variable representation
(26). During the proof of part (B), we will choose $\mathcal{E}$ so that
$\texttt{Sig}_{x}^{2}\leq(B^{2}-1)(\log(p-s_{x}))/8n$, which will ensure that
$\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})\subset\mathcal{P}_{\texttt{Sig}}(r,s_{x},s_{y},\mathcal{B})$
as well.

Note that for $\mathbb{P}\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})$, $U$
corresponds to $\alpha$, and hence $D(U)=D(\alpha)$. Therefore for the proof
of both parts, it suffices to show that for any decoder $\widehat{D}_{\alpha}$
of $D(\alpha)$,

$\inf_{\widehat{D}_{\alpha}}\sup_{\mathbb{P}\in\mathcal{M}(s_{x},s_{y},\mathcal{E})}\mathbb{P}\Big{(}\widehat{D}_{\alpha}\neq
D(\alpha)\Big{)}>1/2.$ (33)

In both of the proofs, our $\mathcal{E}$ will be a finite set. Our goal is to
choose $\mathcal{E}$ so that $\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})$ is
structurally rich enough to garuantee (33), yet lends itself to easy
computations. The guidance for choosing $\mathcal{E}$ comes from our main
technical tool for this proof, which is Fano’s inequality. We use the verson
of Fano’s inequality in Yatracos (1988) (Fano’s Lemma). Applied to our
problem, this inequality yields

$\inf_{\widehat{D}_{\alpha}}\sup_{\mathbb{P}\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})}\mathbb{P}\Big{(}\widehat{D}_{\alpha}\neq
D(\alpha)\Big{)}\geq
1-\dfrac{\frac{\sum_{\mathbb{P}_{1},\mathbb{P}_{2}\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})}KL(\mathbb{P}_{1}^{n}|\mathbb{P}_{2}^{n})}{|\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})|^{2}}+\log
2}{\log(|\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})|-1)},$ (34)

where $\mathbb{P}_{n}$ denotes the product measure corresponding to $n$ i.i.d.
observations from $\mathbb{P}$. We also have the following result for product
measures,
$KL(\mathbb{P}_{1}^{n}|\mathbb{P}_{2}^{n})=nKL(\mathbb{P}_{1}|\mathbb{P}_{2})$.
Moreover, when
$\mathbb{P}_{1},\mathbb{P}_{2}\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})$
with left singular vectors $\alpha_{1}$ and $\alpha_{2}$, respectively,

$KL(\mathbb{P}_{1}|\mathbb{P}_{2})=\log\frac{\text{det}(\Sigma_{2})}{\text{det}(\Sigma_{1})}-(p+q)+Tr(\Sigma_{2}^{-1}\Sigma_{1}),$

where $\text{det}(\Sigma_{1})=\text{det}(\Sigma_{2})=1-\rho^{2}$ by Lemma 32,
and

$-(p+q)+Tr(\Sigma_{2}^{-1}\Sigma_{1})=\frac{2\rho^{2}}{1-\rho^{2}}\bigg{(}1-(\alpha_{1}^{T}\alpha_{2})\|\beta_{0}\|_{2}^{2}\bigg{)}$

by Lemma 33. Noting $\alpha_{1}$, $\alpha_{2}$, and $\beta_{0}$ are unit
vectors, we derive
$KL(\mathbb{P}_{1}|\mathbb{P}_{2})=\rho^{2}(\|\alpha_{1}-\alpha_{2}\|^{2}_{2})/(1-\rho^{2})$.
Therefore, in our case, (34) reduces to

$\inf_{\widehat{D}_{\alpha}}\sup_{\mathbb{P}\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})}\mathbb{P}\Big{(}\widehat{D}_{\alpha}\neq
D(\alpha)\Big{)}\geq
1-\dfrac{n\rho^{2}\sup_{\alpha_{1},\alpha_{2}\in\mathcal{E}}\|\alpha_{1}-\alpha_{2}\|^{2}/(1-\rho^{2})+\log
2}{\log(|\mathcal{E}|-1)}.$ (35)

Thus, to ensure the right hand side of (35) is non-negligible, the key is to
choose $\mathcal{E}$ so that the $\alpha$’s in $\mathcal{E}$ are close in
$l_{2}$ norm, but $|\mathcal{E}|$ is sufficiently large. Note that the above
ensures that distinguishing the $\alpha$’s in $\mathcal{E}$ is difficult.

###  C.1 Proof of part (A)

Note that our main job is to choose $\mathcal{E}$ and $\rho$ suitably. Let us
denote

$\alpha_{0}=(\underbrace{1/\sqrt{s}_{x},\ldots,1/\sqrt{s}_{x}}_{s_{x}\text{
many }},\underbrace{0,\ldots,0}_{p-s_{x}\text{ many }}).$

We generate a class of $\alpha$’s by replacing one of the $1/\sqrt{s}_{x}$’s
in $\alpha_{0}$ by $0$, and one of the zero’s in $\alpha_{0}$ by
$1/\sqrt{s}_{x}$. A typical $\alpha$ obtained this way looks like

$\alpha=\Big{(}\underbrace{1/\sqrt{s}_{x},\ldots,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{0}},\ldots
1/\sqrt{s}_{x}}_{s_{x}\text{ many
}},\underbrace{0,\ldots,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\mathbf{1/\sqrt{s_{x}}}},\ldots,0}_{p-s_{x}\text{
many }}\Big{)}.$

Let $\mathcal{E}$ be the class, which consists of $\alpha_{0}$, and all such
resulting $\alpha$’s. Note that $|\mathcal{E}|=s_{x}(p-s_{x})$, and
$\alpha_{1},\alpha_{2}\in\mathcal{E}$ satisfy

$\|\alpha_{1}-\alpha_{2}\|^{2}_{2}\leq\|\alpha_{1}-\alpha_{0}\|^{2}_{2}+\|\alpha_{2}-\alpha_{0}\|^{2}_{2}\leq
4s_{x}^{-1}.$

Because $p>s_{x}>1$, we have $\log(s_{x}(p-s_{x})-1)\geq\log(p-s_{x})$.
Therefore, (35) leads to

$\inf_{\widehat{D}_{\alpha}}\sup_{\mathbb{P}\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})}\mathbb{P}\Big{(}\widehat{D}_{\alpha}\neq
D(\alpha)\Big{)}\geq 1-\frac{4\rho^{2}ns_{x}^{-1}/(1-\rho^{2})+\log
2}{\log(p-s_{x})},$

which is bounded below by $1/2$ whenever

$s_{x}>\frac{8\rho^{2}n}{(1-\rho^{2})\\{\log(p-s_{x})-\log 4\\}},\ \text{
which follows if }\ s_{x}>\frac{16\rho^{2}n}{(1-\rho^{2})\log(p-s_{x})}$

because $4=\sqrt{16}<\sqrt{p-s_{x}}$. To get the best bound on $s_{x}$, we
choose the value of $\rho$ which minimizes $\rho^{2}/(1-\rho^{2})$ for
$\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$, that is
$\rho=1/\mathcal{B}$. Plugging in $\rho=1/\mathcal{B}$, the proof follows.

###  C.2 Proof of part (B)

Suppose each $\alpha\in\mathcal{E}$ is of the following form

$\alpha=\Big{(}\underbrace{b,\ldots,b}_{s_{x}-1\text{ many
}},\underbrace{0,\ldots,0,{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}z},0,\ldots,0}_{p-s_{x}+1\text{
many }}\Big{)}.$

We fix $z\in(0,1)$, and hence $b=\sqrt{(1-z^{2})/(s_{x}-1)}$ is also fixed. We
will choose the value of $\rho$ and $z$ later so that
$\mathcal{P}_{\texttt{Sig}}(r,s_{x},s_{y},\mathcal{B})\supset\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})$.
Since $z$ is fixed, such an $\alpha$ can be chosen in $p-s_{x}+1$ ways.
Therefore $|\mathcal{E}|=p-s_{x}+1$. Also note that for
$\alpha,\alpha^{\prime}\in\mathcal{E}$,
$\|\alpha-\alpha^{\prime}\|_{2}^{2}\leq 2z^{2}$. Therefore (35) implies

$\inf_{\widehat{D}_{\alpha}}\sup_{P\in\mathcal{M}(s_{x},s_{y},\rho,\mathcal{E})}P\Big{(}\widehat{D}_{\alpha}\neq
D(\alpha)\Big{)}\geq 1-\dfrac{2n\rho^{2}z^{2}/(1-\rho^{2})+\log
2}{\log(p-s_{x})},$ (36)

which is greater than $1/2$ whenever

$z^{2}<\frac{1-\rho^{2}}{4n\rho^{2}}\log(\frac{p-s_{x}}{4}),\quad\text{which
holds if }\quad z^{2}=\frac{1-\rho^{2}}{8n\rho^{2}}\log(p-s_{x})$

because $16<p-s_{x}$. To get the best bound on $z$, we choose the value of
$\rho$ for $\mathbb{P}\in\mathcal{P}(r,s_{x},s_{y},\mathcal{B})$ which
maximizes $(1-\rho^{2})/\rho^{2}$, that is $\rho=1/\mathcal{B}$. Thus (33) is
satisfied when $\rho=1/\mathcal{B}$, and $\mathcal{E}$ corresponds to
$z^{2}=(\mathcal{B}^{2}-1)\log(p-s_{x})/(8n)$. Since the minimal signal
strength $\texttt{Sig}_{x}$ for any
$\mathbb{P}\in\mathcal{M}(s_{x},s_{y},\mathcal{B}^{-1},\mathcal{E})$ equals
$\min(z,b)\leq z$, we have
$\mathcal{P}_{\texttt{Sig}}(r,s_{x},s_{y},\mathcal{B})\supset\mathcal{M}(s_{x},s_{y},\mathcal{B}^{-1},\mathcal{E})$,
which completes the proof.

##  D Proof of Theorem 10

We first introduce some notations and terminologies that are required for the
proof. For $w\in\mathbb{Z}^{m}$, and $x\in\mathbb{R}^{m}$, we denote
$w!=\prod_{i=1}^{m}w_{i}!$ and $x^{w}=\prod_{i=1}^{m}x_{i}^{w_{i}}$. In low-
degree polynomial literature, when $w\in\mathbb{Z}^{m}$, the notation $|w|$ is
commonly used to denote the sum $\sum_{i=1}^{m}w_{i}$ for sake of simplicity.
We also follow the above convention. Here the notation $|\cdot|$ should not be
confused with the absolute value of real numbers. Also, for any function
$f:\mathbb{R}^{m}\mapsto\mathbb{R}$, $w\in\mathbb{Z}^{m}$, and
$t=(t_{1},\ldots,t_{m})$, we denote

$\partial_{t}^{w}f(t)=\frac{\partial^{|w|}}{\partial
t_{1}^{w_{1}}\ldots\partial t_{r}^{w_{r}}}f(t).$

We will also use the shorthand notation $\\\ E_{\pi}$ to denote
$\mathbb{E}_{\alpha\sim\pi_{x},\beta\sim\pi_{y}}$ sometimes.

Our analysis relies on the Hermite polynomial, which we will discuss here very
briefly. For a detailed account on the Hermite polynomials, see Chapter V of
Szegö (1939). The univariate Hermite polynomials of degree $k$ will be denoted
by $h_{k}$. For $k\geq 0$, the univariate Hermite polynomials
$h_{k}:\mathbb{R}\mapsto\mathbb{R}$ are defined recursively as follows:

$\displaystyle h_{0}(x)=1,\quad h_{1}(x)=xh_{0}(x),\quad\ldots,\quad
h_{k+1}(x)=xh_{k}(x)-h_{k}^{\prime}(x).$

The normalized univariate Hermite polynomials are given by
$\widehat{h}_{k}(x)=h_{k}(x)/\sqrt{k!}$. The univariate Hermite polynomials
form an orthogonal basis of $L_{2}(N(0,1))$. For $w\in\mathbb{Z}^{m}$, the
$m$-variate Hermite polynomials are given by
$H_{w}(y)=\prod_{i=1}^{m}h_{w_{i}}(y_{i})$, where $y\in\mathbb{R}^{m}$. The
normalized version $\widehat{H}_{w}$ of $H_{w}$ equals $H_{w}/\sqrt{w!}$. The
polynomials $\widehat{H}_{w}$’s form an orthogonal basis of
$L_{2}(N_{m}(0,I_{m}))$. We denote by $\Pi_{n}^{\leq D_{n}}$ the linear span
of all $n(p+q)$-variate Hermite polynomials of degree at most $D_{n}$. Since
$\mathbb{L}_{n}^{\leq D_{n}}$ is the projection of $\mathbb{L}_{n}$ on
$\Pi^{\leq D_{n}}$, it then follows that

$\|\mathbb{L}_{n}^{\leq
D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}=\sum_{\begin{subarray}{c}w\in\mathbb{Z}^{n(p+q)}\\\
|w|\leq
D_{n}\end{subarray}}\langle\mathbb{L}_{n},\widehat{H}_{w}\rangle_{L^{2}(\mathbb{Q}_{n})}^{2}.$
(37)

From now on, the degree-index vector $w$ of $\widehat{H}_{w}$ or $H_{w}$ will
be assumed to lie in $\mathbb{Z}^{n(p+q)}$. We will partition $w$ into $n$
components, which gives $w=(w_{1},\ldots,w_{n})$, where
$w_{i}\in\mathbb{Z}^{p+q}$ for each $i\in[n]$. Clearly, $i$ here corresponds
to the $i$-th observation. We also separate each $w_{i}$ into two parts
$w_{i}^{x}\in\mathbb{Z}^{p}$ and $w_{i}^{y}\in\mathbb{Z}^{q}$ so that
$w_{i}=(w_{i}^{x},w_{i}^{y})$. We will also denote
$w^{x}=(w^{x}_{1},\ldots,w^{x}_{n})$, and
$w_{y}=(w^{y}_{1},\ldots,w^{y}_{n})$. Note that $w^{x}\in\mathbb{Z}^{np}$ and
$w^{y}\in\mathbb{Z}^{nq}$, but $w\neq(w^{x},w^{y})$ in general, although
$|w|=|w^{x}|+|w^{y}|$.

Now we state the main lemmas which yields the value of $\|\mathbb{L}_{n}^{\leq
D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}$. The first lemma, proved in Appendix
G.3, gives the form of the inner products
$\langle\mathbb{L}_{n},\widehat{H}_{w}\rangle_{L^{2}(\mathbb{Q}_{n})}$.

###### Lemma 24

Suppose $w$ is as defined above and $\mathbb{L}_{n}$ is as in (15). Then it
holds that

$\langle\mathbb{L}_{n},\widehat{H}_{w}\rangle_{L^{2}(\mathbb{Q}_{n})}^{2}=\begin{cases}\frac{\mathcal{B}^{-|w|}}{w!}\bigg{\\{}\mathbb{E}_{\pi}\bigg{[}1\\{\|\alpha\|_{2}\|\beta\|_{2}<\mathcal{B}\\}\alpha^{\sum_{i=1}^{n}w^{x}_{i}}\beta^{\sum_{i=1}^{n}w^{y}_{i}}\bigg{]}\bigg{\\}}^{2}\bigg{(}\prod_{i=1}^{n}{|w^{x}_{i}|!}\bigg{)}^{2}&\\\
\quad\text{ if }\quad|w^{x}_{i}|=|w^{y}_{i}|\text{ for all }i\in[n],&\\\
0&o.w.\end{cases}$

Here the priors $\pi_{x}$ and $\pi_{y}$ are the Rademacher priors defined in
(12).

Our next lemma uses Lemma 24 to give the form of $\|\mathbb{L}_{n}^{\leq
D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}$. This lemma uses replicas of $\alpha$
and $\beta$. Suppose $\alpha_{1},\alpha_{2}\sim\pi_{x}$ and
$\beta_{1},\beta_{2}\sim\pi_{y}$ are all independent Rademacher priors, where
$\pi_{x}$ and $\pi_{y}$ are defined as in (12). We overload notation, and use
$\mathbb{E}_{\pi}$ to denote the expectation under $\alpha_{1}$, $\alpha_{2}$,
$\beta_{1}$, and $\beta_{2}$.

###### Lemma 25

Suppose $W$ is the indicator function of the event
$\\{\|\alpha_{1}\|_{2}\|\beta_{1}\|_{2}<\mathcal{B},\
\|\alpha_{2}\|_{2}\|\beta_{2}\|_{2}<\mathcal{B}\\}$. Then For any
$D_{n}\in\mathbb{N}$,

$\|\mathbb{L}_{n}^{\leq
D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}=\mathbb{E}_{\pi}\bigg{[}W\sum_{d=0}^{\left\lfloor
D_{n}/2\right\rfloor}{d+n-1\choose
d}\bigg{\\{}\mathcal{B}^{-2}(\alpha_{1}^{T}\alpha_{2})(\beta_{1}^{T}\beta_{2})\bigg{\\}}^{d}\bigg{]}.$

The proof of Lemma 25 is also deferred to Appendix G.3. We remark in passing
that the negative binomial series expansion yields

$(1-x)^{-n}=\sum_{d=0}^{\infty}{n+d-1\choose d}x^{d},\quad\text{for }|x|<1,$
(38)

whose $D_{n}$-th order truncation equals

$\bigg{(}(1-x)^{-n}\bigg{)}^{\leq D_{n}}=\sum_{d=0}^{D_{n}}{n+d-1\choose
d}x^{d}.$

Note that $W$ is non-zero if and only if
$\|\alpha_{1}\|_{2}\|\beta_{1}\|_{2}<\mathcal{B}$ and
$\|\alpha_{2}\|_{2}\|\beta_{2}\|_{2}<\mathcal{B}$, which, by Cauchy Schwarz
inequality, implies

$|(\alpha_{1}^{T}\alpha_{2})(\beta_{1}^{T}\beta_{2})|<\mathcal{B}^{2}.$

Thus $|\mathcal{B}^{-2}(\alpha_{1}^{T}\alpha_{2})(\beta_{1}^{T}\beta_{2})|<1$
when $W=1$. Hence Lemma 25 can also be written as

$\|\mathbb{L}_{n}^{\leq
D_{n}}\|^{2}_{L_{2}(\mathbb{Q}_{n})}=\mathbb{E}_{\pi}\left[W\bigg{\\{}\bigg{(}1-\mathcal{B}^{-2}(\alpha_{1}^{T}\alpha_{2})(\beta_{1}^{T}\beta_{2})\bigg{)}^{-n}\bigg{\\}}^{\leq\left\lfloor
D_{n}/2\right\rfloor}\right].$

Now we are ready to prove Theorem 10.

Proof [Proof of Theorem 10] Our first task is to get rid of $W$ from the
expression of $\|\mathbb{L}_{n}^{\leq D_{n}}\|_{L_{2}(\mathbb{Q}_{n})}$ in
Lemma 25. However, we can not directly bound $W$ by one since the term
$(\alpha_{1}^{T}\alpha_{2})^{d}(\beta_{1}^{T}\beta_{2})^{d}W$ may be negative
for odd $d\in\mathbb{N}$. We claim that
$\mathbb{E}[(\alpha_{1}^{T}\alpha_{2})^{d}(\beta_{1}^{T}\beta_{2})^{d}W]=0$ if
$d\in\mathbb{N}$ is odd. To see this, first we write

$\displaystyle\mathbb{E}\Big{[}(\alpha_{1}^{T}\alpha_{2})^{d}(\beta_{1}^{T}\beta_{2})^{d}W\Big{]}=\mathbb{E}\Big{[}\mathbb{E}\Big{[}(\alpha_{1}^{T}\alpha_{2})^{d}W\Big{|}\beta_{1},\beta_{2}\Big{]}(\beta_{1}^{T}\beta_{2})^{d}\Big{]}.$
(39)

Note that