plain

# On Statistical Inference with High Dimensional Sparse CCA

N. Laha  N. Huey  B. Coull  R. Mukherjee<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Harvard
University, 677 Huntington Ave, Boston, MA 02115, , , ,

###### Abstract

We consider asymptotically exact inference on the leading canonical
correlation directions and strengths between two high dimensional vectors
under sparsity restrictions. In this regard, our main contribution is the
development of a loss function, based on which, one can operationalize a one-
step bias-correction on reasonable initial estimators. Our analytic results in
this regard are adaptive over suitable structural restrictions of the high
dimensional nuisance parameters, which, in this set-up, correspond to the
covariance matrices of the variables of interest. We further supplement the
theoretical guarantees behind our procedures with extensive numerical studies.

###### keywords:

Sparse Canonical Correlation Analysis; Asymptotically Valid Confidence
Intervals; One-Step Bias Correction; High Dimensional Nuisance Parameters.

##  1 Introduction

Statistical analyses of biomedical applications require methods which can
handle complex data structures. In particular, to understand the relationship
between potentially high dimensional variables, formal and systematic
Exploratory Data Analysis (EDA) is often an important first step. Key examples
in this regard include, but are not limited to, eQTL mapping studies (Witten
et al., 2009; Chen et al., 2012), epigenetic studies (Holm et al., 2010; Sofer
et al., 2012; Hu et al., 2017, 2016), and in general studies involving
integration of multiple biological data such as genetic markers, gene
expressions, and disease phenotypes (Kang et al., 2013; Lin et al., 2013). Of
critical relevance in each of these examples is that of understanding
relationships between possibly high dimensional variables of interest.  In
this regard, linear relationships are the simplest, most intuitive, and lend
themselves to easy interpretations. Subsequently, a large volume of
statistical literature has been devoted to exploring linear relationships
through variants of the classical statistical toolbox of Canonical Correlation
Analysis (CCA) (Hotelling, 1992). Our focus in this paper pertains to some
fundamental inferential questions in the context of high dimensional CCA.

To formally set up the inferential questions in the CCA framework, we consider
i.i.d. data $(X_{i},Y_{i})_{i=1}^{n}\sim\mathbb{P}$ on two random vectors
$X\in\mathbb{R}^{p}$ and $Y\in\mathbb{R}^{q}$ with joint covariance matrix

$\Sigma=\begin{bmatrix}{\Sigma}_{x}&{\Sigma}_{xy}\\\
{\Sigma}_{yx}&{\Sigma}_{y}\\\ \end{bmatrix}.$

The first canonical correlation $\rho_{0}$ is defined as the maximum possible
correlation between two linear combinations of $X$ and $Y$. More specifically,
consider the following optimization problem:

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

The maximum value attained in (1) is ${\rho}_{0}$, and the solutions to (1)
are commonly referred as the first canonical directions, which we will denote
by $\alpha_{0}$ and $\beta_{0}$, respectively.  This paper considers inference
on $\alpha_{0}$, $\beta_{0}$, and associated quantities of interest. In most
scientific problems, the first canonical correlation coefficient is of prime
interest as it summarizes the “maximum linear association” between $X$ and $Y$
and thereby motivating our choice of inferential target.

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, 1962) and references therein for further details. These
classical results have been thereafter heavily used to provide statistical
inference (i.e. asymptotically valid hypotheses tests, confidence intervals
and P-values) across a vast canvas of disciplines such as psychology,
agriculture, oceanography and others. However, modern surge in interests for
CCA, often being motivated by data from high throughput biological
experiments, requires re-thinking several aspects of the traditional theory
and methods. In particular, in most modern data examples, the number of
samples is typically comparable to or much smaller than the number of
variables in the study – rendering the classical CCA inconsistent and
inadequate without further structural assumptions Cai et al. (2018); Ma et al.
(2020); Bao et al. (2019). A natural structural constraint that has gained
popularity in this regard, is that of sparsity i.e. the phenomenon of an
(unknown) few collection of variables being related to each other rather than
contributions to the associations from the whole collection of high
dimensional components. The framework of Sparse Canonical Correlation Analysis
(SCCA) (Witten et al., 2009) has thereafter been developed to target such low
dimensional structures, and to subsequently provide consistent estimation in
the context of high dimensional CCA. Although such structured CCA problems
have witnessed a renewed enthusiasm from both theoretical and applied
communities, most papers have heavily focused on key aspects of estimation (in
suitable norms) and relevant scalable algorithms – se e.g. Chen et al. (2013);
Gao et al. (2015, 2017); Ma et al. (2020); Mai & Zhang (2019). However,
asymptotically valid inference is yet to be explored systematically in the
context of SCCA. In particular, none of the existing estimation methods for
SCCA lend themselves to uncertainty quantification, i.e. inference on
$\alpha_{i}$ $(i=1,\ldots,p)$, $\beta_{j}$ $(j=1,\ldots,q)$, or $\rho$. This
is indeed not surprising, since being based on penalized procedures, existing
estimators are asymptotically biased, super-efficient for estimating $0$
coordinates, and not tractable in terms of estimating underlying asymptotic
distribution Leeb & Pötscher (2005, 2006, 2008); Pötscher & Leeb (2009).
Therefore, construction of asymptotically valid confidence intervals for
$\alpha_{i}$, $\beta_{j}$’s or $\rho$ is not straightforward. In absence of
such intervals, bootstrap or permutation tests are typically used in practice
(Witten et al., 2009). However, these methods are often empirically justified
and even then might suffer from subtle pathological issues that underlie
standard re-sampling techniques in penalized estimation framework Chatterjee &
Lahiri (2010, 2011, 2013). This paper is motivated by taking a step in
resolving these fundamental issues with inference on SCCA.

###  1.1 Main contribution

The main results of this paper is the construction of asymptotically valid
confidence intervals for $\sqrt{\rho_{0}}\alpha_{0}$ and
$\sqrt{\rho_{0}}\beta_{0}$. Our method is based on a one-step bias-correction
performed on preliminary estimators of the canonical directions. The resulting
bias-corrected estimators have an asymptotic linear influence function type
expansion (see e.g. Tsiatis (2007) for asymptotic influence function
expansions) with $\sqrt{n}$-scaling (see Theorem 4.1 and Proposition 4.3)
under suitable sparsity conditions on the truth. This representation is
subsequently exploited to build confidence intervals for a variety of relevant
lower dimensional functions of the top canonical directions; see Corollary 4.2
and Corollary 4.4 and the discussions that follow. Finally, we will show that
the entire de-biased vector is asymptotically equivalent to a high dimensional
Gaussian vector in a suitably uniform sense; see Proposition 4.3, which
enables the control of familywise error rate.

The bias correction procedure crucially relies on a novel representation of
$\sqrt{\rho_{0}}\alpha_{0}$ and $\sqrt{\rho_{0}}\beta_{0}$ as the unique
maximizers (up to a sign flip) of a smooth objective (see Lemma 3.2), which
may be of independent interest. The uniqueness criteria is indispensable here
since otherwise a crucial local convexity convexity property (see Lemma 3.4),
which we fundamentally exploit to deal with high dimensionality of the
problem, is not guaranteed. We also discuss why the commonly used
representations of the top canonical correlations is difficult to work with
owing to either the lack of such local convexity properties, or the
flexibility of its form to offer a non-cumbersome derivation of the one-step
bias correction. We elaborate on these subtleties in Section 3.2 for details.

Further, we pay special attention to adapt to underlying sparsity structures
of the marginal precision matrices ($\Sigma_{x}^{-1},\Sigma_{y}^{-1}$) of the
high dimensional variables ($X,Y$) under study – which serve as high
dimensional nuisance parameters in the problem. Consequently, our construction
of asymptotically valid confidence intervals for top canonical correlation
strength and directions are agnostic over the structures (e.g. sparsity of the
precision matrices of $X$ and $Y$) of these complex nuisance parameters. The
de-biasing procedure can be implemented using our R package de.bias.CCA
available at _https://github.com/nilanjanalaha/de.bias.CCA_.

Finally, we supplement our methods for inference with suitable constructions
of initial estimators of canonical correlation directions as well as nuisance
parameters under suitable sparsity assumptions. The construction of these
estimators, although motivated by existing ideas, requires careful
modifications to tackle inference on the first canonical correlation strength
and directions – while treating remaining directions as nuisance parameters.

##  2 Mathematical Formalism

In this section we collect some assumptions and notation that will be used
throughout the rest of the paper.

###  2.1 Structural Assumptions

Throughout this paper, we will assume that $X$ and $Y$ are centered sub-
Gaussian random vectors 111 see Vershynin (2010) for more details. with joint
covariance matrix $\Sigma$ as described above. We will let ${\Sigma}_{xy}$ to
have a fixed rank $r\geq 1$ (implying that apart from $\rho_{0}$, there are
$r-1$ additional canonical correlations Anderson, 2003). Since the cross-
covariance matrix ${\Sigma}_{xy}$ has rank $r$, it it can be shown that (cf.
Chen et al., 2013; Gao et al., 2017)

${\Sigma}_{xy}={\Sigma}_{x}U\Lambda V^{T}{\Sigma}_{y},$ (2)

where $U=[u_{1}\ldots u_{r}]$ and $V=[v_{1}\ldots v_{r}]$ are $p\times r$ and
$q\times r$ dimensional matrices satisfying $U^{T}{\Sigma}_{x}U=I_{r}$ and
$V^{T}{\Sigma}_{y}V=I_{r}$, respectively. The $\Lambda$ in (2) is a diagonal
matrix, whose diagonal entries are the canonical correlations, i.e.

${\rho}_{0}=\Lambda_{1}\geq\Lambda_{2}\geq\ldots\geq\Lambda_{r}>0.$

In this regard, the matrices $U$ and $V$ need not be unique unless the
canonical correlations, i.e. the $\Lambda_{i}$’s, are all unique. Indeed, we
will at the least require uniqueness of $\alpha_{0}$ and $\beta_{0}$, since
otherwise they are not even identifiable. To that end, we will make the
following assumption that is common in the literature since it grants
uniqueness of $\alpha_{0}$ and $\beta_{0}$ up to a sign flip (cf. Chen et al.,
2013; Gao et al., 2017; Mai & Zhang, 2019). [Eigengap Assumption] There exists
$\epsilon_{0}\in(0,1)$ so that ${\rho}_{0}-\Lambda_{1}>\epsilon_{0}$ for all
$n$.  Note that Assumption 2.1 also implies that ${\rho}_{0}$ stays bounded
away from zero. We will further assume that ${\Sigma}_{x}$ and ${\Sigma}_{y}$
are positive definite and bounded in operator norm. [Bounded eigenvalue
Assumption] There exists $M>0$ such that the eigenvalues of ${\Sigma}_{x}$ and
${\Sigma}_{y}$ are bounded below by $M^{-1}$ and bounded above by $M$.  This
regularity assumption is also common in the literature of SCCA (Gao et al.,
2017, 2015; Mai & Zhang, 2019; Laha & Mukherjee, 2021).

###  2.2 Notation

We will denote the set of all positive integers by $\mathbb{N}$. For a matrix
$A$, we denote its $j$th column by $A_{j}$. Also, let $\Lambda_{max}(A)$ and
$\Lambda_{min}(A)$ denote the largest and smallest eigenvalue of $A$,
respectively. We denote the gradient of a function $f$ by $\dot{f}$ or
$\gradient f$, where we reserve the notation $\gradient^{2}f$ for the hessian.
The $i$th element of any vector $v$ is denoted by $v_{i}$. We use the notation
$\|\cdot\|_{p}$ to denote the usual $l_{p}$ norm of a vector for any
$p\in\mathbb{N}$. For a matrix $A\in\mathbb{R}^{p\times q}$, $\norm{A}_{F}$
and $\norm{A}_{op}$ will denote the Frobenius and the operator norm,
respectively. We denote by $|A|_{\infty}$ the elementwise supremum of $A$.
Throughout the paper, $C$ will be used to denote a positive constant whose
value may change from line to line.

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 denote the set of the indices of the non-zero rows in $U$ and $V$ by
$S_{U}$ and $S_{V}$, respectively. We let $s_{U}$ and $s_{V}$ be the
cardinalities of $S_{U}$ and $S_{V}$ and use $s=s_{U}+s_{V}$ to denote the
total sparsity. We further denote by $s_{x}$ and $s_{y}$ the number of nonzero
elements of $\alpha_{0}$ and $\beta_{0}$, respectively. The supports of
$\alpha_{0}$ and $\beta_{0}$ will be similarly be denoted by $S_{x}$ and
$S_{y}$, respectively. We will discuss the precise requirements on these
sparsities, and the necessities of such assumptions in detail in Section 4.1.

Our method requires initial estimators of $\alpha_{0}$, $\beta_{0}$, and
$\rho_{0}$. We let $\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$ be the
initial estimators of $\alpha_{0}$ and $\beta_{0}$, respectively. Also, 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. The estimate $\widehat{\rho}_{n}$ of
${\rho}_{0}$ is

$\widehat{\rho}_{n}=\frac{\widehat{\alpha}_{n}^{T}\widehat{\Sigma}_{n,xy}\widehat{\beta}_{n}}{(\widehat{\alpha}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{\alpha}_{n})^{1/2}(\widehat{\beta}_{n}^{T}\widehat{\Sigma}_{n,y}\widehat{\beta}_{n})^{1/2}}.$
(3)

The quantity $\widehat{\rho}_{n}$ may not be positive for any
$\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$. Therefore, mostly we will
use $|\widehat{\rho}_{n}|$ as an estimate of ${\rho}_{0}$. Finally, for the
sake of simplicity, we let $\lambda$ denote the term

$\lambda=\bigg{(}\dfrac{\log(p\vee q)}{n}\bigg{)}^{1/2}.$ (4)

##  3 Methodology

In this section we discuss the intuitions and details of our main proposed
methodology that we will analyze in later sections. The discussions are
divided across three main subsections. The first Subsection 3.1 presents the
driving intuitions behind obtaining general de-biased estimators of generic
parameters of interest that can be defined through generic optimization
framework. Subsequently, Subsection 3.2 translates this intuition to a working
principle in the context of SCCA. In particular, we design a suitable
optimization criterion which allows a principled application of the general
de-biasing method and additionally lends itself to rigorous theoretical
analyses. Finally, our last Subsection 3.3 elaborates on the benefit of
designing this specific optimization objective function over other possible
choices of optimization problems for defining the leading canonical
directions.

###  3.1 The Debiasing Method in General

We first discuss the simple intuition behind reducing the bias of estimators
defined through estimating equations. To that end, suppose we are interested
in estimating $\theta_{0}\in\mathbb{R}^{p}$, which minimizes the function
$f:\mathbb{R}^{p}\mapsto\mathbb{R}$. If $f$ is smooth, then $\theta_{0}$
solves the equation $\dot{f}(\theta)=0$. Suppose $\theta$ is in a small
neighborhood of $\theta_{0}$. the Taylor series expansion of $f(\theta)$
around $\theta_{0}$ yields
$\dot{f}(\theta)-\dot{f}(\theta_{0})=\gradient^{2}f(\overline{\theta})(\theta-\theta_{0})$,
where $\overline{\theta}\in\mathbb{R}^{p}$ lies on the line segment joining
$\theta_{0}$ and $\theta$. If $f$ has finitely many global minimums, then $f$
can not be flat at $\theta_{0}$. In that case, $f$ is strongly convex at some
neighborhood of $\theta_{0}$. Therefore $\gradient^{2}f(\overline{\theta})$ is
positive definite, leading to
$\theta_{0}=\theta-(\gradient^{2}f(\overline{\theta}))^{-1}\dot{f}(\theta)$.
Suppose $\widehat{\theta}_{n}$ and $\widehat{\Phi}_{n}$ are reliable
estimators of $\theta_{0}$ and $(\gradient^{2}f(\theta_{0}))^{-1}$,
respectively. Correcting the first order bias of $\widehat{\theta}_{n}$ then
yields the de-biased estimator
$\widehat{\theta}^{db}_{n}=\widehat{\theta}_{n}-\widehat{\Phi}_{n}\dot{f}(\widehat{\theta}_{n})$.
Thus, to find a bias-corrected estimator of $\theta_{0}$, it suffices to find
a smooth function which is minimized at $\theta_{0}$ and has at most finitely
many global minima. This simple intuition is the backbone of our strategy.

######  Remark 3.1 (Positive definiteness of $\gradient^{2}f(\theta_{0})$).

The positive definiteness of $\gradient^{2}f(\theta_{0})$ is important because
most existing methods for estimating the inverse of a high dimensional matrix
requires the matrix to be positive definite. These methods proceed via
estimating the columns of $\Sigma^{-1}$ separately through a quadratic
optimization step. Unless the original matrix is positive definite, these
intermediate optimization problems are unbounded. Therefore, the algorithms
are likely to diverge even with enough observations. For more details, see
Section 1 of Janková & van de Geer (2018) (see also Section 2.1 of Yuan,
2010).

###  3.2 The Debiasing Method for SCCA

To operationalize the intuition described above in Section 3.1, we begin with
a lemma which represents ${\rho}_{0}^{1/2}\alpha_{0}$ and
${\rho}_{0}^{1/2}\beta_{0}$ as the unique minimizers (upto a sign flip) of a
smooth objective function. We defer the proof of Lemma 3.2 to Supplement 14.

######  Lemma 3.2.

For any $C>0$, we have

$\pm({\rho}_{0}^{1/2}\alpha_{0},{\rho}_{0}^{1/2}\beta_{0})=\operatorname*{arg\,min}_{x\in\mathbb{R}^{p},y\in\mathbb{R}^{q}}h(x,y).$

where
$h(x,y)=(1-C/2)(x^{T}{\Sigma}_{x}x)(y^{T}{\Sigma}_{y}y)+C(x^{T}{\Sigma}_{x}x)^{2}/4+C(y^{T}{\Sigma}_{y}y)^{2}/4-2x^{T}{\Sigma}_{xy}y$.

The proof of Lemma 3.2 hinges on a seminal result on low rank matrix
approximation dating back to Eckart & Young (1936), which implies that for any
matrix $A$ with singular value decomposition
$\sum_{i=1}^{r}\Lambda_{i}\tilde{u}_{i}\tilde{v}_{i}^{T}$,

$\sum_{i=1}^{k}\Lambda_{i}\tilde{u}_{i}\tilde{v}_{i}^{T}=\operatorname*{arg\,min}_{B\in\mathcal{M}_{k}}\|A-B\|_{F}^{2}\quad(k=1,\ldots,r),$
(5)

where $\mathcal{M}_{k}$ is the set of all $p\times q$ matrices with rank $k$.
Our main inferential method for leading canonical directions builds on Lemma
3.2, and consequently, corrects for the bias of estimating
$x^{0}={\rho}_{0}^{1/2}\alpha_{0}$ and $y^{0}={\rho}_{0}^{1/2}\beta_{0}$ using
preliminary plug-in estimators from literature. It is worth noting that we
focus on the the leading canonical directions up to a multiplicative factor
since from our inferential point of view, this quantity is enough to explore
the nature of projection operators onto these directions. in particular, for
the sake of constructing tests for no-signal such as
$H_{0}:(\alpha_{0})_{i}=0$ it is equivalent to the test $H_{0}:x^{0}_{i}=0$.

######  Remark 3.3.

Suppose $h$ is as in Lemma 3.2. It can be shown that the other stationary
points of $h(x,y)$, to be denoted by $(\tilde{x}_{i},\tilde{y}_{i})$,
correspond to the canonical pairs with correlations $\Lambda_{i}$, $(i\geq
2)$. Moreover, the Hessian of $h(x,y)$ at $(\tilde{x}_{i},\tilde{y}_{i})$ has
both positive and negative eigenvalues, indicating that the function is
neither concave nor convex at these points. Therefore, all these stationary
points are saddle points. Consequently, any minimum of $h(x,y)$ is a global
minimum – irrespective of the choice of $C>0$.

Now note that

$\displaystyle{\partialderivative{h}{x}}(x,y)=$ $\displaystyle\
(2-C)(y^{T}{\Sigma}_{y}y){\Sigma}_{x}x+C(x^{T}{\Sigma}_{x}x){\Sigma}_{x}x-2{\Sigma}_{xy}y,$
$\displaystyle{\partialderivative[2]{h}{x}}(x,y)=$ $\displaystyle\
(2-C)(y^{T}{\Sigma}_{y}y){\Sigma}_{x}+C(x^{T}{\Sigma}_{x}x){\Sigma}_{x}+2C{\Sigma}_{x}xx^{T}{\Sigma}_{x},$
$\displaystyle{\partialderivative{h}{x}{y}}(x,y)=$ $\displaystyle\
2(2-C){\Sigma}_{x}xy^{T}{\Sigma}_{y}-2{\Sigma}_{xy},$ (6)

and hence by symmetry, the Hessian $H(x,y)$ of $h$ at $(x,y)$ is given by

$\begin{bmatrix}(2-C)(y^{T}{\Sigma}_{y}y){\Sigma}_{x}+C(x^{T}{\Sigma}_{x}x){\Sigma}_{x}&2(2-C){\Sigma}_{x}xy^{T}{\Sigma}_{y}-2{\Sigma}_{xy}\\\
+2C{\Sigma}_{x}xx^{T}{\Sigma}_{x}&\\\
2(2-C){\Sigma}_{y}yx^{T}{\Sigma}_{x}-2{\Sigma}_{yx}&(2-C)(x^{T}{\Sigma}_{x}x){\Sigma}_{y}+C(y^{T}{\Sigma}_{y}y){\Sigma}_{y}\\\
&+2C{\Sigma}_{y}yy^{T}{\Sigma}_{y}\end{bmatrix}.$

At this point we note the flexibility of our approach in choosing $C$ so as to
being able to work with a relatively amenable form of the Hessian and its
inverse that we need to estimate. We subsequently set $C=2$ so that the
estimation of the cross term ${\Sigma}_{x}xy^{T}{\Sigma}_{y}$ can be avoided.
In particular, when $x^{0}={{\rho}_{0}^{1/2}}\alpha_{0}$ and
$y^{0}={{\rho}_{0}^{1/2}}\beta_{0}$, then
$(x^{0})^{T}{\Sigma}_{x}x^{0}=(y^{0})^{T}{\Sigma}_{y}(y^{0})={\rho}_{0}$. We
denote the Hessian in this case as

$H^{0}=H(x,y):=2{\rho}_{0}\begin{bmatrix}{\Sigma}_{x}+2{\Sigma}_{x}\alpha_{0}\alpha_{0}^{T}{\Sigma}_{x}&-{\Sigma}_{xy}/{\rho}_{0}\\\
-{\Sigma}_{yx}/{\rho}_{0}&{\Sigma}_{y}+2{\Sigma}_{y}\beta_{0}\beta_{0}^{T}{\Sigma}_{y}\end{bmatrix}.$
(7)

A plug-in estimator $\widehat{H}_{n}(x,y)$ of $H^{0}$ is given by

$\widehat{H}_{n}(x,y)=2\begin{bmatrix}(x^{T}\widehat{\Sigma}_{n,x}x)\widehat{\Sigma}_{n,x}+2\widehat{\Sigma}_{n,x}xx^{T}\widehat{\Sigma}_{n,x}&-\widehat{\Sigma}_{n,xy}\\\
-\widehat{\Sigma}_{n,yx}&(y^{T}\widehat{\Sigma}_{n,y}y)\widehat{\Sigma}_{n,y}+2\widehat{\Sigma}_{n,y}yy^{T}\widehat{\Sigma}_{n,y}\end{bmatrix}.$

Because our $h$ is a sufficiently well-behaved function, it possesses a
positive definite Hessian at the minima $\pm(x^{0},y^{0})$, thereby
demonstrating the crucial strong convexity property mentioned in Remark 3.1.
This property of $H^{0}$ is the content of our following lemma, the proof of
which can be found in Supplement 14.

######  Lemma 3.4.

Under Assumptions 2.1 and 2.1, the matrix $H^{0}$ defined in (7) is positive
definite with minimum eigenvalue $\Lambda_{min}(H^{0})\geq
2({\rho}_{0}-\Lambda_{2})/M$ where $M$ is as in Assumption 2.1.

Lemma 3.2 and Lemma 3.4 subsequently allows us to constructed de-biased
estimators of the leading canonical directions as follows. Suppose
$\widehat{x}_{n}={|\widehat{\rho}_{n}|}^{1/2}\widehat{\alpha}_{n}$ and
$\widehat{y}_{n}={|\widehat{\rho}_{n}|}^{1/2}\widehat{\beta}_{n}$ are
estimators of $x^{0}$ and $y^{0}$, where $\widehat{\alpha}_{n}$ and
$\widehat{\beta}_{n}$ are the preliminary estimators of $\alpha_{0}$ and
$\beta_{0}$, and $\widehat{\rho}_{n}$ is as defined in (3). Our construction
of de-biased estimators in SCCA now relies on two objects: (a) estimators of
$\partial h(\widehat{x}_{n},\widehat{y}_{n})/\partial x$ and $\partial
h(\widehat{x}_{n},\widehat{y}_{n})/\partial y$, which are simply given by

$\displaystyle{\partialderivative{\widehat{h}_{n}}{x}}(\widehat{x}_{n},\widehat{y}_{n})=$
$\displaystyle\
2(\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{x}_{n})\widehat{\Sigma}_{n,x}\widehat{x}_{n}-2\widehat{\Sigma}_{n,xy}\widehat{y}_{n},$
$\displaystyle{\partialderivative{\widehat{h}_{n}}{\widehat{y}_{n}}}(\widehat{x}_{n},\widehat{y}_{n})=$
$\displaystyle\
2(\widehat{y}_{n}^{T}\widehat{\Sigma}_{n,y}\widehat{y}_{n})\widehat{\Sigma}_{n,y}\widehat{y}_{n}-2\widehat{\Sigma}_{n,yx}\widehat{x}_{n},$
(8)

and (b) an estimator $\widehat{\Phi}_{n}$ of $\Phi^{0}$ – the inverse of
$H^{0}$. Construction of such an estimator is can be involved and to tackle
this we develop a version of the Node-wise Lasso algorithm (see Supplement
10.4 for details) popularized in recent research van de Geer et al. (2014).
Following the intuitions discussed in Section 3.1, we can then complete the
construction of the de-biased estimators, whose final form writes as

$\displaystyle\begin{bmatrix}\widehat{x}^{db}_{n}\\\
\widehat{y}^{db}_{n}\end{bmatrix}=\begin{bmatrix}\widehat{x}_{n}\\\
\widehat{y}_{n}\end{bmatrix}-\widehat{\Phi}_{n}^{T}\begin{bmatrix}\partialderivative{\widehat{h}_{n}}{x}\\\
\partialderivative{\widehat{h}_{n}}{y}\end{bmatrix}.$ (9)

In Supplement 11, we will discuss how our proposed method connects to the
broader scope of de-biased inference in high dimensional problems. In regard
to the targets of our estimators, we note that if $\widehat{\alpha}_{n}$
estimates $\alpha_{0}$, then $\widehat{x}^{db}_{n}$ also estimates $x^{0}$.
However, if $\widehat{\alpha}_{n}$ approximates $-\alpha_{0}$ instead, then
$\widehat{x}^{db}_{n}$ instead approximates $-x^{0}$. The similar phenomenon
can be observed for $\widehat{\beta}_{n}$ as well. Our theoretical analyses of
these estimators will be designed accordingly.

At this time, we are also ready to construct a de-biased estimator of
${\rho}_{0}^{2}$. To that end, suppose $\widehat{x}_{n}$ and $\widehat{y}_{n}$
are such that $\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,xy}\widehat{y}_{n}\geq
0$. Note that if that is not the case, we can always switch $\widehat{x}_{n}$
to $-\widehat{x}_{n}$ so that
$\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,xy}\widehat{y}_{n}\geq 0$. Our
estimator of ${\rho}_{0}^{2}$ can then be constructed as
$\widehat{\rho}_{n}^{2,db}=\min(1,|\widehat{\rho}_{n}^{2,\text{raw}}|)$, where

$\widehat{\rho}_{n}^{2,\text{raw}}=\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,xy}\widehat{y}^{db}_{n}+(\widehat{x}^{db}_{n})^{T}\widehat{\Sigma}_{n,xy}\widehat{y}_{n}-\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,xy}\widehat{y}_{n}.$

Before moving onto the theoretical properties of our proposed methods, we make
a slight relevant digression by noting that there are many ways to formulate
the optimization program in (1) so that $\pm(\alpha_{0},\beta_{0})$ can be
characterized as the global optimizer. We therefore close this current section
with a discussion on why the particular formulation in Lemma 3.2 particularly
useful for our purpose.

###  3.3 Subtleties with Other Representations of $\alpha_{0}$ and $\beta_{0}$

Indeed, the most intuitive approach to characterize
$\pm(\alpha_{0},\beta_{0})$ is to see it as the maximizer of the constrained
maximization problem (1). This leads to the Lagrangian

$\displaystyle
L(\alpha,\beta,l_{1},l_{2})=-\alpha^{T}{\Sigma}_{xy}\beta+l_{1}(\alpha^{T}{\Sigma}_{x}\alpha-1)+l_{2}(\beta^{T}{\Sigma}_{y}\beta-1),$
(10)

where $l_{1}$ and $l_{2}$ are the Lagrange multipliers. Denoting
$\theta=(\alpha,\beta,l_{1},l_{2})$, it can be verified that since
$\theta_{0}=(\alpha_{0},\beta_{0},{\rho}_{0}/2,{\rho}_{0}/2)$ is a stationary
point of (1), $\theta_{0}$ also solves $\dot{L}(\theta)=0$. Using the first
order Taylor series expansion of $L$, one can subsequently show that any
$\theta$ in a small neighborhood of $\theta_{0}$ has the approximate expansion

$\theta-\theta_{0}\approx\ddot{L}(\theta_{0})^{-1}\dot{L}(\theta).$

If we then replace $\theta$ by an estimator of $\theta_{0}$, one can use the
above expansion to estimate the first order bias of this estimator provided
$\ddot{L}(\theta_{0})$ is suitably nice and estimable. However, by strong max-
min property (cf. Section 5.4.1 Boyd et al., 2004), $L$ satisfies

$\sup_{l_{1},l_{2}\in\mathbb{R}}\inf_{\alpha\in\mathbb{R}^{p},\beta\in\mathbb{R}^{q}}L(\alpha,\beta,l_{1},l_{2})=\inf_{\alpha\in\mathbb{R}^{p},\beta\in\mathbb{R}^{q}}\sup_{l_{1},l_{2}\in\mathbb{R}}L(\alpha,\beta,l_{1},l_{2}),$
(11)

which implies $(\alpha_{0},\beta_{0},{\rho}_{0}/2,{\rho}_{0}/2)$ is a saddle
point of $L$. Thus $\ddot{L}(\theta_{0})$ fails to be positive definite. In
fact, any constrained optimization program fails to provide a Lagrangian with
positive definite hessian, and thus violates the requirements outlined in
Section 3.1. We have already pointed out in Remark 3.1 that statistical tools
for efficient estimation of the inverse of a high dimensional matrix is scarce
unless the matrix under consideration is positive definite. Therefore, we
refrain from using the constrained optimization formulation in (1) for the de-
biasing procedure.

For any $C>0$, the function

$f:(\alpha,\beta)\mapsto-\frac{\alpha^{T}{\Sigma}_{xy}\beta}{(\alpha^{T}{\Sigma}_{x}\alpha)^{1/2}(\beta^{T}{\Sigma}_{y}\beta)^{1/2}}+C(\alpha^{T}{\Sigma}_{x}\alpha-1)^{2}+C(\beta^{T}{\Sigma}_{y}\beta-1)^{2},$

however, is a valid choice for the $f$ outlined in Subsection 3.1 since its
only global minimizers are $\pm(\alpha_{0},\beta_{0})$, which also indicates
strong convexity at $\pm(\alpha_{0},\beta_{0})$. However, the gradient and the
Hessian of this function takes a complicated form. Therefore, establishing
asymptotic results for the de-biased estimator based on this $f$ is
significantly more cumbersome than its counterpart based on the $h$ in Lemma
3.2. Hence, we refrain from using this objective function for our de-biasing
procedure as well.

##  4 Asymptotic Theory for the De-biased Estimator

In this section we establish theoretical properties of our proposed estimators
under a high dimensional sparse asymptotic framework. To set up our main
theoretical results, we first present assumptions on sparsities of the true
canonical directions and desired conditions on initial estimators of
$\alpha_{0},\beta_{0},\Phi^{0}$ in Subsection 4.1. The construction of
estimators with these desired properties are discussed in Appendices 1 and 2.
Subsequently, we present the main asymptotic results and its implications for
construction of confidence intervals of relevant quantities of interest in
Subsection 4.2.

###  4.1 Assumptions on $\widehat{\alpha}_{n}$, $\widehat{\beta}_{n}$, and
$\widehat{\Phi}_{n}$

For the de-biasing procedure to be successful, it is important that
$\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$ are both $l_{1}$ and $l_{2}$
consistent for $\alpha_{0}$ and $\beta_{0}$ with suitable rates of
convergence. In particular, we will require them to satisfy the following
condition. [Preliminary estimator condition] The preliminary estimators
$\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$ of $\alpha_{0}$ and
$\beta_{0}$ satisfy the followings for some $\kappa\in[1/2,1]$,
$s=s_{U}+s_{V}$, and $\lambda$ as defined in (4):

$\inf_{w\in\\{\pm
1\\}}\|w\widehat{\alpha}_{n}-\alpha_{0}\|_{2}+\inf_{w\in\\{\pm
1\\}}\|w\widehat{\beta}_{n}-\beta_{0}\|_{2}=O_{p}(s^{\kappa}\lambda),$

and

$\inf_{w\in\\{\pm
1\\}}\|w\widehat{\alpha}_{n}-\alpha_{0}\|_{1}+\inf_{w\in\\{\pm
1\\}}\|w\widehat{\beta}_{n}-\beta_{0}\|_{1}=O_{p}(s^{\kappa+1/2}\lambda).$

We present discussions regarding the necessity of the rates presented above as
well as the motivation behind the exponent $\kappa\in[1/2,1]$ in Supplement
12. Moreover, we also discuss the construction of estimators satisfying
Condition 4.1 in Supplement 9. Our method for developing these initial
estimators is motivated by the recent results in Gao et al. (2017), who
jointly estimate $U$ and $V$ up to an orthogonal rotation with desired $l_{2}$
guarantees. However, our situation is somewhat different since we need to
estimate $\alpha_{0}$ and $\beta_{0}$ up to a sign flip, which might not be
obtained from the joint estimation of all the directions up to orthogonal
rotation. This is an important distinction since the remaining directions act
as nuisance parameters in our set up. It turns out that the asymptotics of the
sign-flipped version requires crucial modification of the arguments of Gao et
al. (2017). The analysis of this modified procedure presented in Supplement 1
in turn allows us to extract both the desired $l_{1}$ and $l_{2}$ guarantees
in the process.

We will also require an assumption on the sparsities $s_{U}$ and $s_{V}$, the
number of nonzero rows of $U$ and $V$, respectively. We present this next
while deferring the discussions on the necessity of such assumptions to
Appendix 12.  [Sparsity Assumption] We assume $s_{U}=o(p)$, $s_{V}=o(q)$, and
$s^{2\kappa}\lambda^{2}=o(n^{-1/2})$ where $s=s_{U}+s_{V}$ and $\kappa$ is as
in Condition 4.1.

Finally, our last condition pertains to the estimator $\widehat{\Phi}_{n}$ on
$\Phi^{0}$. Most methods for estimating precision matrices can be adopted to
estimate $\Phi^{0}$ using an estimator of $H^{0}$. However, care is needed
since $\widehat{\Phi}_{n}$ needs to satisfy some rates of convergence for the
de-biased estimators in (9) to be $\sqrt{n}$-consistent. We collect this
condition below. [Inverse hessian Conditions] The estimator
$\widehat{\Phi}_{n}$ satisfies

$\max_{1\leq j\leq
p+q}\|(\widehat{\Phi}_{n})_{j}-\Phi^{0}_{j}\|_{1}=O_{p}(s^{\kappa+1/2}\lambda),$

and

$\max_{1\leq j\leq
p+q}\|(\widehat{\Phi}_{n})_{j}-\Phi^{0}_{j}\|_{2}=O_{p}(s^{\kappa}\lambda),$

where $\kappa$ is as in Condition 4.1.  We defer the discussion on the
construction of $\widehat{\Phi}_{n}$ to Appendix 10, where, in particular, we
will show that the a nodewise Lasso type estimator, which appeals to the ideas
in van de Geer et al. (2014), satisfies Condition 4.1.

###  4.2 Theoretical Analyses

In what follows, we only present the results on inference for $\alpha_{0}$.
Parallel results for $\beta_{0}$ can be obtained similarly. Before stating the
main theorem, we introduce a few additional notation. We partition the
$i^{\mathrm{th}}$ column of $\Phi^{0}$ comfortably w.r.t. the dimensions of
$X$ and $Y$ as $\Phi^{0}_{i}=(\Phi^{0}_{i,1},\Phi^{0}_{i,2})$ where
$\Phi^{0}_{i,1}\in\mathbb{R}^{p}$ and $\Phi^{0}_{i,2}\in\mathbb{R}^{q}$. We
subsequently define the random variable

$\displaystyle\mathcal{Z}(i)=$ $\displaystyle\
[{\rho}_{0}(\Phi^{0}_{i,1})^{T}+\\{(\Phi^{0}_{i,1})^{T}{\Sigma}_{x}x^{0}\\}(x^{0})^{T}]XX^{T}x^{0}+[{\rho}_{0}(\Phi^{0}_{i,2})^{T}+\\{(\Phi^{0}_{i,2})^{T}{\Sigma}_{y}y^{0}\\}(y^{0})^{T}]YY^{T}y^{0}$
$\displaystyle\
-(\Phi^{0}_{i,1})^{T}XY^{T}y^{0}-(x^{0})^{T}XY^{T}\Phi^{0}_{i,2},$ (12)

and its associated variance as

$\displaystyle\sigma_{i}^{2}=\text{var}(\mathcal{Z}(i)).$ (13)

Since $X$ and $Y$ are sub-Gaussian, it can be shown that all moments of
$\mathcal{Z}(i)$, and in particular, the $\sigma^{2}_{i}$’s are finite under
Assumption 2.1. Indeed, we show the same through the proof of Theorem 4.1 .
Finally define

$\displaystyle\mathcal{L}=\begin{bmatrix}\mathcal{L}_{(1)}\\\
\mathcal{L}_{(2)}\end{bmatrix}=2\Phi^{0}\begin{bmatrix}{\rho}_{0}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})x^{0}-(\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy})y^{0}+((x^{0})^{T}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})x^{0}){\Sigma}_{x}x^{0}\\\
{\rho}_{0}(\widehat{\Sigma}_{n,y}-{\Sigma}_{y})y^{0}-(\widehat{\Sigma}_{n,yx}-{\Sigma}_{yx})x^{0}+((y^{0})^{T}(\widehat{\Sigma}_{n,y}-{\Sigma}_{y})y^{0}){\Sigma}_{y}y^{0}\end{bmatrix}$
(14)

With this we are ready to state the main theorem of this paper.

######  Theorem 4.1 (Asymptotic representation of $\widehat{x}^{db}_{n}$).

Suppose $\mathcal{L}_{1}$ is as defined in (14), $\widehat{\alpha}_{n}$ and
$\widehat{\beta}_{n}$ satisfy Condition 4.1, and $\widehat{\Phi}_{n}$
satisfies Condition 4.1. Then under Assumption 2.1, 2.1, and 4.1, the
estimator $\widehat{x}^{db}_{n}$ defined in (9) can be expanded as either

$\quad\widehat{x}^{db}_{n}=x^{0}-\mathcal{L}_{(1)}+\text{rem},\quad\text{or}\quad\widehat{x}^{db}_{n}=-x^{0}-\mathcal{L}_{(1)}+\text{rem},$

where $\|\text{rem}\|_{\infty}=O_{p}(s^{2\kappa}\lambda^{2})$ with $s$ and
$\lambda$ as defined in Assumption 4.1 and (4), respectively.

A few remarks are in order about the statement and implications of Theorem
4.1. First, we not that under Assumption 4.1,
$\|\text{rem}\|_{\infty}=o_{p}(n^{-1/2})$. The importance of Theorem 4.1
subsequently lies in the fact that it establishes the equivalence between
$\widehat{x}_{n}^{db}$ and the more tractable random vector $\mathcal{L}$
under Assumption 4.1. In particular, one immediately can derive a simple yet
relevant corollary about the asymptotic normal nature of the distributions of
our de-biased estimators.

######  Corollary 4.2.

Under the set up of Theorem 4.1, for any $i=1,\ldots,p$, the following
assertions hold:

  * 1.

If $\alpha_{0,i}\neq 0$, then
$n^{1/2}\Big{(}(\widehat{x}^{db}_{n,i})^{2}-(x^{0}_{i})^{2}\Big{)}$ converges
in distribution to a centered Gaussian random variable with variance
$16\sigma_{i}^{2}(x^{0}_{i})^{2}$.

  * 2.

If $\alpha_{0,i}=0$, then $n(\widehat{x}^{db}_{n,i})^{2}$ converges in
distribution to a central Chi-squared random variable with degrees of freedom
one and scale parameter $4\sigma_{i}^{2}$.

The proof of Corollary 4.4 is deferred to the appendix. Before proceeding, it
is worth mentioning that the decision to provide inference on
$(\widehat{x}^{db}_{n,i})^{2}$ instead of $\widehat{x}^{db}_{n,i}$ is driven
by the fact that the former is unaffected by the sign flip of
$\widehat{x}^{db}_{n,i}$, which, unbeknown to us, can be centered at either
$x^{0}_{i}$ or $-x^{0}_{i}$. However, a result on $\widehat{x}^{db}_{n,i}$ can
also be derived under the set up of Theorem 4.1 and one has

$\sqrt{n}\Big{(}\widehat{x}^{db}_{n,i}-x^{0}_{i}\Big{)}\to_{d}N(0,4\sigma_{i}^{2})\quad\text{or}\quad\sqrt{n}\Big{(}\widehat{x}^{db}_{n,i}+x^{0}_{i}\Big{)}\to_{d}N(0,4\sigma_{i}^{2})\quad(i=1,\ldots,p).$
(15)

Moreover, we note that, often the inference on $(x^{0}_{i})^{2}$ suffices
since in practice the sign of $x^{0}_{i}$ is typically of little interest. As
a specific example, testing $H_{0}:x^{0}_{i}=0$, is equivalent to testing
$H_{0}:(x^{0}_{i})^{2}=0$. More importantly one of the central objects of
interest in low dimensional representations obtained through SCCA is the
projection operators onto the leading canonical directions. It is easy to see
that for this operator it is sufficiently to understand the squared
$(x^{0}_{i})^{2}$ and the cross terms $x^{0}_{i}x^{0}_{j}$ respectively. We
will also present asymptotic characterization of estimators for the cross-
terms $x^{0}_{i}x^{0}_{j}$. However, we first present a somewhat uniform
nature of the joint asymptotic normal behavior for the entire vector
$\hat{x}_{n}^{db}$. To this end, we verify in our next proposition that if
$\log p=o(n^{-1/7})$, then the convergence in (15) is uniform across
$i=1,\ldots,p$ while restricted to sets of suitably nice nature.

######  Proposition 4.3.

Let $\mathcal{A}_{p}$ be the set of all hyperrectangles in $\mathbb{R}^{p}$
and let $\Sigma_{p}$ the covariance matrix of the $p$-variate random vector
$(\mathcal{Z}(1),\ldots,\mathcal{Z}(p))$. Assume the set up of Theorem 4.1,
$\inf_{1\leq i\leq p}\sigma^{2}_{i}>c$ for some $c>0$, and that
$\log(p+q)=o(n^{-1/7})$. Then as $n\to\infty$, either

$\sup_{A\in\mathcal{A}_{p}}\bigg{|}P\Big{(}n^{1/2}(\widehat{x}^{db}_{n}-x^{0})\in
A\Big{)}-P\Big{(}2\mathbb{X}\in A\Big{)}\bigg{|}\to 0,$

or

$\sup_{A\in\mathcal{A}_{p}}\bigg{|}P\Big{(}n^{1/2}(\widehat{x}^{db}_{n}+x^{0})\in
A\Big{)}-P\Big{(}2\mathbb{X}\in A\Big{)}\bigg{|}\to 0,$

where $\mathbb{X}$ is a random vector distributed as $N_{p}(0,\Sigma_{p})$.

Proposition 4.3 can in turn be used, as promised earlier, to infer on the non-
diagonal elements of the matrix $x^{0}(x^{0})^{T}$. This is the content of our
next corollary – the proof of which can be found in Supplement 17.

######  Corollary 4.4.

Consider the set up of Proposition 4.3. Suppose $\Sigma_{p}$ is positive
definite. Let $i,j\in[p]$, and $i\neq j$. Denote by $\sigma_{ij}$ the
covariance between $\mathcal{Z}(i)$ and $\mathcal{Z}(j)$, where
$\mathcal{Z}(i)$’s are as defined in (4.2). Then the following assertions
hold:

  1. 1.

Suppose $x^{0}_{i}x^{0}_{j}\neq 0$. Then

$n^{1/2}\Big{(}(\widehat{x}_{n})_{i}(\widehat{x}_{n})_{j}-x^{0}_{i}x^{0}_{j}\Big{)}\to_{d}N\Big{(}0,4\\{(x^{0}_{i})^{2}\sigma^{2}_{j}+(x^{0}_{j})^{2}\sigma^{2}_{i}+2x^{0}_{i}x^{0}_{j}\sigma_{ij}\\}\Big{)}.$

  2. 2.

Suppose $x^{0}_{i}x^{0}_{j}=0$. Then

$n(\widehat{x}_{n})_{i}(\widehat{x}_{n})_{j}\to_{d}\mathbb{Z}_{i}\mathbb{Z}_{j},$

where $\mathbb{Z}_{i}\sim N(0,\sigma^{2}_{i})$, $\mathbb{Z}_{j}\sim
N(0,\sigma^{2}_{j})$, and
$\text{cov}(\mathbb{Z}_{i},\mathbb{Z}_{j})=\sigma_{ij}$.

Here once again we observe that the de-biased estimators of
$x^{0}_{i}x^{0}_{j}$ have different asymptotics depending on whether
$x^{0}_{i}x^{0}_{j}=0$ or not – which parallels the behavior of the de-biased
estimators of the diagonal elements we demonstrated earlier through Corollary
4.2.

######  Remark 4.5.

Proposition 4.3 can also used to simultaneously test the null hypotheses
$H_{0}:(\alpha_{0})_{i}=0$ $(i=1,\ldots,p)$. The uniform convergence in
Proposition 4.3 can be used o justify multiple hypothesis testing for the
coordinates of $x^{0}$ – whenever the coreesponding p-values are defined
through rectangular rejection regions based on $\hat{x}_{n}^{db}$. To this
end, one can use standard methods like Benjamini and Hochberg (BH) and
Benjamini and Yekutieli (BY) procedures for FDR control. The simultaneous
testing procedure can thereby also be connected to variable selection
procedures. However, we do not pursue it here since specialized methods are
available for the latter in SCCA context (Laha & Mukherjee, 2021).

The proof of Proposition 4.3, which can be found in Supplement 17, relies on a
Berry-Esseen type result. The lower bound requirement on the
$\sigma^{2}_{i}$’s is typical for such Berry Esseen type theorems – see
e.g.Chernozhukov et al. (2017). To check whether this assumptions actually can
hold in specific examples, we provide Corollary 4.6 below to establish the
validity of $\inf_{1\leq i\leq p}\sigma^{2}_{i}>c$ for some $c>0$ when $(X,Y)$
is jointly Gaussian. The proof of Corollary 4.6 can be found in Supplement 17.

######  Corollary 4.6.

Suppose $X,Y$ are jointly Gaussian and ${\rho}_{0}$ is bounded away from zero
and one. Further suppose $\log(p+q)=o(n^{-1/7})$. Then under the set up of
Theorem 4.3, the assertion of $\inf_{1\leq i\leq p}\sigma^{2}_{i}>c$ for some
$c>0$ used in Proposition 4.3 holds.

We end our discussions regarding the inference of $x^{0}$ with a method for
consistent estimation of the $\sigma_{i}^{2}$’s. indeed, this will allow us to
develop tests for the hypotheses $H_{0}:x_{i}^{0}=0$ or build confidence
interval for $(x_{i}^{0})^{2}$. To this end we partition
$(\widehat{\Phi}_{n})_{i}=(\widehat{\Phi}_{i,1},\widehat{\Phi}_{i,2})$ where
$\widehat{\Phi}_{i,1}\in\mathbb{R}^{p}$, and
$\widehat{\Phi}_{i,2}\in\mathbb{R}^{q}$. Because
$\sigma_{i}^{2}=\text{var}(\mathcal{Z})$, for $i=1,\ldots,p$, it can be shown
that a consistent estimator is given by the variance of pseudo-observations
$\\{\widehat{Z}_{j}(i)\\}_{j=1}^{n}$ $(i=1,\ldots,p+q)$, which are defined by

$\displaystyle\widehat{Z}_{j}(i)=$ $\displaystyle\
[\widehat{\rho}_{n}\widehat{\Phi}_{i,1}^{T}+\\{\widehat{\Phi}_{i,1}^{T}\widehat{\Sigma}_{n,x}\widehat{x}_{n}\\}(\widehat{x}_{n})^{T}]X_{j}X_{j}^{T}\widehat{x}_{n}-\widehat{\Phi}_{i,1}^{T}X_{j}Y_{j}^{T}\widehat{y}_{n}$
$\displaystyle\
+[\widehat{\rho}_{n}\widehat{\Phi}_{i,2}^{T}+\\{\widehat{\Phi}_{i,2}^{T}\widehat{\Sigma}_{n,y}\widehat{y}_{n}\\}(\widehat{y}_{n})^{T}]Y_{j}Y_{j}^{T}\widehat{y}_{n}-(\widehat{x}_{n})^{T}X_{j}Y_{j}^{T}\widehat{\Phi}_{i,2}.$

Our final result pertains to the asymptotic distribution of
$\widehat{\rho}_{n}^{2,db}$.

######  Theorem 4.7.

Suppose $s^{2\kappa+1/2}\lambda^{2}=o(n^{-1/2})$ and ${\rho}_{0}<1$. Then
under the set-up of Theorem 4.1,

$n^{1/2}(\widehat{\rho}_{n}^{2,db}-{\rho}_{0}^{2})\to_{d}N(0,\sigma^{2}_{\rho}),$

where
$\sigma^{2}_{\rho}=\text{var}({\rho}_{0}(X^{T}x^{0})^{2}+{\rho}_{0}(Y^{T}y^{0})^{2}-2(X^{T}x^{0})(Y^{T}y^{0}))$.
In particular, when the observations are Gaussian,
$\sigma^{2}_{\rho}={\rho}_{0}^{2}(1-{\rho}_{0}^{2})^{2}$.

A few remarks are in order regarding content of Theorem 4.7. First, one can
$\sigma_{\rho}^{2}$ is consistently

$\widehat{\sigma}^{2}_{\rho}=\sum_{j=1}^{n}\frac{(\widehat{x}_{n}^{T}X_{j})^{2}(\widehat{y}_{n}^{T}Y_{j})^{2}}{n}-\widehat{\rho}_{n}^{4},$

and thereby use Theorem 4.7 to create asymptotically valid confidence
intervals for leading canonical signal strength. Further note that Theorem 4.7
requires stricter condition on $s$ compared to Theorem 4.1. Although we have
not explored the sharpness of this assumption, one can find similar stricter
sparsity requirement in Janková & van de Geer (2018) while demonstrating
$n^{1/2}$-consistency of a de-biased estimator for the largest eigenvalue in
the sparse PCA problem. Finally, the value of $\sigma^{2}_{\rho}$ in the
Gaussian case matches that of the parametric MLE of ${\rho}_{0}^{2}$ under the
Gaussian model (Anderson, 2003, p.505). Such agreement is generally observed
in case of the de-biased estimators, e.g. the de-biased estimator of the
principal eigenvalue (Janková & van de Geer, 2018).

##  5 Numerical Experiments

###  5.1 Preliminaries

In this section we explore aspects of finite sample behavior of the methods
discusses in earlier sections. Further numerical experiments are collected in
Supplement 8.1 where we compare the bias of our method to popular SCCA
alternatives. We start with some preliminary discussions on the choice for the
set-up, initial estimators, and tuning parameters.

Set Up: The set-ups under which we will conduct our comparisons can be
described through specifying the nuisance parameters (marginal covariance
matrices of $X$ and $Y$) along with the strength ($\rho$), sparsity, rank of
$\Sigma_{xy}$, and the joint distribution of $X,Y$. For the marginal marginal
covariance matrices of $X$ and $Y$, motivated by previously studied cases in
the literature (Mai & Zhang, 2019; Gao et al., 2017) we shall consider two
cases as follows:

Identity.

This will correspond to the case where ${\Sigma}_{x}={\Sigma}_{y}=I_{p}$

Sparse-inverse.

This will correspond to the case where ${\Sigma}_{x}={\Sigma}_{y}$ is the
correlation matrix obtained from $\Sigma_{0}$, where $\Sigma_{0}=\Omega^{-1}$,
and $\Omega$ is a sparse matrix with the form

$\Omega_{ij}=1_{\\{i=j\\}}+0.5\times 1_{\\{|i-j|=1\\}}+0.4\times
1_{\\{|i-j|=2\\}},\quad i,j\in[p].$

Analogous to Mai & Zhang (2019) and Gao et al. (2017), we shall also take
${\Sigma}_{xy}=\rho_{0}{\Sigma}_{x}\alpha_{0}\beta_{0}^{T}{\Sigma}_{y}$ to be
a rank one matrix, where we consider the canonical vectors $\alpha_{0}$ and
$\beta_{0}$ with sparsity $2$ as follows:

$\displaystyle\alpha_{*}=(1,1,0,\ldots,0)^{T},$
$\displaystyle\quad\beta_{*}=(1,1,0,\ldots,0)^{T},\quad\alpha_{0}=\frac{\alpha_{*}}{\sqrt{\alpha_{*}^{T}\Sigma_{x}\alpha_{*}}},$
$\displaystyle\quad\beta_{0}=\frac{\beta_{*}}{\sqrt{\beta_{*}^{T}\Sigma_{y}\beta_{*}}}.$

The canonical correlation ${\rho}_{0}$ depicts the signal strength in our set
up. We will explore three different values for the $\rho_{0}$: 0.2, 0.5, and
0.9, which will be referred as the small, medium, and the high signal strength
settings, respectively. The joint distribution of $X,Y$ is finally taken to be
Gaussian with mean $0$. Also, throughout we set the $(p,q,n)$ combination to
be $(80,80,500)$, $(300,200,500)$, and $(600,200,500)$, which correspond to
$p+q$ being small, moderate, and moderately high, respectively. Finally, we
will always consider $N=1000$ Monte Carlo samples.

Initial Estimators and Tuning Parameters: We construct the preliminary
estimators using the modified COLAR algorithm (see Algorithm 1). For the rank
one case, the latter coincides with Gao et al. (2017)’s COLAR estimator.
Recall that throughout we set the $(p,q,n)$ combination to be $(80,80,500)$,
$(300,200,500)$, and $(600,200,500)$. One of the reasons we do not accommodate
higher $p$ and $q$ because the COLAR algorithm, as it is, does not scale well
with $p$ and $q$ 222This was also noted by Mai & Zhang (2019).  Also, we do
not consider smaller values of $n$ since it is expected that de-biasing
procedures generally require $n$ to be at least moderately large (see e.g.
Janková & van de Geer (2018)).

In our proposed methods, tuning parameters arise from two sources: (a)
estimation of the preliminary estimators and (b) precision matrix estimation.
To implement the modified COLAR algorithm, we mostly follow the code for COLAR
provided by the authors Gao et al. (2017). The COLAR penalty parameters,
$\lambda_{1}$ and $\lambda_{2}$, were left as specified in the COLAR code,
namely $\lambda_{1}=0.55\\{\log(p)/n\\}^{1/2}$ and
$\lambda_{2}=[\\{1+\log(p)\\}/n]^{1/2}$. The tolerance level was fixed at
$10^{-4}$ with a fixed maximum of 200 iterations for the first step of the
COLAR algorithm. Next consider the tuning strategy for the nodewise lasso
algorithm (Algorithm 2), which involves the lasso penalty parameter
$\lambda_{j}^{nl}$ and the parameter $B_{j}$ ($j=1,\ldots,p+q$). Theorem 10.3
proposes the choice $\lambda_{j}^{nl}=C\cdot\sqrt{{\log(p+q)}/{n}}$ for all
$j\in[p+q]$. In our simulations, the parameter $C$ is empirically determined
to minimize
$|\widehat{\Phi}_{n}\widehat{H}(\widehat{x}_{n},\widehat{y}_{n})-I_{p+q}|_{\infty}$.
For the settings $(80,80,500)$ and $(300,200,500)$, this parameter is set at
$40$ and $50$ for the identity and sparse inverse cases, respectively. For the
moderately high $p+q$ setting, this parameter is set at $20$. The nodewise
lasso parameter $B_{j}$ is taken to be $10/\lambda_{j}$, which is in line with
Janková & van de Geer (2018), who recommends taking $B_{j}\approx
1/\lambda_{j}$.

Targets of Inference: We present our results for the $1^{st}$ and the
$20^{\text{nd}}$ element of $x^{0}$. The former stands for a typical non-zero
element, where the latter represents a typical zero element. For each element,
we compute confidence intervals for $(x^{0}_{i})^{2}$, and test the null
$H_{0}:$ $|x^{0}_{i}|=0$ $(i=1,20)$. For the latter, we use a
$\chi^{2}$-squared test based on the asymptotic null distribution of
$(\widehat{x}^{db}_{n})^{2}$ given in part two of Corollary 4.2. As mentioned
earlier, this test is equivalent to testing $H_{0}:(\alpha_{0})_{i}=0$. The
construction of the confidence intervals, which we discuss next, is a little
more subtle.

We construct two types of confidence interval. For any $i\in[p]$, the first
confidence interval, which will be referred as the ordinary interval from now
on, is given by

$\Big{(}\max\\{0,(\widehat{x}^{db}_{n,i})^{2}-l_{\text{CI},i}\\},(\widehat{x}^{db}_{n,i})^{2}+l_{\text{CI},i}\Big{)},\quad\text{where}\quad
l_{\text{CI},i}=4z_{0.975}|\widehat{x}^{db}_{n,i}|\widehat{\sigma}_{i}/\sqrt{n}.$
(16)

Here $z_{0.975}$ is the $0.975^{\text{th}}$ quantile of the standard Gaussian
distribution. Corollary 4.2 shows that the asymptotic coverage of the above
confidence interval is $95\%$ when $x^{0}_{i}\neq 0$. For $x_{i}^{0}=0$,
however, the above confidence interval can have asymptotic coverage higher
than $0.95\%$. To see why, note that $(\widehat{x}^{db}_{n,i})^{2}=O_{p}(1/n)$
by Corollary 4.2 in this case. Since both the length and the center of the
ordinary interval depends on $(\widehat{x}^{db}_{n,i})^{2}$, the coverage can
suffer greatly if $(\widehat{x}^{db}_{n,i})^{2}$ underestimates
$(x^{0}_{i})^{2}$. Therefore, we construct another confidence interval by
relaxing the length of the ordinary intervals. This second interval, to be
referred as the conservative interval from now on, is obtained by simply
substituting the $\widehat{x}^{db}_{n,i}$ in the standard deviation term
$l_{\text{CI}}$ in (16) by $\max(|\widehat{x}^{db}_{n,i}|,1)$. Clearly, the
conservative interval can have potentially higher coverage than $95\%$, which
motivates our nomenclature.

###  5.2 Results

We divide the presentation of our results on coordinates with and without
signal, followed by discussions about issues regarding distinctions between
asymptotic and finite sample considerations of our method.

Inference when there is no signal: If $x^{0}_{i}=0$, both confidence intervals
(CI) exhibit high coverage, often exceeding $95\%$, across all settings; see
Figures [$x^{0}_{20}$ plots] in Supplement 8. This is unsurprising in view of
the discussion in the previous paragraph. The conservative confidence
intervals have substantially larger length, which is understandable because
the ratio between the ordinary and the conservative CI length is $O_{p}(1/n)$
in this case. Also, the length of the confidence intervals generally decrease
as the signal strength increases, as expected. The rejection frequency of the
tests (the type I error in this scenario), generally stays below $0.05$,
especially at medium to high signal strength.

Inference when there is signal: When $x^{0}_{i}\neq 0$, the ordinary intervals
exhibit poor coverage at the low and medium signal strength regardless of the
underlying covariance matrix structure, although the performance seems to be
worse for sparse inverse matrices. Figure 1 entails that this underperformance
is due to the underestimation of small signals $(\widehat{x}^{db}_{1})^{2}$,
which is tied to the high negative bias of the preliminary estimator in these
cases; see the histograms in Figure 5. This issue will be discussed in more
detail in Supplement 8.1. Figure 1 also implies that if $(x^{0}_{i})^{2}$ is
small, the confidence intervals crowd near the origin. Also at the high signal
strength, the coverage of the ordinary intervals fail to reach the desired
$95\%$ level.

The relaxation of the ordinary confidence interval length, which leads to the
conservative intervals, substantially improve the coverage, with the
improvement being dramatic at low signal. In the latter case, the conservative
intervals enjoy high coverage, which is well over $95\%$ for moderate or
higher $p,q$. In this case, in general, the relaxation results in a four-fold
or higher increase in the confidence interval length. As signal strength
increases, the increase in the confidence interval length gets smaller, and
consequently, the increase in the coverage slows down. This is unsurprising
noting the ratio between the length of the conservative and the ordinary
interval is proportional to $\widehat{\rho}_{n}^{-1}$. One should be cautious
with the relaxation, however, because it may lead to inclusion of not only the
true signal, as desired, but also zero. This can be clearly seen in the medium
signal strength case of the sparse inverse matrix; compare the middle column
of Figure 1 (b) with that of Figure 2 (b). The inclusion of origin does not
bring any advantage for the relaxed intervals in the no-signal case either,
because as discussed earlier, in the latter case the ordinary intervals are
themselves efficient, with the relaxed versions hardly making any improvement.

Discussion on Asymptotics: The performance of the confidence intervals improve
if $(n,p,q)$ increase. See for example the illustration in Figure 6 in
Supplement 8.2 where the triplet has been doubled. Interestingly, the
asymptotics successfully kicks in for the corresponding tests as soon as the
signal strength reaches the medium level. The test attains power higher than
$0.673$ at the medium signal strength, and the perfect power of one at high
signal strength. This phenomenon is the result of the super-efficiency of the
de-biased estimator at $x^{0}_{i}=0$, as elicited by Corollary 4.2. Since the
test exploits the knowledge of this faster convergence under the null, it has
better precision than the confidence interval, which is oblivious to this
fact. In many situations, the test may get rejected but the confidence
intervals, even the ordinary one, may include zero. During implementation, if
one faces such a situation, they should conclude that either the signal
strength is too small or the sample size is not sufficient for the confidence
intervals to be too precise.

Discussions on Performance of De-biased SCCA: We conclude that since the de-
biased estimators work on sparse estimators which are super efficient at zero,
the inference does not face any obstacle if the true signal $x^{0}_{i}=0$. In
presence of signal, the tests are generally reliable if the signal strength is
at least moderate. In contrast, the ordinary confidence intervals, which are
blindly based on Corollary 4.2, struggle whenever the initial COLAR estimators
incur a bias too large for the de-biasing step to overcome. This is generally
observed at low to medium signal strength. The conservative intervals can
solve this problem partially at the cost of increased length. At present, the
$l_{1}$ and $l_{2}$ guarantees as required by Condition 4.1 are only available
for COLAR type estimators. The performance of the ordinary confidence
intervals may improve if one can construct a SCCA preliminary estimator with
similar strong theoretical garuantees, but better empirical performance in
picking up small signal. Searching for a different SCCA preliminary estimator
is important for another reason – COLAR is not scalable to ultra high
dimension. This problem occurs because COLAR relies on semidefinite
programming, whose scalability issues are well noted (Dey et al., 2018).

(a) Ordinary confidence intervals for identity matrix

(b) Ordinary confidence intervals for sparse inverse matrix

Figure 1: Ordinary confidence intervals for $(x^{0}_{1})^{2}$

(a) Conservative confidence intervals for identity matrix

(b) Conservative confidence intervals for sparse inverse matrix

Figure 2: Conservative confidence intervals for $(x^{0}_{1})^{2}$

##  6 Real Data Application

The physiological functions in human bodies are controlled by complex
pathways, whose deregulation lead to myriad diseases. Therefore it is
important to understand the interaction between different factors
participating in these biological pathways, such as proteins, genes etc. We
consider two important pathways: (a) Cytokine-cytokine receptor interaction
pathway and (b) Adipocytokine signalling pathway. Cytokines are released in
response to inflammation in the body, and pathway (a) is thus related to viral
infection, cell-growth, differentiation, and cancer progression (Lee & Rhee,
2017). Pathway (b) is involved in fat metabolism and insulin resistance, thus
playing a vital role in diabetes (Pittas et al., 2004). We wish to study the
linear interaction between the group of genes and proteins that are involved
in these pathways. To that end, we use the Microarray and proteomic datasets
analysed by Lee et al. (2011), which are originally from the National Cancer
Institute, and available at http://discover.nci.nih.gov/cellminer/.

The dataset contains sixty human cancer cell lines. We use $59$ of the sixty
observations because one has missing microarray information. Although the
microarray data has information on many genes, we considered only those
involved in pathways (a) and (b), giving $p=230$ and $62$ miRNAs,
respectively. To this end, we use https://www.genscript.com/ to get the list
of genes participating in these pathways. The dataset contains $q=94$
proteins. We center and scale all variables prior to our analysis.

Figure 12 indicates that most genes and proteins have negligible correlation,
which hints that only a handful of genes and proteins share linear
interactions in the pathways under concern – thus supporting the possibility
of $\alpha_{0}$ and $\beta_{0}$ being low dimensional. On the other hand,
Figure 11 hints at the existence of low dimensional structures in the variance
matrices of both the genes and the proteins. However, it seems unlikely that
they are totally uncorrelated among themselves, which questions the
applicability of popular methods only suited for diagonal variance matrices,
e.g. PMA (Witten et al., 2009).

Apart from the de-biased estimators, we also look into the SCCA estimates of
the leading canonical covariates using Mai & Zhang (2019), Gao et al. (2017),
Witten et al. (2009), and Wilms & Croux (2015)’s methods. The first three
methods were implemented as discussed in Supplement 8.1. To apply Wilms &
Croux (2015)’s methods, we used the code provided by the authors with the
default choice of tuning parameters. Among these methods, only Witten et al.
(2009)’s method requires ${\Sigma}_{x}$ and ${\Sigma}_{y}$ to be diagonal. For
these methods, we say a gene or protein is “detected” if the corresponding
loading, i.e. the estimated $(\widehat{\alpha}_{n})_{i}$ or
$(\widehat{\beta}_{n})_{i}$, is nonzero.

We construct confidence intervals, both ordinary and conservative, and test
the null that $x^{0}_{i}=0$ or $y^{0}_{j}=0$ for each $i\in[p]$ and $j\in[q]$,
as discussed in Section 5. We apply the false discovery rate corrections of
Benjamini and Hochberg (BH) as well as Benjamini and Yekutieli (BY), the
latter of which does not assume independent P-values. Table 1 tabulate the
number of detections by the above-mentioned methods. Even after false
discovery rate adjustment, most discoveries seem to include zero in the
confidence intervals. We discussed this situation in Section 5, where it was
indicated that the former can occur if the signal strength is small or the
sample size is insufficient. To be conservative, we consider only those genes
and proteins whose ordinary interval excludes zero. These discoveries are
reported in Tables 2 and 3 along with the confidence intervals. The pictorial
representation of the confidence intervals can be found in Figure 9 and
Figures 10 in Supplement 8.2.

Using Gene Ontology toolkit available at _http://geneontology.org/_ , we
observe that our discovered from pathway (a) are mainly involved in biological
processes like positive regulation of gliogenesis and molecular function like
growth factor activity, where the selected proteins play a role in regulating
membrane assembly, enzyme function, and other cellular functions. Gene
Ontology toolkit also entails that the discovered genes from pathway (b) are
involved in positive regulation of cellular processes, and molecular function
like growth factor activity. The only discovered gene in pathway (b) is ANXA2,
which, according to UNIPORT at _https://www.uniport.org_ , is a membrane-
binding protein involved in RNA binding and host-virus infection.

Variable | Mai & Zhang | Wilms & Croux | Gao et al. | Witten et al. | DB+BH | DB+BY
---|---|---|---|---|---|---
| Pathway (a)
Genes | 2 (2) | 1 (1) | 3 (3) | 41 (5) | 13 | 6
Proteins | 4 (3) | 1 (1) | 7 (5) | 13 (5) | 36 | 22
| Pathway (b)
Genes | 2 (1) | 1 (1) | 4 (3) | 11 (2) | 8 | 5
Proteins | 7 (1) | 1 (1) | 9 (1) | 12 (1) | 22 | 2
Table 1: Number of detections: number of non-zero loadings in different SCCA estimators and number of detections by our tests (DB) after Benjamini and Hochberg (BH) and Benjamini and Yekutieli (BY) false discovery rate correction. For the SCCA estimators, size of their intersection with DB+BY are given in parentheses. Gene |  $p$-value* | 95% CI | Relaxed CI | Discovered by
---|---|---|---|---
CLCF1 | 2.0E-07 | (0.055, 0.39) | (0, 0.58) | Witten et al.
EGFR | 8.8E-09 | (0.11, 0.58) | (0,0.74) |  Mai & Zhang, Witten et al.,
|  |  |  | Gao et al.
LIF | 1.6E-05 | (0.022, 0.45) | (0, 0.68) |  Witten et al., Gao et al.
PDGFC | 1.4E-07 | (0.094, 0.64) | (0, 0.82) | Witten et al.
TNFRSF12A | 7.8E-11 | (0.15, 0.60) | (0.01, 0.75) |  Mai & Zhang, Witten et al.,
|  |  |  |  Gao et al., Wilms & Croux
Protein |  $p$-value* | 95% CI | Relaxed CI | Discovered by
ANXA2 | 1.3E-15 | (0.13, 0.38) | (0.01, 0.51) |  Mai & Zhang, Witten et al.,
|  |  |  |  Gao et al., Wilms & Croux
CDH2 | 5.1E-09 | (0.22, 1.1) | (0.12, 1.23) |  Mai & Zhang, Witten et al.,
|  |  |  | Gao et al.
FN1 | 4.2E-07 | (0.96, 7.6) | (0.96, 7.6) | none
GTF2B | 6.7E-05 | (0.034, 4.0) | (0.034, 4.0) | none
KRT20 | 1.2E-05 | (0.015, 0.27) | (0, 0.48) | none
MVP | 2.6E-05 | (0.021, 0.59) | (0, 0.82) | Witten et al.
Table 2: Discovered genes and protein from pathway (a). The confidence
intervals are obtained using the methods described in Section 5. The P-values
are the original P-values before false discovery rate correction.
*All genes and proteins were also detected by Benjamini and Yekutieli method. Gene |  $p$-value* | 95% CI | Relaxed CI | Discovered by
---|---|---|---|---
ACSL5 | 2.9E-05 | (0.014, 0.45) | (0, 0.68) | none
RXRG | 4.1E-10 | (0.073, 0.32) | (0, 0.47) |  Wilms & Croux, Gao et al.,
|  |  |  | Mai & Zhang
TNFRSF1B | 1.1E-09 | (0.49, 2.2) | (0.49, 2.2) | none
Protein |  $p$-value* | 95% CI | Relaxed CI | Discovered by
ANXA2 | 2.7E-74 | (1.1, 1.7) | (1.1, 1.7) | none
Table 3: Discovered genes and protein from pathway (b). The confidence
intervals are obtained using the methods described in Section 5. The P-values
are the original P-values before false discovery rate correction.
*All genes and proteins were also detected by Benjamini and Yekutieli method.

##  7 Acknowledgements

Rajarshi Mukherjee’s research was partially supported by NSF Grant
EAGER-1941419 and NIH Grant NIH/NIEHS P42ES030990.

## References

  * Anderson (2003) Anderson, T. (2003).  An Introduction to Multivariate Statistical Analysis.  Wiley Series in Probability and Statistics. Wiley.
  * Anderson (1962) Anderson, T. W. (1962).  An introduction to multivariate statistical analysis.  Tech. rep., Wiley New York.
  * Bao et al. (2019) Bao, Z., Hu, J., Pan, G. & Zhou, W. (2019).  Canonical correlation coefficients of high-dimensional gaussian vectors: Finite rank case.  Ann. Statist. 47, 612–640.
  * Bellec & Zhang (2019) Bellec, P. C. & Zhang, C.-H. (2019).  De-biasing the lasso with degrees-of-freedom adjustment.  arXiv preprint arXiv:1902.08885 .
  * Billingsley (2008) Billingsley, P. (2008).  Probability and measure.  John Wiley & Sons.
  * Boyd et al. (2004) Boyd, S., Boyd, S. P. & Vandenberghe, L. (2004).  Convex optimization.  Cambridge university press.
  * Bühlmann & Van De Geer (2011) Bühlmann, P. & Van De Geer, S. (2011).  Statistics for high-dimensional data: methods, theory and applications.  Springer Science & Business Media.
  * Cai et al. (2011) Cai, T., Liu, W. & Luo, X. (2011).  A constrained $l_{1}$ minimization approach to sparse precision matrix estimation.  Journal of the American Statistical Association 106, 594–607.
  * Cai et al. (2017) Cai, T. T., Guo, Z. et al. (2017).  Confidence intervals for high-dimensional linear regression: Minimax rates and adaptivity.  The Annals of statistics 45, 615–646.
  * Cai et al. (2018) Cai, T. T., Zhang, A. et al. (2018).  Rate-optimal perturbation bounds for singular subspaces with applications to high-dimensional statistics.  The Annals of Statistics 46, 60–89.
  * Chatterjee & Lahiri (2010) Chatterjee, A. & Lahiri, S. (2010).  Asymptotic properties of the residual bootstrap for lasso estimators.  Proceedings of the American Mathematical Society 138, 4497–4509.
  * Chatterjee & Lahiri (2011) Chatterjee, A. & Lahiri, S. N. (2011).  Bootstrapping lasso estimators.  Journal of the American Statistical Association 106, 608–625.
  * Chatterjee & Lahiri (2013) Chatterjee, A. & Lahiri, S. N. (2013).  Rates of convergence of the adaptive lasso estimators to the oracle distribution and higher order refinements by the bootstrap.  The Annals of Statistics 41, 1232–1259.
  * Chen et al. (2013) Chen, M., Gao, C., Ren, Z. & Zhou, H. H. (2013).  Sparse cca via precision adjusted iterative thresholding.  arXiv preprint arXiv:1311.6186 .
  * Chen et al. (2012) Chen, X., Han, L. & Carbonell, J. (2012).  Structured sparse canonical correlation analysis.  In Artificial intelligence and statistics. PMLR.
  * Chen et al. (2020) Chen, Y., Chi, Y., Fan, J. & Ma, C. (2020).  Spectral methods for data science: A statistical perspective.  arXiv preprint arXiv:2012.08496 .
  * Chernozhukov et al. (2017) Chernozhukov, V., Chetverikov, D., Kato, K. et al. (2017).  Central limit theorems and bootstrap in high dimensions.  Annals of Probability 45, 2309–2352.
  * Dey et al. (2018) Dey, S. S., Mazumder, R. & Wang, G. (2018).  A convex integer programming approach for optimal sparse pca.  arXiv preprint arXiv:1810.09062 .
  * Eckart & Young (1936) Eckart, C. & Young, G. (1936).  The approximation of one matrix by another of lower rank.  Psychometrika 1, 211–218.
  * Friedman et al. (2008) Friedman, J., Hastie, T. & Tibshirani, R. (2008).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9, 432–441.
  * Gao et al. (2015) Gao, C., Ma, Z., Ren, Z., Zhou, H. H. et al. (2015).  Minimax estimation in sparse canonical correlation analysis.  The Annals of Statistics 43, 2168–2197.
  * Gao et al. (2017) Gao, C., Ma, Z., Zhou, H. H. et al. (2017).  Sparse cca: Adaptive estimation and computational barriers.  The Annals of Statistics 45, 2074–2101.
  * Holm et al. (2010) Holm, K., Hegardt, C., Staaf, J., Vallon-Christersson, J., Jönsson, G., Olsson, H., Borg, Å. & Ringnér, M. (2010).  Molecular subtypes of breast cancer are associated with characteristic dna methylation patterns.  Breast cancer research 12, 1–16.
  * Hotelling (1992) Hotelling, H. (1992).  Relations between two sets of variates.  In Breakthroughs in statistics. Springer, pp. 162–190.
  * Hu et al. (2016) Hu, W., Lin, D., Calhoun, V. D. & Wang, Y.-p. (2016).  Integration of snps-fmri-methylation data with sparse multi-cca for schizophrenia study.  In 2016 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE.
  * Hu et al. (2017) Hu, W., Lin, D., Cao, S., Liu, J., Chen, J., Calhoun, V. D. & Wang, Y.-P. (2017).  Adaptive sparse multiple canonical correlation analysis with application to imaging (epi) genomics study of schizophrenia.  IEEE Transactions on Biomedical Engineering 65, 390–399.
  * Janková & Van De Geer (2016) Janková, J. & Van De Geer, S. (2016).  Confidence regions for high-dimensional generalized linear models under sparsity.  arXiv preprint arXiv:1610.01353 .
  * Janková & van de Geer (2017) Janková, J. & van de Geer, S. (2017).  Honest confidence regions and optimality in high-dimensional precision matrix estimation.  Test 26, 143–162.
  * Janková & van de Geer (2018) Janková, J. & van de Geer, S. (2018).  De-biased sparse pca: Inference and testing for eigenstructure of large covariance matrices.  arXiv preprint arXiv:1801.10567 .
  * Javanmard & Montanari (2014) Javanmard, A. & Montanari, A. (2014).  Confidence intervals and hypothesis testing for high-dimensional regression.  The Journal of Machine Learning Research 15, 2869–2909.
  * Kang et al. (2013) Kang, M., Zhang, B., Wu, X., Liu, C. & Gao, J. (2013).  Sparse generalized canonical correlation analysis for biological model integration: a genetic study of psychiatric disorders.  In 2013 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). IEEE.
  * Laha & Mukherjee (2021) Laha, N. & Mukherjee, R. (2021).  On support recovery with sparse cca: Information theoretic and computational limits.  arXiv preprint arXiv:2108.06463 .
  * Lee & Rhee (2017) Lee, M. & Rhee, I. (2017).  Cytokine signaling in tumor progression.  Immune network 17, 214.
  * Lee et al. (2011) Lee, W., Lee, D., Lee, Y. & Pawitan, Y. (2011).  Sparse canonical covariance analysis for high-throughput data.  Statistical Applications in Genetics and Molecular Biology 10.
  * Leeb & Pötscher (2005) Leeb, H. & Pötscher, B. M. (2005).  Model selection and inference: Facts and fiction.  Econometric Theory 21, 21–59.
  * Leeb & Pötscher (2006) Leeb, H. & Pötscher, B. M. (2006).  Can one estimate the conditional distribution of post-model-selection estimators?  The Annals of Statistics 34, 2554–2591.
  * Leeb & Pötscher (2008) Leeb, H. & Pötscher, B. M. (2008).  Sparse estimators and the oracle property, or the return of hodges’ estimator.  Journal of Econometrics 142, 201–211.
  * Lin et al. (2013) Lin, D., Zhang, J., Li, J., Calhoun, V. D., Deng, H.-W. & Wang, Y.-P. (2013).  Group sparse canonical correlation analysis for genomic data integration.  BMC bioinformatics 14, 1–16.
  * Ma et al. (2020) Ma, Z., Li, X. et al. (2020).  Subspace perspective on canonical correlation analysis: Dimension reduction and minimax rates.  Bernoulli 26, 432–470.
  * Ma et al. (2013) Ma, Z. et al. (2013).  Sparse principal component analysis and iterative thresholding.  The Annals of Statistics 41, 772–801.
  * Mai & Zhang (2019) Mai, Q. & Zhang, X. (2019).  An iterative penalized least squares approach to sparse canonical correlation analysis.  Biometrics .
  * Mazumder & Hastie (2012) Mazumder, R. & Hastie, T. (2012).  The graphical lasso: New insights and alternatives.  Electronic journal of statistics 6, 2125.
  * Meinshausen et al. (2006) Meinshausen, N., Bühlmann, P. et al. (2006).  High-dimensional graphs and variable selection with the lasso.  The annals of statistics 34, 1436–1462.
  * Mitra et al. (2016) Mitra, R., Zhang, C.-H. et al. (2016).  The benefit of group sparsity in group inference with de-biased scaled group lasso.  Electronic Journal of Statistics 10, 1829–1873.
  * Neykov et al. (2018) Neykov, M., Ning, Y., Liu, J. S., Liu, H. et al. (2018).  A unified theory of confidence regions and testing for high-dimensional estimating equations.  Statistical Science 33, 427–443.
  * Ning et al. (2017) Ning, Y., Liu, H. et al. (2017).  A general theory of hypothesis tests and confidence regions for sparse high dimensional models.  Annals of statistics 45, 158–195.
  * Pittas et al. (2004) Pittas, A. G., Joseph, N. A. & Greenberg, A. S. (2004).  Adipocytokines and insulin resistance.  The Journal of Clinical Endocrinology & Metabolism 89, 447–452.
  * Pötscher & Leeb (2009) Pötscher, B. M. & Leeb, H. (2009).  On the distribution of penalized maximum likelihood estimators: The lasso, scad, and thresholding.  Journal of Multivariate Analysis 100, 2065–2082.
  * Rao & Bhimasankaram (2000) Rao, A. & Bhimasankaram, P. (2000).  Linear Algebra.  Texts and Readings in Mathematics. Hindustan Book Agency.
  * Sofer et al. (2012) Sofer, T., Maity, A., Coull, B., Baccarelli, A. A., Schwartz, J. & Lin, X. (2012).  Multivariate gene selection and testing in studying the exposure effects on a gene set.  Statistics in biosciences 4, 319–338.
  * Tsiatis (2007) Tsiatis, A. (2007).  Semiparametric theory and missing data.  Springer Science & Business Media.
  * van de Geer et al. (2014) van de Geer, S., Bühlmann, P., Ritov, Y. & Dezeure, R. (2014).  On asymptotically optimal confidence regions and tests for high-dimensional models.  Ann. Statist. 42, 1166–1202.
  * Vershynin (2010) Vershynin, R. (2010).  Introduction to the non-asymptotic analysis of random matrices.  arXiv preprint arXiv:1011.3027 .
  * Vershynin (2018) Vershynin, R. (2018).  High-dimensional probability: An introduction with applications in data science, vol. 47.  Cambridge university press.
  * Wang et al. (2016) Wang, T., Berthet, Q., Samworth, R. J. et al. (2016).  Statistical and computational trade-offs in estimation of sparse principal components.  The Annals of Statistics 44, 1896–1930.
  * Wilms & Croux (2015) Wilms, I. & Croux, C. (2015).  Sparse canonical correlation analysis from a predictive point of view.  Biometrical Journal 57, 834–851.
  * Witten et al. (2009) Witten, D. M., Tibshirani, R. & Hastie, T. (2009).  A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis.  Biostatistics 10, 515–534.
  * Yu et al. (2015) Yu, Y., Wang, T. & Samworth, R. J. (2015).  A useful variant of the davis–kahan theorem for statisticians.  Biometrika 102, 315–323.
  * Yuan (2010) Yuan, M. (2010).  High dimensional inverse covariance matrix estimation via linear programming.  Journal of Machine Learning Research 11, 2261–2286.
  * Yuan & Zhang (2013) Yuan, X.-T. & Zhang, T. (2013).  Truncated power method for sparse eigenvalue problems.  Journal of Machine Learning Research 14, 899–925.
  * Zhang & Zhang (2014) Zhang, C.-H. & Zhang, S. S. (2014).  Confidence intervals for low dimensional parameters in high dimensional linear models.  Journal of the Royal Statistical Society: Series B: Statistical Methodology , 217–242.

Supplement to “On Statistical Inference with High Dimensional Sparse CCA”

##  8  Extra Simulations

###  8.1 Bias estimation

This section compares the elementwise bias of our de-biased CCA estimator with
other commonly used sparse CCA estimators. We use the same simulation settings
as in Section 5. Also, the tuning parameters for the de-biased estimators are
kept exactly as in Section 5. As competitors, we choose the COLAR estimator of
Gao et al. (2017), and the SCCA of Mai & Zhang (2019) and Witten et al.
(2009). Since we are in the rank one setting, the COLAR estimator coincides
with the modified COLAR estimator, which is our preliminary estimator, and has
already been discussed in Section 5. The SCCA of Mai & Zhang (2019) is
computed using the code provided by the authors, where we set the penalty
parameters lambda.alpha and lambda.beta to be $\log(p)/n$ and $\log(q)/n$,
respectively. Witten et al. (2009)’s method is implemented using the R package
PMA with $l_{1}$ penalty, using the default tuning parameters. Finally, we
consider $N=1000$ Monte Carlo replications as before.

Table 4 and Table 5 tabulate the absolute bias and the standard deviation of
$|\widehat{x}_{i}|$ for $i=1$ and $20$, respectively, estimated using the 1000
Monte Carlo samples. Recall from Section 5 that $x^{0}_{1}$ is nonzero but
$x^{0}_{20}$ is zero.

Bias in the estimation of $|x^{0}_{1}|$: Table 4 entails that the de-biased
estimators of $x^{0}_{1}$ almost always outperform the remaining estimators in
terms of the absolute bias, and the difference is more prominent when the
signal strength is small. The only exception is the high signal strength
setting, where sometimes the bias of the initial COLAR estimator is so small
such that the de-biasing step does not lead to further improvement. The bias
of the de-biased estimator and COLAR, in general, is close, and they exhibit
the same pattern. A sharp decrease in the bias of the COLAR and the de-biased
estimator can be observed at signal strength $0.5$ and $0.9$, respectively,
for identity and sparse inverse matrix. The QQ plots in Figure 7 and the
histograms in Figure 5 also reveal that the de-biased estimators attain
asymptotic normality at these signal strength. These observations explain why
the ordinary confidence intervals in Section 5, which rely on Corollary 4.2,
have poor coverage at lower signal strength in the above cases. In the sparse
inverse case, the bias of Witten et al. (2009)’s estimator stays substantially
high, and increases with the signal strengths for high $p$, $q$. This is
unsurprising because Witten et al. (2009)’s method is best suited for diagonal
covariance matrices.

Bias in the estimation of $|x^{0}_{20}|$: In this case, the SCCA estimators
have much smaller bias than our de-biased estimator, which is expected because
sparse estimators would generally set this co-ordinate to zero. However, as
the QQ plots in Figure 8 indicate, the de-biased estimator attains asymptotic
normality pretty quickly, even at low signal strength, while the initial COLAR
estimator stays quite non-normal unless the signal strength is high. This
observation explains the satisfactory performance of the confidence intervals
for $x^{0}_{20}$.

Method | Identity | Sparse Inverse
---|---|---
| $n=500,\ p=80,\ q=80$
| $\rho_{0}=0\cdot 2$ | $\rho_{0}=0\cdot 5$ | $\rho_{0}=0\cdot 9$ | $\rho_{0}=0\cdot 2$ | $\rho_{0}=0\cdot 5$ | $\rho_{0}=0\cdot 9$
PMA | 28 (18) | 9.1 (15) | 2.9 (3.2) | 29 (10) | 43 (22) | 31 (23)
Mai & Zhang | 28 (15) | 5.0 (6.3) | 2.5 (3.1) | 30 (10) | 37 (37) | 54 (6.4)
COLAR | 28 (19) | 5.0 (6.4) | 2.2 (2.8) | 29 (25) | 7.3 (9.6) | 3.0 (2.6)
Db | 21 (23) | 4.5 (5.6) | 2.1 (2.6) | 25 (24) | 6.0 (8.7) | 2.3 (2.5)
| $n=500,\ p=300,\ q=200$
PMA | 29 (16) | 33 (29) | 8.5 (16) | 29 (13) | 47 (21) | 57 (36)
Mai & Zhang | 31 (9.3) | 5.2 (6.4) | 2.5 (3.1) | 30 (5.6) | 39 (40) | 54 (7.4)
COLAR | 30 (15) | 5.3 (6.6) | 2.2 (2.7) | 29 (27) | 42 (38) | 3.2 (2.5)
Db | 24 (25) | 4.7 (5.8) | 2.1 (2.6) | 26 (26) | 39 (40) | 2.3 (2.4)
| $n=500,\ p=600,\ q=200$
PMA | 30 (15) | 40 (30) | 17 (27) | 42 (20) | 67 (34) | 87 (49)
Mai & Zhang | 31 (8.3) | 4.7 (6.1) | 2.4 (3.0) | 43 (9.7) | 44 (34) | 3.8 (4.6)
COLAR | 31 (11) | 5.0 (6.3) | 2.1 (2.7) | 43 (20) | 42 (49) | 3.9 (4.2)
Db | 25 (26) | 4.3 (5.4) | 2.0 (2.5) | 37 (37) | 35 (40) | 3.4 (3.7)
Table 4: Table of the estimated bias of $|\widehat{x}_{1}|$. The standard deviation estimate is given in the parentheses. The bias and the standard error is estimated from 1000 Monte Carlo samples. All entries are scaled by $10^{-2}$. Here PMA: Penalized Multivariate Analysis Witten et al. (2009); Db: The de-biased estimator. Method | Identity | Sparse Inverse
---|---|---
| $n=500,\ p=80,\ q=80$
| $\rho_{0}=0\cdot 2$ | $\rho_{0}=0\cdot 5$ | $\rho_{0}=0\cdot 9$ | $\rho_{0}=0\cdot 2$ | $\rho_{0}=0\cdot 5$ | $\rho_{0}=0\cdot 9$
PMA | 2.5 (8.3) | 1.6 (4.0) | 1.8 (3.2) | 2.9 (9.3) | 2.5 (8.5) | 1.1 (5.0)
Mai & Zhang | 1.3 (6.7) | 0.05 (0.5) | 0 (0) | 0.98 (5.2) | 0.32 (2.7) | 0 (0)
COLAR | 1.0 (6.2) | 0.02 (0.34) | 0 (0) | 0.59 (4.2) | 0.06 (1.4) | 0 (0)
Db | 7.6(10) | 4.4 (5.5) | 1.7 (2.1) | 6.9 (9.2) | 4.6 (5.9) | 1.7 (2.1)
| $n=500,\ p=300,\ q=200$
PMA | 1.7 (6.1) | 1.6 (5.2) | 1.5 (3.3) | 1.9 (6.5) | 1.9 (6.5) | 1.9 (6.2)
Mai & Zhang | 0.40 (3.6) | 0.01 (0.24) | 0 (0) | 0.31 (2.6) | 0.07 (0.91) | 0 (0)
COLAR | 0.31 (3.5) | 0 (0.07) | 0 (0) | 0.16 (2.2) | 0.07 (1.5) | 0 (0)
Db | 6.8 (8.9) | 4.2 (5.3) | 1.6 (2.0) | 6.3 (8.0) | 5.0 (6.5) | 1.6 (2.1)
| $n=500,\ p=600,\ q=200$
PMA | 1.3 (4.1) | 1.3 (3.9) | 1.2 (3.2) | 1.8 (5.9) | 1.8 (6.0) | 1.8 (5.9)
Mai & Zhang | 0.28 (2.8) | 0.02 (0.3) | 0 (0) | 0.15 (1.2) | 0.12 (1.1) | 0 (0)
COLAR | 0.17 (2.3) | 0 (0.05) | 0 (0) | 0.17 (1.9) | 0.09 (1.4) | 0 (0)
Db | 6.3 (8.0) | 4.2 (5.3) | 1.6 (2.0) | 6.1 (7.7) | 5.1 (6.6) | 1.7 (2.1)
Table 5: Table of the estimated bias of $|\widehat{x}_{20}|$. The standard
deviation estimate is given in the parentheses. The bias and the standard
error is estimated from 1000 Monte Carlo samples. All entries are scaled by
$10^{-2}$. PMA: Penalized Multivariate Analysis Witten et al. (2009); Db: The
de-biased estimator.

###  8.2 Extra plots: simulation

(a) Ordinary confidence intervals for identity matrix

(b) Ordinary confidence intervals for sparse inverse matrix

Figure 3: Ordinary confidence intervals for $(x^{0}_{20})^{2}$

(a) Conservative confidence intervals for identity matrix

(b) Conservative confidence intervals for sparse inverse matrix

Figure 4: Conservative confidence intervals for $(x^{0}_{20})^{2}$

(a) Identity matrix

(b) Sparse inverse matrix

Figure 5: Histograms of $(\widehat{x}^{0}_{1})^{2}$: the estimates were
centered by $(x^{0}_{1})^{2}$ and scaled by
$4|x^{0}_{1}|\sigma_{i}n^{-1/2}$, where $\sigma_{i}$ is as in Theorem 4.1.
Preliminary estimates in blue and de-biased versions in red. A standard normal
curve is imposed.

(a) Ordinary intervals

(b) Conservative intervals

Figure 6: Effect of doubling $(p,q,n)$: note that the coverage of both
ordinary and conservative confidence intervals increase. Here the underlying
covariance matrices are taken to be identity. The coverage and the rejection
probability of the tests are calculated using 1000 Monte Carlo samples.

(a) QQ-plots for identity matrix

(b) QQ-plots for sparse inverse matrix

Figure 7: QQ plots for $\widehat{x}_{1}^{2}$.

(a) QQ-plots for identity matrix

(b) QQ-plots for sparse inverse matrix

Figure 8: QQ plots for $\widehat{x}_{20}^{2}$.

###  8.3 Extra plot: data application

(a) Confidence intervals for gene measurements: first half of the genes

(b) Confidence intervals for gene measurements: second half of the genes

(c) Confidence intervals for protein measurements

Figure 9: Confidence intervals for pathway (a) Cytokine-Cytokine receptor
interaction pathway.

(a) Confidence intervals for gene measurements

(b) Confidence intervals for protein measurements

Figure 10: Confidence intervals for pathway (b), i.e. Adipocytokine signal
pathway.

(a) Genes in pathway (a)

(b) Genes in pathway (b)

(c) Proteins

Figure 11: Variance plot (in absolute values) of genes and proteins: here
darker color means higher correlation. The color scale is given to the right.

(a) Pathway (a): first 100 genes

(b) Pathway (b)

(c) Pathway (a): last 131 genes

Figure 12: Covariance plot: plotting the absolute value of covariance between
genes and proteins. Here proteins are in the X axis and genes are in the Y
axis. Darker color is associated with higher correlation. The color scale is
given to the right. Because pathway (a) has 231 genes in comparison to only 94
proteins, it is split into two parts for easier representation.

##  9  Modified COLAR Estimators

###  9.1 Additional notation

We will need a few additional notation for presenting our method for
estimating $\alpha_{0}$ and $\beta_{0}$. The trace inner product between two
matrices $A,B\in\mathbb{R}^{p\times q}$ is defined by $\langle
A,B\rangle=\tr(A^{T}B)$. For any matrix $A$, $A_{-j}$ will denote the matrix
obtained by deleting the $j$th column of $A$. Further, we let $A_{-j,j}$
denote the vector derived from $A_{j}$ by deleting its $j$th element. We
denote by $A_{-j,-j}$ the matrix obtained by deleting the $j$th column and
$j$th row of $A$. We let $A_{j,-j}$ be the the vector obtained by deleting the
$j$th element of the $j$th row of $A$. Also, for any symmetric matrix $A$, we
let $A^{1/2}$ denote the matrix $PD^{1/2}P^{T}$ where $PDP^{T}$ is a spectral
decomposition of $A$. Moreover, for $j\in\mathbb{N}$, we will use
$e_{j}\in\mathbb{R}^{n}$ to represent a unit vector with a one in $j$th
position and $0$’s elsewhere with $n$ determined by context. Finally, we
define the class of orthogonal matrices

$\mathcal{O}(p,r)=\bigg{\\{}U\in\mathbb{R}^{p\times
r}:U^{T}U=I_{r}\bigg{\\}}.$ (17)

###  9.2 Modified COLAR (Gao et al., 2017) algorithm

Our algorithm runs in three stages and is motivated by Gao et al. (2017). The
main difference between Gao et al. (2017)’s algorithm and ours is that we set
the parameter $r$ in the COLAR algorithm from Gao et al. (2017) to a working
and possibly misspecified value of $1$. This step results in some differences
in the analytic details as well as some simplifications of the original
algorithm from Gao et al. (2017). To set up the algorithm, we split the data
into two equal parts indexed by $\\{0,1\\}$, and calculate the empirical
estimators ($\widehat{\Sigma}_{n,xy}^{(0)},\widehat{\Sigma}_{n,x}^{(0)}$) and
($\widehat{\Sigma}_{n,xy}^{(1)},\widehat{\Sigma}_{n,x}^{(1)}$) based on the
sub-samples $(X^{(0)},Y^{(0)})$ and $(X^{(1)},Y^{(1)})$, respectively. Here
the superscripts refers to the specific sub-sample used to calculate the
empirical estimators.

The first stage of our method produces a good preliminary estimator of
$\alpha_{0}$ and $\beta_{0}$, which is used by the second stage to produce
improved estimators. This step is similar to Gao et al. (2017) but we present
the details for the sake of completeness as well as the ease of using the
notation from this stage of the algorithm while developing the analytic
justification of the method. The first stage solves a convex relaxation of (1)
to obtain a set of preliminary estimators $\widehat{\alpha}_{n}^{(0)}$ and
$\widehat{\beta}_{n}^{(0)}$. The convex relaxation hinges on the idea that (1)
can be written as a convex program after the change of variable
$F=\alpha\beta^{T}$. To see this, first note that since
$\alpha^{T}\widehat{\Sigma}_{n,xy}^{(0)}\beta=tr(\widehat{\Sigma}_{n,yx}^{(0)}F)$,
the objective function in (1) can be written as a linear functional of $F$.
The feasible set of (1) can also be written in terms of $F$ but it is not a
convex set in general. However, since the objective is linear in $F$, if we
replace the feasible set with its convex hull $\mathcal{G}$, the solutions
remain unchanged. The latter follows since a linear function is always
maximized at the boundary of any convex set. From Gao et al. (2017) it then
follows that $\mathcal{G}$ takes the form

$\displaystyle\mathcal{G}=\bigg{\\{}F\in\mathbb{R}^{p\times q}\ :\
\|\widehat{\Sigma}_{n,x}^{0,1/2}F\widehat{\Sigma}_{n,y}^{0,1/2}\|_{*}\leq
1,\quad\|\widehat{\Sigma}_{n,x}^{0,1/2}F\widehat{\Sigma}_{n,y}^{0,1/2}\|_{op}\leq
1\bigg{\\}},$ (18)

and our optimization problem reduces to maximizing $tr({\Sigma}_{yx}F)$ with
respect to $F\in\mathcal{G}$. To obtain sparse solutions, we also add an
$l_{1}$ penalty to the objective function. Therefore, at the end, the first
stage solves

$\displaystyle\mathmakebox[width("$\underset{\displaystyle
F\in\mathcal{G}}{\mathrm{subject~{}to}}$")][l]{\underset{\displaystyle
F\in\mathcal{G}}{\mathrm{maximize}}}\quad
tr(\widehat{\Sigma}_{n,yx}^{(0)}F)-\lambda_{1}\|F\|_{1},\hfil\hfil\hfil\hfil$
(19)

where $\lambda_{1}$ is a tuning parameter and
$\|F\|_{1}=\sum_{i=1}^{q}\sum_{j=1}^{p}|F_{ij}|$ is the vector $l_{1}$ norm of
the matrix $F$. We take $\lambda_{1}=C\lambda$, where $\lambda$ is as defined
in (4), and $C>0$ is some constant, whose value will be chosen later. The
above optimization program gives an estimate $\widehat{F}_{n}$ of
$F_{0}=\alpha_{0}\beta_{0}^{T}$. The first pair of left and right singular
vectors of $\widehat{F}_{n}$ give the preliminary estimators of $\alpha$ and
$\beta$, which we denote by $\widehat{\alpha}_{n}^{(0)}$ and
$\widehat{\beta}_{n}^{(0)}$, respectively.

Our second stage is where we differ from Gao et al. (2017). This modified
second stage improves upon the preliminary estimators and estimates
$\alpha_{0}$ and $\beta_{0}$ up to a sign flip. This is crucial, since
estimating all the canonical directions simultaneously, as developed in Gao et
al. (2017), does not lend itself to identifying the first directions only up
to a sign flip. Henceforth, we will only consider the estimation of
$\alpha_{0}$ because the estimation of $\beta_{0}$ will be similar. To obtain
an improved estimator of $\alpha_{0}$, the second stage solves

$\displaystyle\mathmakebox[width("$\underset{\displaystyle
x\in\mathbb{R}^{p}}{\mathrm{subject~{}to}}$")][l]{\underset{\displaystyle
x\in\mathbb{R}^{p}}{\mathrm{minimize}}}\quad
x^{T}\widehat{\Sigma}_{n,x}^{(1)}x-2x^{T}\widehat{\Sigma}_{n,xy}^{(1)}\widehat{\beta}_{n}^{0}+\lambda_{2}\|x\|_{1}\hfil\hfil\hfil\hfil$
(20)

where $\lambda_{2}$ is a penalizing parameter. We will take
$\lambda_{2}=C\lambda$ for $\lambda$ defined in (4) and some constant $C>0$
whose value will be chosen later. When the observations are centered, i.e.
$\widehat{\Sigma}_{n,x}^{(1)}=(X^{(1)})^{T}X^{(1)}$ and
$\widehat{\Sigma}_{n,xy}^{(1)}=(X^{(1)})^{T}Y^{(1)}$, (20) can be re-written
as

We will denote the solution to (20) by $\tilde{x}_{n}$. Gao et al. (2017) uses
a group lasso penalty instead of the $l_{1}$ penalty in (20). These panalties
are, however, equivalent when $x$ is a vector, as in our case. Had we been
estimating more than one leading canonical vector, as in Gao et al. (2017)’s
case, $x$ would be a matrix, and the group lasso penalty is no longer
equivalent to the $l_{1}$ penalty.

The third stage is the normalization step, which simply sets

$\widehat{\alpha}_{n}=\begin{cases}\tilde{x}_{n}\Big{(}(\tilde{x}_{n})^{T}\widehat{\Sigma}_{n,x}\tilde{x}_{n}\Big{)}^{-1/2}&(\tilde{x}_{n})^{T}\widehat{\Sigma}_{n,x}\tilde{x}_{n}>0\\\
0&o.w.\end{cases}$ (21)

It will be later shown in Lemma 19.3 that the quadratic form
$(\tilde{x}_{n})^{T}\widehat{\Sigma}_{n,x}\tilde{x}_{n}$ is non-zero and
$\widehat{\alpha}_{n}=\tilde{x}_{n}\Big{(}(\tilde{x}_{n})^{T}\widehat{\Sigma}_{n,x}\tilde{x}_{n}\Big{)}^{-1/2}$
with probability tending to one. Our third stage is slightly different from
Gao et al. (2017), who used the sample covariance matrix from a third part of
the data to normalize $\tilde{x}_{n}$, where we use the full covariance matrix
$\widehat{\Sigma}_{n,x}$. Since we want to estimate only $\alpha_{0}$ instead
of the whole matrix $U$ as in Gao et al. (2017), normalization is simpler in
our case, which circumvents the need of the stage final data splitting. We
remark on passing that we could use Gao et al. (2017)’s third step as well,
and the asymptotics would remain the same. However, we avoid three way data
splitting because unnecessary data splitting may not be beneficial in finite
sample. For convenience, we list the modified COLAR algorithm in Algorithm 1.
Somewhere write the full form of COLAR in main text; as well as in supplement.

Input:

$\widehat{\Sigma}_{n,x}^{(i)}$, $\widehat{\Sigma}_{n,xy}^{(i)}$,
$\widehat{\Sigma}_{n,y}^{(i)}$ $(i=0,1)$, $\lambda_{1}$ and $\lambda_{2}$.

Stage 1:

1\. Solve the convex program
$\widehat{F}_{n}=\operatorname*{arg\,max}_{F\in\mathcal{G}}\bigg{\\{}tr(\widehat{\Sigma}_{n,yx}^{(0)}F)-\lambda_{1}\|F\|_{1}\bigg{\\}}$
where $\mathcal{G}$ is as in (18).

2\. Obtain the first pair of singular vectors
$\widehat{\alpha}_{n}^{(0)}\in\mathbb{R}^{p}$ and
$\widehat{\beta}_{n}^{(0)}\in\mathbb{R}^{q}$ of $\widehat{F}_{n}$.

Stage 2:

Solve the convex program
$\tilde{x}_{n}=\operatorname*{arg\,min}_{x\in\mathbb{R}^{p}}\\{x^{T}\widehat{\Sigma}_{n,x}^{(1)}x-2x^{T}\widehat{\Sigma}_{n,xy}^{(1)}\widehat{\beta}_{n}^{0}+\lambda_{2}\|x\|_{1}\\}$

Stage 3: Set

$\widehat{\alpha}_{n}=\tilde{x}_{n}\Big{(}(\tilde{x}_{n})^{T}\widehat{\Sigma}_{n,x}\tilde{x}_{n}\Big{)}^{-1/2}$

Output: $\widehat{\alpha}_{n}$

Algorithm 1 Modified COLAR (Gao et al., 2017) algorithm

###  9.3 Asymptotic properties of the COLAR estimator

In Section 4, we noted that $n^{1/2}$-consistency of the de-biased estimators
requires some restrictions on the $l_{1}$ and $l_{2}$ errors of the
preliminary estimators of $\alpha_{0}$ and $\beta_{0}$, which are satisfied by
our $\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$.

######  Theorem 9.1.

Suppose Assumption 2.1 and Assumption 2.1 hold. Further suppose
$s=s_{U}+s_{V}$ satisfies $s\lambda\to 0$, where $\lambda$ is as in (4), and
$s=o(p)$. Then there exist $C_{1}$, $C_{2}>0$ such that for
$\lambda_{1}=C\lambda$ with $C>C_{1}$, and $\lambda_{2}=C^{\prime}\lambda$
with $C^{\prime}>C_{2}$, the estimators $\widehat{\alpha}_{n}$ and
$\widehat{\beta}_{n}$ defined in (21) satisfy Condition 4.1 with

$\displaystyle\kappa=\begin{cases}1/2&\text{if }r=1\\\
1&\text{o.w.}\end{cases}$ (22)

Note that the sparsity requirement on $s$ is $s\lambda=o(1)$, which is a
weaker condition than our Assumption 4.1 that requires
$s^{2\kappa}\lambda^{2}=o(n^{-1/2})$. Fact 1 indicates that Assumption 4.1
implies $s\lambda=o(1)$.

######  Remark 9.2.

When $r=1$, the proof of Theorem 9.1 implies that a slightly stronger result
holds than that implied by Condition 4.1. More explicitly, the $l_{1}$ and
$l_{2}$ errors of $\widehat{\alpha}_{n}$ depends only on $s_{U}$, and not on
$s_{V}$. Similarly, the asymptotics of $\widehat{\beta}_{n}$ depend only on
$s_{V}$. To be more precise,

$\inf_{w\in\\{\pm
1\\}}\|w\widehat{\alpha}_{n}-\alpha_{0}\|_{1}=O_{p}(s_{U}\lambda),\quad\inf_{w\in\\{\pm
1\\}}\|w\widehat{\beta}_{n}-\beta_{0}\|_{1}=O_{p}(s_{V}\lambda)$

and

$\inf_{w\in\\{\pm
1\\}}\|w\widehat{\alpha}_{n}-\alpha_{0}\|_{2}=O_{p}(s_{U}^{1/2}\lambda),\quad\inf_{w\in\\{\pm
1\\}}\|w\widehat{\beta}_{n}-\beta_{0}\|_{2}=O_{p}(s_{V}^{1/2}\lambda).$

The above result is substantially sharper than that implied by the statement
of Condition 4.1 if $s_{U}\ll s_{V}$ or vice versa.

The optimal value of $\lambda_{1}$ and $\lambda_{2}$ rely on $C_{1}$ and
$C_{2}$, which depend unknown quantities like $M$ in Assumption 2.1.
Therefore, cross-validation may be required to choose these tuning parameters
efficiently. According to Gao et al. (2017), there is scope of improving the
algorithm so that it adapts to the unknown $M$. However, it is beyond the
scope of the current paper.

######  Remark 9.3 (Chen et al. (2013)’s estimators).

Although Chen et al. (2013) uses an iterative thresholding type method to
estimate $\alpha_{0}$ and $\beta_{0}$ upto a sign flip and the resulting
estimators attain the minimax rate in $l_{2}$ norm, they consider the rank one
model. It is remains unknown whether their method continues to work similarly
for $r>1$ case while estimating purely the leading canonical directions up to
a sign flip. Here we discuss one potential roadblock on its straightforward
extension to the general $r>1$ case. The theoretical guarantees of Chen et al.
(2013)’s method rely heavily on the initial estimators. To obtain these
initial estimators, they apply singular value decomposition on a suitably
chosen estimator of ${\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}$. When
$r=1$, the matrix ${\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}$ has rank
one, and its leading singular vectors are proportional to $\alpha_{0}$ and
$\beta_{0}$. Therefore the above-mentioned initialization method works.
However, when $r>1$, unless ${\Sigma}_{x}$ and ${\Sigma}_{y}$ are identity,
the leading singular vectors of
${\Sigma}_{x}^{-1}{\Sigma}_{xy}{\Sigma}_{y}^{-1}$ are no longer proportional
to $\alpha_{0}$ and $\beta_{0}$. Therefore, the idea behind the initialization
method of Chen et al. (2013) ceases to work for $r>1$.

##  10  Nodewise Lasso Estimator

###  10.1 The main algorithm

The nodewise lasso algorithm was first implemented by Meinshausen et al.
(2006), who used the name graphical lasso. Meinshausen et al. (2006) showed
that the $p\times p$ dimensional precision matrix can be estimated by
regressing each of the $p$ variables against the other; see also the nodewise
regression of van de Geer et al. (2014). Although originally invented for
precision matrix estimation, the basic idea of nodewise lasso applies to the
inversion of any real symmetric matrix; cf. Janková & van de Geer (2018).
Asymptotic garuantees, however, can be established only if the input matrix
consistently estimates a positive definite matrix. Also, for the most part,
the asymptotics of the nodewise lasso algorithm is case-specific, which is to
say that the convergence results solely depend on the matrix to be inverted,
which is $\widehat{H}_{n}(\widehat{x}_{n},\widehat{y}_{n})$ in our case. For
the sake of completeness, we include this algorithm in our paper; see
Algorithm 2.

Input:

$A\in\mathbb{R}^{m\times m}$ where $m\in\mathbb{N}$, positive penalty
parameters $(\lambda^{nl}_{j},B_{j})$, $j=1,\ldots,m$.

for $j=1,\dots m:$

  * NL1.

Compute any stationary point $\widehat{\eta}_{j}$ of the minimization program

$\displaystyle\mathmakebox[width("$\underset{\displaystyle\eta_{j}\in\mathbb{R}^{p+q-1},\|\eta_{j}\|_{1}\leq
B_{j}}{\mathrm{subject~{}to}}$")][l]{\underset{\displaystyle\eta_{j}\in\mathbb{R}^{p+q-1},\|\eta_{j}\|_{1}\leq
B_{j}}{\mathrm{minimize}}}\quad\eta_{j}^{T}A_{-j,-j}\eta_{j}-2A_{-j,j}^{T}\eta_{j}+\lambda^{nl}_{j}\|\eta_{j}\|_{1},\hfil\hfil\hfil\hfil$
(23)

where we remind the readers that $A_{-j,-j}$ is the matrix obtained by
deleting the $j$th row and the $j$th column of the matrix $A$, and $A_{-j,j}$
is the vector obtained by deleting the $j$th element of $A_{j}$.

  * NL2.

Compute the estimator of the noise-level

$\widehat{\tau}_{j}^{2}=\widehat{\Gamma}_{j}^{T}A\widehat{\Gamma}_{j}+\dfrac{1}{2}\lambda^{nl}_{j}\|\widehat{\eta}_{j}\|_{1},$
(24)

where

$\widehat{\Gamma}_{j}=(-(\widehat{\eta}_{j})_{1},\ldots,-(\widehat{\eta}_{j})_{j-1},1,-(\widehat{\eta}_{j})_{j+1},\ldots,-(\widehat{\eta}_{j})_{m})$
(25)

Set
$\overline{A}=[\widehat{\Gamma}_{1}/\widehat{\tau}_{1}^{2},\ldots,\widehat{\Gamma}_{m}/\widehat{\tau}_{m}^{2}]$

Output: $\overline{A}$

Algorithm 2 Non-convex Nodewise Lasso

A couple of remarks are in order. First, the $l_{1}$ penalty in (23) is
introduced to enforce a sparse solution. Second, the constraint
$\|\eta_{j}\|_{1}\leq B_{j}$ ensures a bounded solution to the problem.
Without this boundedness condition, the optimization problem (23) can become
unbounded since $A_{-j,-j}$ is potentially singular. Third, there is no
guarantee that $\overline{A}$ will be symmetric when $A$ is symmetric.
Therefore, we have to compute the full matrix $\overline{A}$ even if $A$ is
known to be symmetric. Finally, notice that Algorithm 2 does not require us to
solve (23), which is possibly non-convex, since a stationary point of (23)
suffices. In Section 10.4, we will demonstrate how to choose the tuning
parameters $\lambda^{nl}_{j}$ and $B_{j}$.

######  Remark 10.1 (Possible other choices of $\widehat{\Phi}_{n}$).

The de-biasing literature borrows nodewise lasso from precision matrix
estimation literature. Other methods for precision matrix estimation, e.g.
Constrained $l_{1}$-minimization for Inverse Matrix Estimation aka CLIME (Cai
et al., 2011), the graphical lasso aka GLASSO (Friedman et al., 2008) etc. may
also be used in place of nodewise lasso to construct $\widehat{\Phi}_{n}$
provided Condition 4.1 is satisfied under realistic structural assumptions. In
this regard, CLIME solves convex optimization problems and has fast
implementation. It has also seen application in context of de-biasing (Neykov
et al., 2018). We conjecture that if the columns of $\Phi^{0}$ are bounded in
$l_{1}$ norm, then the CLIME estimator satisfies the desired Condition 4.1 as
well. This requirement, however, is stricter than that of the nodewise lasso;
see Assumption 10.4 in Section 10.4. To keep our discussions focused, we
refrain from further discussion on the asymptotics of CLIME here. Similar to
CLIME, GLASSO also has fast implementation and is widely used in precision
matrix estimation. However, the current literature lacks results supporting
its consistency. There is, instead, some evidence against its asymptotic
convergence to the precision matrix, at least in $l_{\infty}$ norm (Mazumder &
Hastie, 2012).

###  10.2 Intuition behind the nodewise lasso algorithm

To provide an intuition behind why Algorithm 2 works, we argue that if the
input matrix $A$ in Algorithm 2 is positive definite, the algorithm outputs
$A^{-1}$ when the penalty parameters $\lambda^{nl}_{j}$’s are set to zero. To
that end, we first invoke a standard linear algebra result (Rao &
Bhimasankaram, 2000, cf.).

######  Lemma 10.2.

Suppose $m\in\mathbb{N}$ and $A$ is an $m\times m$ positive definite matrix.
Then $A_{-j,-j}$ is invertible for $j=1,\ldots,m$. Moreover,

$(A^{-1})_{j,j}=\dfrac{1}{A_{j,j}-A_{j,-j}^{T}A_{-j,-j}A_{-j,j}}$
$(A^{-1})_{-j,j}=-(A^{-1})_{j,j}(A_{-j,-j})^{-1}A_{-j,j}.$

Defining $\eta_{j}=(A_{-j,-j})^{-1}A_{-j,j}$, we note that

$\operatorname*{arg\,min}_{\eta\in\mathbb{R}^{p-1}}(\eta^{T}A_{-j,-j}\eta-2A_{-j,j}^{T}\eta)=\eta_{j},\quad
j=1,\ldots,m.$

Also, in parallel with (2), we define

$\Gamma_{j}=(-(\eta_{j})_{1},\ldots,-(\eta_{j})_{j-1},1,-(\eta_{j})_{j+1},\ldots,-(\eta_{j})_{r-1}),$

. Then it follows that

$\tau_{j}^{2}=\Gamma_{j}^{T}A\Gamma_{j}=A_{j,j}-A_{j,-j}^{T}A_{-j,-j}A_{-j,j}\stackrel{{\scriptstyle(a)}}{{=}}1/(A^{-1})_{j,j},$
(26)

where (a) follows from Lemma 10.2. Applying Lemma 10.2 again, we can show that
the $j$th element of the output matrix $\overline{A}$ equals

$\Gamma_{j}/\tau_{j}^{2}=(A^{-1})_{j}.$ (27)

###  10.3 Nodewise lasso for our case

In this section, we discuss the finite sample properties of our nodewise lasso
estimator $\widehat{\Phi}_{n}$. We begin with some implications of the
discussion in Section 10.2 for the special case when the input matrix
$A=H^{0}$. First, note that, in this case, for $j=1,\ldots,p+q$,

$\eta_{j}^{0}=\operatorname*{arg\,min}_{\eta\in\mathbb{R}^{p-1}}\Big{(}\eta^{T}H^{0}_{-j,-j}\eta-2(H^{0}_{-j,j})^{T}\eta\Big{)}$

satisfies

$\eta_{j}^{0}=(H^{0}_{-j,-j})^{-1}H^{0}_{-j,j}.$ (28)

Moreover, (26) implies for $j=1,\ldots,p+q$,

$(\tau^{0}_{j})^{2}=H^{0}_{j,j}-(H^{0}_{j,-j})^{T}H^{0}_{-j,-j}H^{0}_{-j,j}$
(29)

satisfies $(\tau^{0}_{j})^{2}=(\Phi^{0}_{j,j})^{-1}$. From (27), it then
follows that

$(\Phi^{0})_{-j,j}=-\eta_{j}^{0}/(\tau^{0}_{j})^{2}.$ (30)

Because $\widehat{\eta}_{j}$ is a stationary point of (23), it satisfies the
KKT condition, which takes the form

$-2A_{j,-j}+2A_{-j,-j}\widehat{\eta}_{j}+\lambda^{nl}_{j}\partial\|\widehat{\eta}_{j}\|_{1}=0,$

where $\lambda^{nl}_{j}$ is as in (23) and
$\partial\|\widehat{\eta}_{j}\|_{1}$ is the partial derivative of the $l_{1}$
norm evaluated at $\widehat{\eta}_{j}$. It then follows that (cf. Section 3.1
of Janková & van de Geer, 2018)

$A_{j}^{T}\widehat{\Gamma}_{j}=\widehat{\tau}_{j}^{2}\quad\text{and}\quad\|A_{-j}^{T}\widehat{\Gamma}_{j}\|_{\infty}\leq\lambda_{j}^{nl}/2$
(31) $A^{T}\overline{A}-I_{p+q}|_{\infty}=O\Big{(}\max_{1\leq j\leq
p+q}\lambda^{nl}_{j}/{\widehat{\tau}}_{j}^{2}\Big{)}.$

###  10.4 Asymptotic properties of the nodewise lasso estimator

In this Section, we will show that the nodewise lasso estimator satisfies
Condition 4.1 under some regulatory conditions. We will go through these
regulatory conditions first.

Recall from (28) the definition of $\eta_{j}^{0}$. We will require the number
of non-zero elements in $\eta_{j}^{0}$, i.e. $\|\eta_{j}^{0}\|_{0}$, to be
small, which is in parallel with Janková & van de Geer (2018). [ Assumption on
the column sparsity of $\Phi^{0}$] $\max_{1\leq j\leq
p+q}\|{\eta^{0}}_{j}\|_{0}=O(s)$, where $s=s_{U}+s_{V}$. Since
$\absolutevalue{\|\eta^{0}_{j}\|_{0}-\|\Phi^{0}_{j}\|_{0}}\leq 1$ by (30), a
restriction on $\|\eta^{0}_{j}\|_{0}$ actually induces a restriction on
$\|\Phi^{0}_{j}\|_{0}$, which explains the nomenclature of Assumption 10.4.

Assumption 10.4 can be hard to decipher, and it may be hard to verify.
Therefore we will now give a sufficient condition for Assumption 10.4. Lemma
21.1 in Supplement 21.1 gives the explicit form of $\Phi^{0}$, which indicates
that

$\|{\eta^{0}}_{j}\|_{0}\leq
s_{U}+s_{V}+\|({\Sigma}_{x})_{j}^{-1}\|_{0}+\|({\Sigma}_{y})_{j}^{-1}\|_{0}\quad(j=1,\ldots,p).$

Therefore, we only require the column sparsities of ${\Sigma}_{x}^{-1}$ and
${\Sigma}_{y}^{-1}$ to be $O(s)$ for Assumption 10.4 is satisfied. This
sparsity requirement is formulated as Condition 10.4.

[ A sufficient condition for Assumption 10.4] The maximum number of non-zero
elements per column of ${\Sigma}_{x}^{-1}$ or ${\Sigma}_{y}^{-1}$ is $O(s)$.
Sparsity restriction on the columns of ${\Sigma}_{x}^{-1}$ and
${\Sigma}_{y}^{-1}$ is more intuitive than sparsity restriction on $\Phi^{0}$.
It is a well known fact that ${\Sigma}_{x}^{-1}$ or ${\Sigma}_{y}^{-1}$ is
sparse if the partial correlation between the $X_{i}$’s or the $Y_{i}$’s are
mostly zero, which may be satisfied when only a few of these variables
interact among themselves. The latter is not unusual in high dimensional
genomic data because genes, proteins etc. generally form clusters. Such
sparsity restrictions are also common in the literature; cf. Bühlmann & Van De
Geer (2011); Janková & van de Geer (2018).

Now we are ready to state the main theorem of this section.

######  Theorem 10.3 (Nodewise lasso Theorem).

Suppose Assumptions 2.1, 2.1, 4.1, 10.4 hold, and the preliminary estimators
$\widehat{x}_{n}$ and $\widehat{y}_{n}$ satisfy Condition 4.1. Further suppose

$\|\eta_{j}^{0}\|_{1}\leq B_{j}\leq C_{T}s^{1/2}\quad(j=1,\ldots,p+q)$

for some $C_{T}>0$. Then there exists an absolute constant $C>0$ depending on
$C_{T}$ such that for

$\lambda^{nl}_{j}=C\lambda,$

the estimator $\widehat{\Phi}_{n}$ obtained by feeding
$\widehat{H}_{n}(\widehat{x}_{n},\widehat{y}_{n})$ to Algorithm 2 satisfies
Condition 4.1. here $\lambda$ is as in (4).

The nodewise lasso estimator $\widehat{\Phi}_{n}$ depends on $\widehat{x}_{n}$
and $\widehat{y}_{n}$ via $\widehat{H}_{n}(\widehat{x}_{n},\widehat{y}_{n})$,
which does not rely on the sign of $\widehat{x}_{n}$ and $\widehat{y}_{n}$.
Hence, the asymptotics of $\widehat{\Phi}_{n}$, unlike the de-biased
estimators, is unaffected by the sign flip of $\widehat{\alpha}_{n}$ and
$\widehat{\beta}_{n}$.

##  11  Connection to Related Literature

The study of asymptotic inference in the context of SCCA naturally connects to
the popular research direction of de-biased/de-sparsified inference in high
dimensional models (Zhang & Zhang, 2014; Javanmard & Montanari, 2014; van de
Geer et al., 2014; Janková & van de Geer, 2018; Zhang & Zhang, 2014; Ning et
al., 2017; Neykov et al., 2018; Janková & van de Geer, 2017; Janková & Van De
Geer, 2016; Cai et al., 2017; Mitra et al., 2016; Bellec & Zhang, 2019). This
line of research, starting essentially from the seminal work of Zhang & Zhang
(2014), more or less follows the general prescription laid out in Section 3.1.
Similar to our case, these methods also often depend on potentially high
dimensional parameters – and thereby require initial good estimators of them.
For example, asymptotically valid confidence interval for the coordinates of a
sparse linear regression vector relies critically on good initial estimators
of the regression vector and nuisance parameter in form of the precision
matrix of the covariates (Zhang & Zhang, 2014; Javanmard & Montanari, 2014;
van de Geer et al., 2014). The construction of a suitable estimating equation
is however somewhat case specific, and can be involved based on the nature of
the high dimensional nuisance parameters. Since SCCA involves a list of high
dimensional nuisance parameters including the covariance matrices
${\Sigma}_{x}$ and ${\Sigma}_{y}$, special attention is required in deriving
our inferential procedures.

Among the above-mentioned methods, our approach bears the greatest resemblance
to the method recently espoused by Janková & van de Geer (2018) in the context
of Sparse Principal Component Analysis (SPCA). However, there are substantial
differences between Janková & van de Geer (2018)’s approach and ours. First,
due to the presence of high dimensional nuisance parameters ${\Sigma}_{x}$ and
${\Sigma}_{y}$, the canonical correlation analysis problem in general is more
complicated than the principal component analysis problem (Gao et al., 2015,
2017). Thus, blindly following Janková & van de Geer (2018) works neither for
the SCCA part nor for the actual de-biasing step. Second, Section 3.1
indicates that the correct choice of the objective function $f$ is crucial to
any de-biasing method. Janková & van de Geer (2018)’s objective function bases
on the well-known fact that the first principal component extraction problem
can be written as an unconstrained Frobenius norm minimization problem. No
such analogue, to the best of our knowledge, was previously available in the
CCA literature. We had to construct a novel objective function whose
unconstrained optimization yields the first canonical directions; see Lemma
3.2. Third, Janková & van de Geer (2018) applies the de-biasing procedure on
some preliminary estimator, similar to us. However, their preliminary
estimators are based on solving a penalized version of the non-convex
principal component analysis optimization problem. To aid the computation, the
authors restrict the search space to a small neighborhood of a consistent
estimator of the first principal component. The said consistent estimator is
found by semi-definite programming. They also show that, any stationary point
of the resulting optimization program consistently estimates the first
principal component. This removes the burden of finding the global minima, but
the program still remains non-convex. Our SCCA method, on the other hand, is
inspired by Gao et al. (2017)’s approach, where the non-convex optimization
part is replaced by a lasso.

##  12 On the Conditions and Assumptions of Section 4

In this section we provide a detailed discussions on assumptions made for the
sake of theoretical developments in Section 4.1.

Discussion on Condition 4.1 First, some remarks are in order regarding the
range of $\kappa\in[1/2,1]$ in Condition 4.1. Theorem 3.2 of Gao et al. (2017)
implies that it is impossible for $\kappa$ to be strictly less than $1/2$
since the minimax rate of the $l_{2}$ error is roughly $s^{1/2}\lambda$ under
Assumption 2.1 and Assumption 2.1. If $\kappa$ is larger, i.e.
$\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$ have slower rates of
convergence, and we pay a price in terms of the sparsity restriction
$s=o(n^{1/(4\kappa)}(\log(p+q))^{-1/(2\kappa)})$ in Assumption 4.1. Supplement
9 shows that estimators satisfying Condition 4.1 with $\kappa=1$ exist. In
fact, most SCCA estimators with theoretical guarantees have $l_{2}$ error
guarantee of $s^{\kappa}\lambda$ with $\kappa\in[1/2,1]$. The interested
reader can refer to Gao et al. (2017, 2015); Chen et al. (2013) and references
therein. Subsequently, in view of the above, we let $\kappa\in[1/2,1]$.

In light of Condition 4.1, indeed $\widehat{\alpha}_{n}$ and
$\widehat{\beta}_{n}$ with faster rate of convergence, i.e. $\kappa=1/2$, is
preferable.  COLAR and Chen et al. (2013)’s estimator attain this minimax rate
when $r=1$. We do not yet know if there are SCCA estimators which attain the
minimax rate for $r>1$ while only estimating the first canonical direction.
For $r>1$, the estimation problem becomes substantially harder because the
remaining $r-1$ canonical directions start acting as high dimensional nuisance
parameters. It is likely that a trade-off between computational and estimation
efficiency arises in presence of these additional nuisance parameters. In
particular, it is plausible that the minimax rate of $\kappa=1/2$ may not be
achievable by polynomial time algorithms in this case. To gather intuition
about this, it is instructive to look at the literature on estimating the
first principal component direction in high dimensions under sparsity. In this
case, to the best of our knowledge, polynomial time algorithms attain the
minimax rate only in the single spike model, or a slightly relaxed version of
the latter. We refer the interested reader to Wang et al. (2016) for more
details.  The algorithms that do succeed to estimate the first principal
component under multiple spikes at the desired minimax rate attempt to solve
the underlying non-convex problem, and hence are not immediately clear to be
polynomial time (Yuan & Zhang, 2013; Ma et al., 2013; Janková & van de Geer,
2018). In this case, Yuan & Zhang (2013) and Ma et al. (2013)’s methods
essentially reduce to power methods that induce sparsity by iterative
thresholding. Chen et al. (2013)’s method tries to borrow this idea in context
of SCCA in the rank one case; see Remark 9.3 for a discussion on the problems
that their method may face in presence of nuisance canonical directions.

Finally for the inferential question, it is natural to consider an extension
of ideas from sparse PCA as developed in (Janková & van de Geer, 2018). When
translated to SCCA, their approach will aim to solve

$\displaystyle\mathmakebox[width("$\underset{\displaystyle
x\in\mathbb{R}^{p},y\in\mathbb{R}^{q}}{\mathrm{subject~{}to}}$")][l]{\underset{\displaystyle
x\in\mathbb{R}^{p},y\in\mathbb{R}^{q}}{\mathrm{minimize}}}\quad\widehat{h}_{n}(x,y)+C\lambda(\|x\|_{1}+\|y\|_{1})\hfil\hfil\hfil\hfil,$
(32)

where $C>0$ is a constant, and

$\widehat{h}_{n}(x,y)=(x^{T}\widehat{\Sigma}_{n,x}x)^{2}/2+(y^{T}\widehat{\Sigma}_{n,y}y)^{2}/2-2x^{T}\widehat{\Sigma}_{n,xy}y.$

We conjecture that for a suitably chosen $C$, the resulting estimators will
satisfy Condition 4.1 with $\kappa=1/2$. However, (32) is non-convex and
solving (32) is computationally challenging for large $p$ and $q$. Analogous
to Janková & van de Geer (2018), one can simplify the problem by searching for
any stationary point of (32) over a smaller feasible set, namely a small
neighborhood of a consistent preliminary estimator of $\alpha_{0}$ and
$\beta_{0}$. However, while this first stage does guarantee a good
initialization, the underlying optimization problem still remains non-convex.
Since the aim of the paper is efficient inference of $\alpha_{0}$ and
$\beta_{0}$ whose computational efficiency is theoretically guaranteed, we
stick with the modified COLAR estimators and refrain from exploring the above-
mentioned route.

Discussion on Assumption 4.1: It is natural to wonder whether the condition
$s^{2\kappa}\lambda^{2}=o(n^{-1/2})$ is at all necessary, especially since it
is much stricter than $s\lambda=o(1)$, which is sufficient for the $l_{2}$
consistency of $\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$ presented in
Theorem 9.1 of Supplement 9. However, current literature on inference in high
dimensional sparse models bears evidence that the restriction
$s\lambda^{2}=o(n^{-1/2})$ might be unavoidable. In fact, this sparsity
requirement is a staple in most de-biasing approaches whose preliminary
estimators are minimax optimal, including sparse principal component analysis
(Janková & van de Geer, 2018) and sparse generalized linear models (van de
Geer et al., 2014; Javanmard & Montanari, 2014). Indeed, in case of sparse
linear regression, Cai et al. (2017) shows that this sparsity is necessary for
adaptive inference. We believe similar results hold for our case as well.
However, further enquiry in that direction is beyond the scope of the present
paper.

Next, it is natural to ask why Assumption 4.1 involves sparsity restriction
not only on $\alpha_{0}$ and $\beta_{0}$, but also on the other columns of $U$
and $V$. This restriction stems from the initial estimation procedure of
$\alpha_{0}$ and $\beta_{0}$. Although we estimate only the first pair of
canonical directions, the remaining canonical directions act as nuisance
parameters. Thus, to efficiently estimate $\alpha_{0}$ and $\beta_{0}$ , we
need to separate the other covariates from $\alpha_{0}$ and $\beta_{0}$.
Therefore, we need to estimate the other covariates’ effect efficiently
enough. Consequently we require some regularity assumptions on these nuisance
parameters as precisely quantified by Assumption 4.1.

Discussion on Condition 4.1: This is a standard assumption in de-biasing
literature in that similar assumptions have appeared in sparse PCA (Janková &
van de Geer, 2018) and sparse generalized linear models literature (van de
Geer et al., 2014) – both of whom use the nodewise lasso algorithm to
construct $\widehat{\Phi}_{n}$. We remark in passing that that Javanmard &
Montanari (2014)’s construction of de-biased lasso does not require the
analogue of $\widehat{\Phi}_{n}$, which is the precision matrix estimator in
their case, to satisfy any condition like Condition 4.1. Instead, it requires
$(\widehat{\Phi}_{n})_{i}^{T}\widehat{\Sigma}_{n,x}(\widehat{\Phi}_{n})_{i}$’s
to be small. It is unknown whether such constructions work in the more
complicated scenario of CCA or PCA.

##  13 Proof Preliminaries

This section states the facts and lemmas that are used repeatedly in the
proofs. The proofs are deferred to Section 22 unless they are very trivial.

First, we derive some results for $\widehat{\alpha}_{n}$ and
$\widehat{\beta}_{n}$ satisfying Condition 4.1.

######  Lemma 13.1.

Suppose $\widehat{\alpha}_{n}$ and $\widehat{\beta}_{n}$ satisfy Condition
4.1. Further suppose Assumption 2.1 and Assumption 4.1 hold. Then

$\|\widehat{\rho}_{n}|-{\rho}_{0}|=O_{p}(s^{\kappa}\lambda).$

Moreover

$\widehat{\alpha}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{\alpha}_{n}-1=O_{p}(s^{\kappa}\lambda)\quad\text{and}\quad\widehat{\beta}_{n}^{T}\widehat{\Sigma}_{n,y}\widehat{\beta}_{n}-1=O_{p}(s^{\kappa}\lambda)$

Recall that we have defined

$\widehat{x}_{n}=|\widehat{\rho}_{n}|^{1/2}\widehat{\alpha}_{n},\quad\widehat{y}_{n}=|\widehat{\rho}_{n}|^{1/2}\widehat{\beta}_{n},\quad
x^{0}=({\rho}_{0})^{1/2}\alpha_{0},\quad y^{0}=({\rho}_{0})^{1/2}\beta_{0}.$

The following lemma gives the rates of $\widehat{x}_{n}$ and $\widehat{y}_{n}$
when Condition 4.1 holds.

######  Lemma 13.2.

Under the set up of Lemma 13.1,

$\inf_{w\in\\{\pm 1\\}}\|w\widehat{x}_{n}-x^{0}\|_{1}+\inf_{w\in\\{\pm
1\\}}\|w\widehat{y}_{n}-y^{0}\|_{1}=O_{p}(s^{\kappa+1/2}\lambda)$

and

$\inf_{w\in\\{\pm 1\\}}\|w\widehat{x}_{n}-x^{0}\|_{2}+\inf_{w\in\\{\pm
1\\}}\|w\widehat{y}_{n}-y^{0}\|_{2}=O_{p}(s^{\kappa}\lambda)$

where $\kappa$ is as defined in (22).

The following lemma entails that
$\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{x}_{n}$ and
$\widehat{y}_{n}^{T}\widehat{\Sigma}_{n,y}\widehat{y}_{n}$ consistently
estimate ${\rho}_{0}^{2}$.

######  Lemma 13.3.

Under the set up of Lemma 13.1, we have

$\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{x}_{n}-{\rho}_{0}=O_{p}(s^{\kappa}\lambda)\quad\widehat{y}_{n}^{T}\widehat{\Sigma}_{n,y}\widehat{y}_{n}-{\rho}_{0}=O_{p}(s^{\kappa}\lambda)$

where $\kappa$ is as defined in(22).

######  Proof 13.4 (of Lemma 13.3).

Noting

$|\widehat{x}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{x}_{n}-{\rho}_{0}|\leq|\widehat{\rho}_{n}(\widehat{\alpha}_{n}^{T}\widehat{\Sigma}_{n,x}\widehat{\alpha}_{n}-1)|+\|\widehat{\rho}_{n}|-{\rho}_{0}|,$

the proof follows from Lemma 13.1 and the fact $|\widehat{\rho}_{n}|\leq 1$.
The proof for $\widehat{y}_{n}$ follows in a similar way.

Now we state an implication of Assumption 4.1.

######  Fact 1.

Suppose $\lambda$ is as in (4). Then Assumption 4.1 implies
$s^{\kappa+1/2}\lambda=o(1)$ and $s\lambda=o(1)$.

Now we state some linear algebra facts.

######  Fact 2.

For any two matrices $A,B\in\mathbb{R}^{p\times q}$, we have

$\|P_{A}-P_{B}\|_{F}^{2}=\text{rank}(A)+\text{rank}(B)-2\text{tr}(P_{A}P_{B}),$

where $P_{A}$ and $P_{B}$ are the projection matrices onto the column spaces
of $A$ and $B$, respectively.

######  Proof 13.5.

Noting $P_{A}^{2}=P_{A}$ and $P_{B}^{2}=P_{B}$, we obtain

$\displaystyle\|P_{A}-P_{B}\|_{F}^{2}=\text{tr}\Big{(}(P_{A}-P_{B})^{T}(P_{A}-P_{B})\Big{)}=\text{tr}(P_{A})+\text{tr}(P_{B})-2\text{tr}(P_{A}P_{B}),$

from which the result follows because for projection matrix $P_{A}$,
$\text{tr}(P_{A})=\text{rank}(A)$.

######  Fact 3.

For any matrix $A\in\mathbb{R}^{r\times r}$, $\|A\|_{F}\leq r^{1/2}\|A\|_{op}$

######  Proof 13.6 (of Fact 3).

Suppose $\varsigma_{i}$’s are the singular values of $A$. Then

$\|A\|_{F}^{2}=\sum_{i=1}^{r}\varsigma_{i}^{2}\leq r\max_{1\leq i\leq
r}\varsigma_{i}^{2}.$

Therefore $\|A\|_{F}\leq r^{1/2}\|A\|_{op}$.

######  Fact 4 (Lemma 2.1.3 of Chen et al. (2020)).

Suppose $A$ and $B$ are two matrices in $\mathbb{R}^{p\times q}$. Then

$2^{-1/2}\|P_{A}-P_{B}\|_{F}\leq\inf_{W\in\mathcal{O}(r,r)}\|AW-B\|_{F}\leq\|P_{A}-P_{B}\|_{F}$

######  Fact 5.

Suppose $x$ and $y\in\mathbb{R}^{p}$. Then

$\|P_{x}-P_{y}\|_{F}\leq 4\inf_{w\in\\{\pm
1\\}}\frac{\|wx-y\|_{2}}{\max(\|x\|_{2},\|y\|_{2})}.$

The next lemma shows that the $l_{1}$ and $l_{2}$ norms of $\alpha_{0}$ and
$\beta_{0}$ are bounded.

######  Lemma 13.7.

Under Assumption 2.1, we have

$\|\alpha_{0}\|_{1},\|\beta_{0}\|_{1}\leq(Ms)^{1/2},$

where $s=s_{U}+s_{V}$. Also

$\|\alpha_{0}\|_{2},\|\beta_{0}\|_{2}\leq M^{1/2},$

where $M$ is as in Assumption 2.1.

######  Proof 13.8.

Since
$\|\alpha_{0}\|^{2}_{2}\Lambda_{min}({\Sigma}_{x})\leq|\alpha_{0}^{T}{\Sigma}_{x}\alpha_{0}|$,
we have $\|\alpha_{0}\|_{2}\leq\sqrt{M}$ by the 2.1 Assumption. Similarly,
$\|\beta_{0}\|_{2}\leq\sqrt{M}.$ Now, Cauchy Schwartz inequality implies

$\|\alpha_{0}\|_{1}\leq\sqrt{s}\|\alpha_{0}\|_{2}=\sqrt{Ms}.$

The same can be proved for $\beta_{0}$, which completes the proof.

Now we state some rate-results for Sub-Gaussian covariance matrices.

######  Lemma 13.9.

Suppose $X\in\mathbb{R}^{n\times p}$ and $Y\in\mathbb{R}^{n\times q}$ are sub-
Gaussian matrices. Then there exists constant $C$ depending only on the
subgaussian parameter of $(X,Y)$ so that

$|\widehat{\Sigma}_{n,x}-{\Sigma}_{x}|_{\infty},|\widehat{\Sigma}_{n,y}-{\Sigma}_{y}|_{\infty},|\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy}|_{\infty}\leq
C\lambda$

with high probability as $n,p,q\to\infty$. Moreover, for any
$v\in\mathbb{R}^{p}$, there exists constant $C$ depending only on the
subgaussian parameter of $(X,Y)$ so that

$\|(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})v\|_{\infty}\leq
C\|v\|_{2}\lambda\quad\text{and}\quad\|(\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy})v\|_{\infty}\leq
C\|v\|_{2}\lambda.$

with high probability as $n,p,q\to\infty$.

######  Lemma 13.10.

Let $v\in\mathbb{R}^{p}$ and the set $S\subset\\{1,\ldots,p\\}$ has
cardinality $s$. Suppose $v$ satisfies the cone condition
$\|v_{S^{c}}\|_{1}\leq C^{\prime}\|v_{S}\|_{1}$ for some constant
$C^{\prime}>0$ and $S$ has cordinality $s$. Then under the set up of Lemma
13.9, there exists $C>0$ depending only on the subgaussian parameter of $X$
and $C^{\prime}$ so that

$|v^{T}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})v|\leq Cs\|v\|_{2}^{2}\lambda,$

with high probability as $n,p,q\to\infty$.

######  Proof 13.11.

Noting
$\|v\|_{1}\leq(C^{\prime}+1)\|v_{S}\|_{1}\leq(C^{\prime}+1)s^{1/2}\|v_{S}\|_{2}$,
we obtain

$\absolutevalue{v^{T}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})v}\leq\|v\|_{1}^{2}|\widehat{\Sigma}_{n,x}-{\Sigma}_{x}|_{\infty}\leq(C^{\prime}+1)^{2}s\|v\|_{2}^{2}|\widehat{\Sigma}_{n,x}-{\Sigma}_{x}|_{\infty}.$

Then the proof follows by Lemma 13.9.

######  Lemma 13.12.

Supppose $X$ is sub-Gaussian and a random vector
$\widehat{z}_{n}\in\mathbb{R}^{p}$ satisfies
$\|\widehat{z}_{n}\|_{1}=O_{p}(s^{1/2})$, where $s$ satisfies
$\lambda^{1/2}s=o(1)$. Then there exists $C$ depending only on the sub-
gaussian parameters of $X$ so that

$|\widehat{z}_{n}^{T}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})\widehat{z}_{n}|\leq
C(s^{1/2}\lambda\|\widehat{z}_{n}\|_{2}^{2}+\lambda\|\widehat{z}_{n}\|_{1})$

with high probability for sufficiently large $n,p$, and $q$.

######  Lemma 13.13.

Suppose $X$ and $Y$ are sub-Gaussian and random vectors
$\widehat{z}_{n}\in\mathbb{R}^{p},\widehat{w}_{n}\in\mathbb{R}^{q}$ satisfy

$\|\widehat{z}_{n}\|_{1}+\|\widehat{w}_{n}\|_{1}=O_{p}(s^{1/2}),$

where $s$ satisfies Assumption 4.1. Then there exists $C$ depending only on
the subgaussian parameters of $X$ and $Y$ so that

$|\widehat{z}_{n}^{T}(\widehat{\Sigma}_{n,xy}-{\Sigma}_{xy})\widehat{w}_{n}|\leq
C\lambda\Big{(}s^{1/2}(\|\widehat{z}_{n}\|_{2}^{2}+\|\widehat{w}_{n}\|_{2}^{2})+(\|\widehat{z}_{n}\|_{1}+\|\widehat{w}_{n}\|_{1})\Big{)}$

with high probability as $n,p,q\to\infty$.

######  Lemma 13.14 (Lemma 8 of Janková & van de Geer).

Suppose $X$ is sub-Gaussian and $z\in\mathbb{R}^{p}$ is a vector with
$\|z\|_{0}=s$. Then

$\sup_{z\in\mathbb{R}^{p}}\frac{z^{T}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})z}{\|z\|^{2}_{2}}=O_{p}((s\log
p/n)^{1/2}).$

######  Lemma 13.15.

Suppose $\widehat{z}_{n}$ is a random vector, possibly depending on $X$, so
that $\|\widehat{z}_{n}-z_{0}\|_{1}=o_{p}(1)$ where $z_{0}$ is a fixed vector
with finite $l_{2}$ norm. Then depending only on the subgaussian parameter of
$X$ so that

$|x^{T}(\widehat{\Sigma}_{n,x}-{\Sigma}_{x})\widehat{z}_{n}|\leq
C\|z_{0}\|_{2}\|x\|_{1}O_{p}(\lambda)$

with high probability as $n,p,q\to\infty$.

Our next result is on multivariate normal distribution. This result quite well
known and can be obtained via straightforward calculation.

######  Fact 6 (Fourth moments of multivariate normal distribution).

Suppose $X\sim N_{p}(0,{\Sigma}_{x})$. Then