Sparse CCA: Adaptive Estimation and Computational Barriers ∗
Chao Gao1, Zongming Ma2 , and Harrison H. Zhou1

arXiv:1409.8565v4 [stat.ME] 4 Apr 2016

1

2

Yale University
University of Pennsylvania

Abstract
Canonical correlation analysis is a classical technique for exploring the relationship
between two sets of variables. It has important applications in analyzing high dimensional datasets originated from genomics, imaging and other fields. This paper considers
adaptive minimax and computationally tractable estimation of leading sparse canonical
coefficient vectors in high dimensions. First, we establish separate minimax estimation
rates for canonical coefficient vectors of each set of random variables under no structural
assumption on marginal covariance matrices. Second, we propose a computationally feasible estimator to attain the optimal rates adaptively under an additional sample size
condition. Finally, we show that a sample size condition of this kind is needed for any
randomized polynomial-time estimator to be consistent, assuming hardness of certain
instances of the Planted Clique detection problem. The result is faithful to the Gaussian models used in the paper. As a byproduct, we obtain the first computational lower
bounds for sparse PCA under the Gaussian single spiked covariance model.
Keywords. Convex programming, group-Lasso, Minimax rates, Computational complexity, Planted Clique, Sparse CCA (SCCA), Sparse PCA (SPCA)

1

Introduction

Canonical correlation analysis (CCA) [23] is a classical and important tool in multivariate
statistics [1, 31]. For two random vectors X ∈ Rp and Y ∈ Rm , at the population level, CCA
finds successive vectors uj ∈ Rp and vj ∈ Rm (called canonical coefficient vectors) that solve
max a′ Σxy b,
a,b

subject to a′ Σx a = b′ Σy b = 1, a′ Σx ul = b′ Σy vl = 0, ∀0 ≤ l ≤ j − 1,

(1)

where Σx = Cov(X), Σy = Cov(Y ), Σxy = Cov(X, Y ), u0 = 0, and v0 = 0. Since our primary
interest lies in the covariance structure among X and Y , we assume that their means are
zeros from here on. Then the linear combinations (u′j X, vj′ Y ) are the j-th pair of canonical
variates. This technique has been widely used in various scientific fields to explore the
relationship between two sets of variables. In practice, one does not have knowledge about
∗

The research of C. Gao and H. H. Zhou is supported in part by NSF Career Award DMS-0645676 and NSF
FRG Grant DMS-0854975. The research of Z. Ma is supported in part by NSF Career Award DMS-1352060.

1

b x,
the population covariance, and Σx , Σy , and Σxy are replaced by their sample versions Σ
b y , and Σ
b xy in (1).
Σ
Recently, there have been growing interests in applying CCA to analyzing high-dimensional
datasets, where the dimensions p and m could be much larger than the sample size n. It has
been well understood that classical CCA breaks down in this regime [25, 4, 18]. Motivated
by genomics, neuroimaging and other applications, people have become interested in seeking
sparse leading canonical coefficient vectors. Various estimation procedures imposing sparsity on canonical coefficient vectors have been developed in the literature, which are usually
termed sparse CCA. See, for example, [42, 43, 33, 22, 27, 38, 3].
The theoretical aspect of sparse CCA has also been investigated in the literature. A useful
model for studying sparse CCA is the canonical pair model proposed in [14]. In particular,
suppose there are r pairs of canonical coefficient vectors (and canonical variates) among the
two sets of variables, then the model reparameterizes the cross-covariance matrix as
Σxy = Σx U ΛV ′ Σy ,

where U ′ Σx U = V ′ Σy V = Ir .

(2)

Here U = [u1 , ..., ur ] and V = [v1 , ..., vr ] collect the canonical coefficient vectors and Λ =
diag(λ1 , . . . , λr ) with 1 > λ1 ≥ · · · ≥ λr > 0 are the ordered canonical correlations. Let
Su = supp(U ) and Sv = supp(V ) be the indices of nonzero rows of U and V . One way to
impose sparsity on the canonical coefficient vectors is to require the sizes of Su and Sv to be
small, namely |Su | ≤ su and |Sv | ≤ sv for some su ≤ p and sv ≤ m. Under this model, Gao
et al. [18] showed that the minimax rate for estimating U and V under the joint loss function
b Vb ′ − U V ′ k2 is
kU
F
ep
em 
1 
r(s
+
s
)
+
s
log
+
s
log
.
(3)
u
v
u
v
nλ2r
su
sv
However, to achieve the rate, Gao et al. [18] used a computationally infeasible and nonadaptive procedure, which requires exhaustive search of all possible subsets with the given
cardinality and the knowledge of su and sv . Moreover, it is unclear from (3) whether the
estimation error of U depends on the sparsity and the ambient dimension of V and vice versa.
The goal of the present paper is to study three fundamental questions in sparse CCA: (1)
What are the minimax rates for estimating the canonical coefficient vectors on the two sets
of variables separately? (2) Is there a computationally efficient and sparsity-adaptive method
that achieves the optimal rates? (3) What is the price one has to pay to achieve the optimal
rates in a computationally efficient way?

1.1

Main contributions

We now introduce the main contributions of the present paper from three different viewpoints
as suggested by the three questions we have raised.
b Vb ′ − U V ′ k2 studied by [18] characterizes the
Separate minimax rates The joint loss kU
F
joint estimation error of both U and V . In this paper, we provide a finer analysis by studying
individual estimation errors of U and V under a natural loss function that can be interpreted
as prediction error of canonical variates. The exact definition of the loss functions is given in
Section 2. Separate minimax rates are obtained for U and V . In particular, we show that the
minimax rate of convergence in estimating U depends only on n, r, λr , p and su , but not on
2

either m or sv . Consequently, if U is sparser than V , then convergence rate for estimating U
can be faster than that for estimating V . Such a difference is not reflected by the joint loss,
since its minimax rate (3) is determined by the slower of the rates of estimating U and V .
Adaptive estimation As pointed out in [14] and [18], sparse CCA is a more difficult
problem than the well-studied sparse PCA. A naive application of sparse PCA algorithm to
sparse CCA can lead to inconsistent results [14]. The additional difficulty in sparse CCA
mainly comes from the presence of the nuisance parameters Σx and Σy , which cannot be
estimated consistently in a high-dimensional regime in general. Therefore, our goal is to
design an estimator that is adaptive to both the nuisance parameters and the sparsity levels.
Under the canonical pair model, we propose a computationally efficient algorithm. The
algorithm has two stages. In the first stage, we propose a convex program for sparse CCA
based on a tight convex relaxation of a combinatorial program in [18] by considering the
smallest convex set containing all matrices of the form AB ′ with both A and B being rankr orthogonal matrices. The convex program can be efficiently solved by the Alternating
Direction Method with Multipliers (ADMM) [16, 11]. Based on the output of the first stage,
we formulate a sparse linear regression problem in the second stage to improve estimation
b and Vb can be obtained via a group-Lasso algorithm [45].
accuracy, and the final estimator U
Under the sample size condition that
n ≥ Csu sv log(p + m)/λ2r

(4)

b and Vb recover the true canonical
for some sufficiently large constant C > 0, we show U
coefficient matrices U and V within optimal error rates adaptively with high probability.
Computational lower bound We require the sample size condition (4) for the adaptive
procedure to achieve optimal rates of convergence. Assuming hardness of certain instances of
the Planted Clique detection problem, we provide a computational lower bound to show that a
condition of this kind is unavoidable for any computationally feasible estimation procedure to
achieve consistency. Up to an asymptotically equivalent discretization which is necessary for
computational complexity to be well-defined, our computational lower bound is established
directly for the Gaussian canonical pair model used throughout the paper.
An analogous sample size condition has been imposed in the sparse PCA literature [24,
29, 13, 37], namely n ≥ Cs2 log p/λ2 where s is the sparsity of the leading eigenvector
and λ the gap between the leading eigenvalue and the rest of the spectrum. Berthet and
Rigollet [7] showed that if there existed a polynomial-time algorithm for a generalized sparse
PCA detection problem while such a condition is violated, the algorithm could be made (in
randomized polynomial-time) into a detection method for the Planted Clique problem in a
regime where it is believed to be computationally intractable. However, both the null and the
alternative hypotheses in the sparse PCA detection problem were generalized in [7] to include
all multivariate distributions whose quadratic forms satisfy certain uniform tail probability
bounds and so the distributions need not be Gaussian or having a spiked covariance structure
[24]. The same remark also applies to the subsequent work on sparse PCA estimation [39].
Hence, the computational lower bound in sparse PCA was only established for such enlarged
parameter spaces. As a byproduct of our analysis, we establish the desired computational
lower bound for sparse PCA in the Gaussian single spiked covariance model.
3

1.2

Organization

After an introduction to notation, the rest of the paper is organized as follows. In Section
2, we formulate the sparse CCA problem by defining its parameter space and loss function.
Section 3 presents separate minimax rates for estimating U and V . Section 4 proposes a
two-stage adaptive estimator that is shown to be minimax rate optimal under an additional
sample size condition. Section 5 shows a condition of this kind is necessary for all randomized
polynomial-time estimator to achieve consistency by establishing new computational lower
bounds for sparse PCA and sparse CCA. Section 6 presents proofs of theoretical results in
Section 4. Implementation details of the adaptive procedure, numerical studies, additional
proofs and technical details are deferred to the supplement [19].

1.3

Notation

For any t ∈ Z+ , [t] denotes the set {1, 2, ..., t}. For any set S, |S| denotes its cardinality
and S c its complement. For any event E, 1{E} denotes its indicator function. For any
a, b ∈ R, ⌈a⌉ denotes the smallest integer no smaller than a, ⌊a⌋ the largest integer
qP no
2
larger than a, a ∨ b = max(a, b) and a ∧ b = min(a, b). For a vector u, ||u|| =
i ui ,
P
P
||u||0 = i 1{ui 6=0} , and ||u||1 = i |ui |. For any matrix A = (aij ) ∈ Rp×k , Ai· denotes its
i-th row and supp(A) = {i ∈ [p] : kAi· k > 0}, the index set of nonzero rows, is called its
support. For any subset J ⊂ [p] × [k], AJ = (aij 1{(i,j)∈J} ) ∈ Rp×k is obtained by keeping all
entries in J and replacing all entries in J c with zeros. We write AJ1 J2 for AJ1 ×J2 and A(J1 J2 )c
for A(J1 ×J2 )c . Note that AJ1 ∗ = AJ1 ×[k] ∈ Rp×k while AJ1 · stands for the corresponding
nonzero submatrix of size |J1 | × k. In addition, PA ∈ Rp×p stands for the projection matrix
onto the column space of A, O(p, k) denotes the set of all p × k orthogonal matrices and
O(k) = O(k, k). Furthermore, σi (A) stands for the i-th largest singular value of A and
σmax (A) = σq
1 (A), σmin (A) = σp∧k (A). The Frobenius norm and the operator norm of A
P
a2 and kAkop = σ1 (A), respectively. The l1 norm and the nuclear norm
are kAkF =
P i,j ij
P
are ||A||1 =P ij |aij | and kAk∗ = i σi (A), respectively. If A is a square matrix, its trace
is Tr(A) = i aii . For two square matrices A and B, we write A  B if B − A is positive
semidefinite. For any positive semi-definite matrix A, A1/2 denotes its principal square root
that is positive semi-definite and satisfies A1/2 A1/2 = A. The trace inner product of two
matrices A, B ∈ Rp×k is hA, Bi = Tr(A′ B). Given a random element X, L(X) denotes its
probability distribution. The symbol C and its variants C1 , C ′ , etc. are generic constants
and may vary from line to line, unless otherwise specified. The symbols P and E stand for
generic probability and expectation when the distribution is clear from the context.

2
2.1

Problem Formulation
Parameter space

Consider a canonical pair model where the observed pairs of measurement vectors (Xi′ , Yi′ )′ ,
i = 1, . . . , n are i.i.d. from a multivariate Gaussian distribution Np+m (0, Σ) where


Σx Σxy
,
Σ=
Σyx Σy
4

with the cross-covariance matrix Σxy satisfying (2). We are interested in the situation where
the leading canonical coefficient vectors are sparse. One way to quantify the level of sparsity
is to bound how many nonzero rows there are in the U and V matrices. This notion of
sparsity has been used previously in both sparse PCA [13, 37] and sparse CCA [18] problems
when one seeks multiple sparse vectors simultaneously.
Recall that for any matrix A, supp(A) collects the indices of nonzero rows in A. Adopting the above notion of sparsity, we define F(su , sv , p, m, r, λ; M ) to be the collection of all
covariance matrices Σ with the structure (2) satisfying
1. U ∈ Rp×r and V ∈ Rm×r with |supp(U )| ≤ su and |supp(V )| ≤ sv ;

2. σmin (Σx ) ∧ σmin (Σy ) ≥ M −1 and σmax (Σx ) ∨ σmax (Σy ) ≤ M ;

(5)

M −1 .

3. λr ≥ λ and λ1 ≤ 1 −

The probability space we consider is

P(n, su , sv , p, m, r, λ; M ) = L((X1′ , Y1′ )′ , ..., (Xn′ , Yn′ )′ ) :
iid

(Xi′ , Yi′ )′ ∼ Np+m (0, Σ)

with Σ ∈ F(su , sv , p, m, r, λ; M ) ,

(6)

where n is the sample size. We shall allow su , sv , p, m, r, λ to vary with n, while M > 1 is
restricted to be an absolute constant.

2.2

Prediction loss

From now on, the presentation of definitions and results will focus on U only since those
b = [b
for V can be obtained via symmetry. Given an estimator U
u1 , . . . , u
br ] of the leading
canonical coefficient vectors for X, a natural way of assessing its quality is to see how well
it predicts the values of the canonical variables U ′ X ⋆ ∈ Rr for a new observation X ⋆ which
b . This
is independent of and identically distributed as the training sample used to obtain U
leads us to consider the following loss function
b, U) =
L(U

b ′ X ⋆ − U ′ X ⋆ k2 ,
E⋆ kW ′ U

(7)

b W − U )′ Σx (U
b W − U )].
Tr[(U

(8)

inf

W ∈O(r)

b , U ) is still a random quantity
where E⋆ means taking expectation only over X ⋆ and so L(U
b . Since L(U
b , U ) is the expected squared error for predicting the
due to the randomness of U
′
⋆
′
⋆
b
canonical variables U X via U X , we refer to it as prediction loss from now on. It is worth
noting that the introduction of an r × r orthogonal matrix W is unavoidable. To see this,
we can simply consider the case where λ1 = · · · = λr = λ in (2), then we can replace the
pair (U, V ) in (2) by (U W, V W ) for any W ∈ O(r). In other words, the canonical coefficient
vectors are only determined up to a joint orthogonal transform. If we work out the E⋆ part
in the definition (7), then the loss function can be equivalently defined as
b, U) =
L(U

inf

W ∈O(r)

b , X ⋆ and Σx in (7) and (8)
By symmetry, we can define L(Vb , V ) by simply replacing U , U
with V , Vb , Y ⋆ and Σy .
5

A related loss function is kPUb −PU k2F measuring the difference between two subspaces. By
b , U ) is a stronger
Proposition 9.2 in the supplementary material [19], the prediction loss L(U
2
b
loss function. That is, kPUb − PU kF ≤ CL(U , U ) for some constant C > 0 only depending on
b , U ) is strictly stronger. To see this, let Σx = Ip , U ∈ O(p, r) and U
b = 2U .
M . Actually, L(U
2
b , U ) = inf W ∈O(r) kΣ1/2
b
Then, kPUb − PU k2F = 0, while L(U
x (U W − U )kF = inf W ∈O(r) (5r −
b , U ), and provide brief
Tr(W )) = r > 0. In this paper, we will focus on the stronger loss L(U
2
remarks on results for kPUb − PU kF .

3

Minimax Rates

We first provide a minimax upper bound using a combinatorial optimization procedure, and
then show that the resulting rate is optimal by further providing a matching minimax lower
bound.
Let (Xi′ , Yi′ )′ ∈ Rp+m, i = 1, . . . , n, be i.i.d. observations following Np+m (0, Σ) for some
Σ ∈ F(su , sv , p, m, r, λ; M ). For notational convenience, we assume the sample size is divisible
by three, i.e., n = 3n0 for some n0 ∈ N.
Procedure To obtain minimax upper bound, we propose a two-stage combinatorial op0
timization procedure. We split the data into three equal size batches D0 = {(Xi′ , Yi′ )′ }ni=1
,
2n0
′
′
′
′
′
′
n
D1 = {(Xi , Yi ) }i=n0 +1 and D2 = {(Xi , Yi ) }i=2n0 +1 , and denote the sample covariance mab (j)
b (j)
b (j)
trices computed on each batch by Σ
x , Σy and Σxy for j ∈ {0, 1, 2}.
b (0) , Vb (0) ) which solves the following program:
In the first stage, we find (U
max

L∈Rp×r ,R∈Rm×r

b (0) R),
Tr(L′ Σ
xy

b (0) L = R′ Σ
b (0) R = Ir , and
subject to L′ Σ
x
y

(9)

|supp(L)| ≤ su , |supp(R)| ≤ sv .

b (1) solving
In the second stage, we further refine the estimator for U by finding U
min

L∈Rp×r

′ b (1) b (0)
b (1)
Tr(L′ Σ
)
x L) − 2 Tr(L Σxy V

(10)

subject to |supp(L)| ≤ su .

b (1) , defined as
The final estimator is a normalized version of U
b =U
b (1) ((U
b (1) )′ Σ
b (2)
b (1) )−1/2 .
U
x U

(11)

The purpose of sample splitting employed in the above procedure is to facilitate the proof.
Theory and discussion The program (9) was first proposed in [18] as a sparsity constrained version of the classical CCA formulation. However, the resulting estimator will have
a convergence rate that involves the sparsity level sv and the ambient dimension m of the V

6

matrix [18, Theorem 1], which is sub-optimal. The second stage in the procedure is thus proposed to further pursue the optimal estimation rates. First, if we were given the knowledge
of V , then the least square solution of regressing V ′ Y ∈ Rr on X ∈ Rp is
U Λ = argmin EkY ′ V − X ′ Lk2
L∈Rp×r

= argmin Tr(L′ Σx L) − 2 Tr(L′ Σxy V ) + Tr(V ′ Σy V )
L∈Rp×r

(12)

= argmin Tr(L′ Σx L) − 2 Tr(L′ Σxy V ),
L∈Rp×r

where the expectation is with respect to the distribution (X ′ , Y ′ )′ ∼ Np+m (0, Σ). The second
equality results from taking expectation over each of the three terms in the expansion of the
square Euclidean norm, and the last equality holds since Tr(V ′ Σy V ) does not involve the
argument to be optimized over. In fact, from the canonical pair model, one can easily derive
a regression interpretation of CCA, V ′ Y = ΛU ′ X + E, where E ∼ N (0, Ir − Λ2 ). Then, (10)
is a least square formulation of the regression interpretation. However, CCA is different from
regression because the response V ′ Y depends on an unknown V . Comparing (10) with (12),
it is clear that (10) is a sparsity constrained version of (12) where the knowledge of V and the
covariance matrix Σ are replaced by the initial estimator Vb (0) and sample covariance matrix
b (1) can be viewed as an estimator of U Λ. Hence,
from an independent sample. Therefore, U
a final normalization step is taken in (11) to transform it to an estimator of U .
We now state a bound for the final estimator (11).
Theorem 3.1. Assume
ep
em 
1
r(su + sv ) + su log
+ sv log
≤c
n
su
sv

(13)

for some sufficiently small constant c > 0. Then there exist constants C, C ′ > 0 only depending on c such that


b , U ) ≤ C su r + log ep ,
L(U
(14)
nλ2
su

with P-probability at least 1 − exp (−C ′ (su + log(ep/su ))) − exp (−C ′ (sv + log(em/sv ))) uniformly over P ∈ P(n, su , sv , p, m, r, λ; M ).

Remark 3.1. The paper assumes that M is a constant. However, it is worth noting that the
minimax upper bound of Theorem 3.1 does not depend on M even if M is allowed to grow
with n. To be specific, assume the eigenvalues of Σx are bounded
in the interval [M1 , M2 ].

b , U ) would still be 1 2 su r + log ep , because the dependence on
The convergence rate of L(U
su
nλ
M1 , M2 has been implicitly built into the prediction
loss.
On
the other hand, a convergence
 1
ep 
M2
2
rate for the loss kPUb − PU kF would be M1 nλ2 su r + log su , with an extra factor of the
condition number of Σx .
Under assumption (13), Theorem 3.1 achieves a convergence rate for the prediction loss
in U that does not depend on any parameter related to V . Note that the probability tail still
′
involves m and sv . However, it can be shown that exp (−C ′ (sv + log(em/sv ))) ≤ m−C /2 ,
and so the corresponding term in the tail probability goes to 0 as long as m → ∞. The
optimality of this upper bound can be justified by the following minimax lower bound.
7

v
Theorem 3.2. Assume that r ≤ su ∧s
2 . Then there exists some constant C > 0 only depending on M and an absolute constant c0 > 0, such that



C
ep 
b
inf sup P L(U , U ) ≥ c0 ∧
su r + log
≥ 0.8,
b P∈P
nλ2
su
U

where P = P(n, su , sv , p, m, r, λ; M ).

By Theorem 3.1 and Theorem 3.2, the rate in (14), whenever it is upper bounded by a
constant, is the minimax rate of the problem.

4

Adaptive and Computationally Efficient Estimation

Section 3 determines the minimax rates for estimating U under the prediction loss. However,
there are two drawbacks of the procedure (9) – (11). One is that it requires the knowledge
of the sparsity levels su and sv . It is thus not adaptive. The other is that in both stages
one needs to conduct exhaustive search over all subsets of given sizes in the optimization
problems (9) and (10), and hence the computation cost is formidable.
In this section, we overcome both drawbacks by proposing a two-stage convex program
approach towards sparse CCA. The procedure is named CoLaR, standing for Convex program
with group-Lasso Refinement. It is not only computationally feasible but also achieves the
minimax estimation error rates adaptively over a large collection of parameter spaces under an
additional sample size condition. The issues related to this additional sample size condition
will be discussed in more detail in the subsequent Section 5.

4.1

Estimation scheme

The basic principle underlying the computationally feasible estimation scheme is to seek tight
convex relaxations of the combinatorial programs (9) – (10). In what follows, we introduce
convex relaxations for the two stages in order. As in Section 3, we assume that the data is
b (j)
b (j)
b (j)
split into three batches D0 , D1 and D2 of equal sizes and for j = 0, 1, 2, let Σ
x , Σy and Σxy
be defined as before.
First stage By the definition of trace inner product, the objective function in (9) can be
b xy R) = hΣ
b xy , LR′ i. Since it is linear in F = LR′ , this suggests treating
rewritten as Tr(L′ Σ
′
LR as a single argument rather than optimizing over L and R separately. Next, the support
size constraints |supp(L)| ≤ su , |supp(R)| ≤ sv imply that the vector ℓ0 norm kLR′ k0 ≤ su sv .
Applying the convex relaxation of ℓ0 norm by ℓ1 norm and including it as a Lagrangian term,
we are led to consider a new objective function
b (0) , F i − ρ||F ||1 ,
max hΣ
xy

F ∈Rp×m

(15)

P
where F serves as a surrogate for LR′ , kF k1 = i∈[p],j∈[m] |Fij | denotes the vector ℓ1 norm
of the matrix argument, and ρ is a penalty parameter controlling sparsity. Note that (15)
is the maximization problem of a concave function, which becomes a convex program if the

8

constraint set is convex. Under the identity F = LR′ , the normalization constraint in (9)
reduces to
1/2
1/2
b (0)
b (0)
(Σ
F (Σ
∈ Or = {AB ′ : A ∈ O(p, r), B ∈ O(m, r)}.
x )
y )

(16)

Cr = {G ∈ Rp×m : kGk∗ ≤ r, kGkop ≤ 1} = conv(Or )

(17)

1/2 F (Σ
1/2 ∈ C where
b (0)
b (0)
Naturally, we relax it to (Σ
x )
y )
r

is the smallest convex set containing Or . The relation (17) is stated in the proof of Theorem
3 of [40]. Combining (15) – (17), we use the following convex program for the first stage in
our adaptive estimation scheme:
max

F ∈Rp×m

b (0) , F i − ρ||F ||1
hΣ
xy

b (0) )1/2 F (Σ
b (0) )1/2 k∗ ≤ r, k(Σ
b (0) )1/2 F (Σ
b (0) )1/2 kop ≤ 1.
subject to k(Σ
x
y
x
y

(18)

Implementation of (18) is discussed in Section 10 in the supplement [19].

Remark 4.1. A related but different convex relaxation was proposed in [37] for the sparse
PCA problem, where the set of all rank r projection matrices (which are symmetric) is relaxed
to its convex hull – the Fantope {P : Tr(P ) = r, 0  P  Ip }. Such an idea is not directly
applicable in the current setting due to the asymmetric nature of the matrices included in
the set Or in (16).
Remark 4.2. The risk of the solution to (18) for estimating U V ′ is sub-optimal compared
to the optimal rates determined in [18]. See Theorem 4.1 below. Nonetheless, it leads to
a reasonable estimator for the subspaces spanned by the first r left and right canonical
coefficient vectors under a sample size condition, which is sufficient for achieving the optimal
estimation rates for U and V in a further refinement stage to be introduced below. Although
it is possible that some other penalty function rather than the ℓ1 penalty in (18) could also
achieve this goal, ℓ1 is appealing due to its simplicity.
Second stage Now we turn to the convex relaxation to (10) in the second stage. By the
discussion following Theorem 3.1, if we view the rows of L as groups, then (10) becomes
a least square problem with a constrained number of active groups. A well-known convex
relaxation for such problems is the group-Lasso [45], where the number of active groups
constraint is relaxed by bounding the sum of ℓ2 norms of the coefficient vector of each group.
b (0) (resp. Vb (0) ) be the matrix consisting of its first r left
b be the solution to (18) and U
Let A
(resp. right) singular vectors. Thus, in the second stage of the adaptive estimation scheme,
we propose to solve the following group-Lasso problem:
b (1) L) − 2 Tr(L′ Σ
b (1) Vb (0) ) + ρu
min Tr(L′ Σ
x
xy

L∈Rp×m

p
X
j=1

kLj· k,

(19)

P
where pj=1 kLj· k is the group sparsity penalty, defined as the sum of the ℓ2 norms of all the
row vectors in L, and ρu is a penalty parameter controlling sparsity. Note that the group
sparsity penalty is crucial, since if one uses an ℓ1 penalty instead, only a sub-optimal rate can
9

b (1) , then our final estimator in the adaptive
be achieved. Suppose the solution to (19) is U
estimation scheme is its normalized version
b =U
b (1) ((U
b (1) )′ Σ
b (2)
b (1) )−1/2 .
U
x U

(20)

As before, sample splitting is only used for technical arguments in the proof. Simulation
results in Section 11 in the supplement [19] show that using the whole dataset repeatedly
in (18) – (20) yields satisfactory performance and the improvement by the second stage is
considerable.

4.2

Theoretical guarantees

b to the convex program (18).
We first state the upper bound for the solution A
Theorem 4.1. Assume that

su sv log(p + m)
,
(21)
λ2
for some sufficiently large constant C1 > 0. Then there p
exist positive constants γ1 , γ2 and
′
C, C only depending on M and C1 , such that when ρ = γ log(p + m)/n for γ ∈ [γ1 , γ2 ],
n ≥ C1

b − U V ′ k2 ≤ C
kA
F

su sv log(p + m)
,
nλ2

with P-probability at least 1 − exp(−C ′ (su + log(ep/su ))) − exp(−C ′ (sv + log(em/sv ))) for any
P ∈ P(n, su , sv , p, m, r, λ; M ).
Note that the error bound in Theorem 4.1 can be much larger than the optimal rate for
joint estimation of U V ′ established in [18]. Nonetheless, under the sample size condition (21),
b is close to U V ′ in Frobenius norm distance. This fact, together with
it still ensures that A
the proposed refinement scheme (19) – (20), guarantees the optimal rates of convergence for
the estimator (20) as stated in the following theorem.
Theorem 4.2. Assume (21) holds for some sufficiently large C1 ≥ 0. p
Then there exist
′ [log(p + m)]/n
constants γ and
γ
only
depending
on
C
and
M
such
that
if
we
set
ρ
=
γ
u
1
p
and ρu = γu′ (r + log p)/n for any γ ′ ∈ [γ, C2 γ] and γu′ ∈ [γu , C2 γu ] for some absolute
constant C2 > 0, there exist a constants C, C ′ > 0 only depending on C1 , C2 and M , such
that
b, U) ≤ C
L(U

su (r + log p)
,
nλ2

with P-probability at least 1 − exp(−C ′ (su + log(ep/su ))) − exp(−C ′ (sv + log(em/sv ))) −
exp(−C ′ (r + log(p ∧ m))) uniformly over P ∈ P(n, su , sv , p, m, r, λ; M ).
Remark 4.3. The result of Theorem 4.2 assumes a constant M . Explicit dependence on the
eigenvalues of the marginal covariance can be tracked even when M is diverging. Assuming
b, U)
the eigenvalues of Σx all lie in the interval [M1 , M2 ], then the convergence rate of L(U


2
3
su (r+log p)
su (r+log p)
M2
M2
and a convergence rate of kPUb − PU k2F would be M
.
would be M
nλ2
nλ2
1
1

M2 2
Compared with Remark 3.1, there is an extra factor M1 , which is also present for the Lasso
error bounds [8, 32]. Evidence has been given in the literature that such an extra factor can
be intrinsic to all polynomial-time algorithms [46].
10

Although both Theorem 4.1 and Theorem 4.2 assume Gaussian distributions, a scrutiny
of the proofs shows that the same results hold if the Gaussian assumption is weakened to
subgaussian. By Theorem 3.2, the rate in Theorem 4.2 is optimal. By Theorem 4.1 and
Theorem 4.2, the choices of the penalty parameters ρ and ρu in (18) and (19) do not depend
on su or sv . Therefore, the proposed estimation scheme (18) – (20) achieves the optimal rate
adaptively over sparsity levels. A full treatment of adaptation to M is beyond the scope of
the current paper, though it seems possible in view of the recent proposals in [6, 35, 12]. A
careful examination of the proofs shows that the dependence of ρ and ρu on M is through
1/2
1/2
kΣx kop kΣy kop and kΣx kop , respectively. When p and m are bounded from above by a
constant multiple of n, we can upper bound the operator norms by the sample counterparts
to remove the dependence of these penalty parameters on M . We conclude this section with
two more remarks.
Remark 4.4. The group sparsity penalty used in the second stage (19) plays an important
role in achieving the optimal rate su (r + log p)/(nλ2 ). Except for the extra λ−2 term, this
convergence rate is a typical one for group Lasso [28]. If we simply use an ℓ1 penalty, then
we will obtain the rate rsu log p/(nλ2 ), which is clearly sub-optimal.
Remark 4.5. Comparing Theorem 3.1 with Theorem 4.2, the adaptive estimation scheme
achieves the optimal rates of convergence for a smaller collection of parameter spaces of
interest due to the more restrictive sample size condition (21). We examine the necessity of
this condition in more details in Section 5 below.

5

Computational Lower Bounds

In this section, we provide evidence that the sample size condition (21) imposed on the
adaptive estimation scheme in Theorems 4.1 and 4.2 is probably unavoidable for any computationally feasible estimator to be consistent. To be specific, we show that for a sequence of
parameter spaces in (5) – (6), if the condition is violated, then any computationally efficient
consistent estimator of sparse canonical coefficients leads to a computationally efficient and
statistically powerful test for the Planted Clique detection problem in a regime where it is
believed to be computationally intractable.
Planted Clique Let N be a positive integer and k ∈ [N ]. We denote by G(N, 1/2) the
Erdős-Rényi graph on N vertices where each edge is drawn independently with probability 1/2, and by G(N, 1/2, k) the random graph generated by first sampling from G(N, 1/2)
and then selecting k vertices uniformly at random and forming a clique of size k on these
vertices. For an adjacency matrix A ∈ {0, 1}N ×N of an instance from either G(N, 1/2) or
G(N, 1/2, k), the Planted Clique detection problem of parameter (N, k) refers to testing the
following hypotheses
H0G : A ∼ G(N, 1/2)

v.s. H1G : A ∼ G(N, 1/2, k).

(22)

It is widely believed that when k = O(N 1/2−δ ), the problem (22) cannot be solved by any
randomized polynomial-time algorithm. In the rest of the paper, we formalize the conjectured
hardness of Planted Clique problem into the following hypothesis.
11

Hypothesis A. For any sequence k = k(N ) such that lim supN →∞
randomized polynomial-time test ψ,

 2
lim inf PH G ψ + PH G (1 − ψ) ≥ .
0
1
N →∞
3

log k
log N

<

1
2

and any

Evidence supporting this hypothesis has been provided in [34, 17]. Computational lower
bounds in several statistical problems have been established by assuming the above hypothesis
and its close variants, including sparse PCA detection [7] and estimation [39] in classes defined
by a restricted covariance concentration condition, submatrix detection [30] and community
detection [21].
Necessity of the sample size condition (21) Under Hypothesis A, the necessity of
condition (21) is supported by the following theorem.
Theorem 5.1. Suppose that Hypothesis A holds and that as n → ∞, p = m satisfying
2n ≤ p ≤ na for some constant a > 1, su = sv , n(log n)5 ≤ cs4u for some sufficiently small
su sv
c > 0, and λ = 7290n(log(12n))
2 . If for some δ ∈ (0, 1),
lim inf
n→∞

(su sv )1−δ log(p + m)
> 0,
nλ2

then for any randomized polynomial-time estimator u
b,
n
1 o 1
> .
lim inf
sup
P L(b
u, u) >
n→∞ P∈P(n,s ,s ,p,m,1,λ;3)
300
4
u v

(23)

(24)

Comparing (21) with (23), we see that subject to a sub-polynomial factor, the condition (21) is necessary to achieve consistent sparse CCA estimation within polynomial time
complexity.
Remark 5.1. The statement in Theorem 5.1 is rigorous only if we assume the computational
complexities of basic arithmetic operations on real numbers and sampling from univariate
continuous distributions with analytic density functions are all Θ(1) [10]. To be rigorous under
the probabilistic Turing machine model [2], we need to introduce appropriate discretization
of the problem and be more careful with the complexity of random number generation.
To convey the key ideas in our computational lower bound construction, we focus on the
continuous case throughout this section and defer the formal discretization arguments to
Section 8 in the supplement [19].
In what follows, we divide the reduction argument leading to Theorem 5.1 into two parts.
In the first part, we show Hypothesis A implies the computational hardness of the sparse
PCA problem under the Gaussian spiked covariance model. In the second part, we show
computational hardness of sparse PCA implies that of sparse CCA as stated in Theorem 5.1.

5.1

Hardness of sparse PCA under Gaussian spiked covariance model

Gaussian single spiked model [24] refers to the distribution Np (0, Σ) where Σ = τ θθ ′ + Ip .
Here, θ is the eigenvector of unit length and τ > 0 is the eigenvalue. Define the following
12

Gaussian single spiked model parameter space for sparse PCA
iid

Q(n, s, p, λ) = {L(W1 , . . . , Wn ) :Wi ∼ Np (0, τ θθ ′ + Ip ),

kθk0 ≤ s, τ ∈ [λ, 3λ]}.

(25)

ep
The minimax estimation rate for θ under the loss kPθb − Pθ k2F is λ+1
nλ2 s log s . See, for instance,
[13]. However, to achieve the above minimax rate via computationally efficient methods such
as those proposed in [9, 29, 13], researchers have required the sample size to satisfy n ≥
2
p
C s λlog
for some sufficiently large constant C > 0. Moreover, no computationally efficient
2
estimator is known to achieve consistency when the sample size condition is violated. As a first
step toward the establishment of Theorem 5.1, we show that Hypothesis A implies hardness
2−δ
p
>0
of sparse PCA under Gaussian spiked covariance model (25) when lim inf n→∞ s nλlog
2
for some δ > 0.
We note that previous computational lower bounds for sparse PCA in [7, 39] cannot be
used here directly because they are only valid for parameter spaces defined via the restricted
covariance concentration (RCC) condition. As pointed out in [39], such parameter spaces
include (but are not limited to) all subgaussian distributions with sparse leading eigenvectors
and the covariance matrices need not be of the spiked form Σ = τ θθ ′ + Ip . Therefore,
the Gaussian single spiked model parameter space defined in (25) only constitutes a small
subset of such RCC parameter spaces. The goal of the present subsection is to establish the
computational lower bound for the Gaussian single spiked model directly.
b 1 , . . . , Wn ) of the leading sparse eigenvector, we
Suppose we have an estimator θb = θ(W
propose the following reduction scheme to transform it into a test for (22). To this end, we
first introduce some additional notation. Consider integers k and N . Define

δN =

k
,
N

ηN =

k
.
45N (log N )2

(26)

For any µ ∈ R, let φµ denote the density function of the N (µ, 1) distribution, and let
φ̄µ =

1
(φµ + φ−µ )
2

(27)

e0
denote the density function of the Gaussian mixture 12 N (µ, 1) +√12 N (−µ,√1). Next, let Φ
be the √
restriction of the N (0, 1) distribution on the interval [−3 log N , 3 log N ]. For any
|µ| ≤ 3 ηN log N , define two probability distributions Fµ,0 and Fµ,1 with densities

−1
[φ̄µ (x) − φ0 (x)] 1{|x|≤3√log N } ,
fµ,0 (x) = M0 φ0 (x) − δN
(28)

−1
(29)
[φ̄µ (x) − φ0 (x)] 1{|x|≤3√log N } ,
fµ,1 (x) = M1 φ0 (x) + δN
R
where the Mi ’s are normalizing constants such that R fµ,i = 1 for i = 0, 1. √It can be
verified that fµ,i are properly defined probability density function when |µ| ≤ 3 ηN log N .
For details, see Lemma 7.4 in the supplement [19].
With the foregoing definition, the proposed reduction scheme can be summarized as
Algorithm 1. Here, the starting point is the adjacency matrix A of the random graph, and
the reduction is well defined for all instances of N ≥ 12n and p ≥ 2n.
13

Algorithm 1: Reduction from Planted Clique to Sparse PCA (in Gaussian Single
Spiked Model)
Input:
1. Graph adjacency matrix A ∈ {0, 1}N ×N ;
2. Estimator θb for the leading eigenvector θ.
Output: A solution to the hypothesis testing problem (22).
e 0 . Set
1 Initialization. Generate i.i.d. random variables ξ1 , . . . , ξ2n ∼ Φ
1/2

µi = ηN ξi ,

2

i = 1, . . . , 2n.

(30)

Gaussianization. Generate two matrices B0 , B1 ∈ R2n×2n where conditioning on the
µi ’s, all the entries are mutually independent satisfying
L((B0 )ij |µi ) = Fµi ,0

and

L((B1 )ij |µi ) = Fµi ,1 .

(31)

Let A0 ∈ {0, 1}2n×2n be the lower–left 2n × 2n submatrix of the matrix A. Generate a
′ ]′ ∈ R2n×p where for each i ∈ [2n], if j ∈ [2n], we set
matrix W = [W1′ , . . . , W2n
Wij = (B0 )ij (1 − (A0 )ij ) + (B1 )ij (A0 )ij .
3

(32)

If 2n < j ≤ p, we let Wij be an independent draw from N (0, 1).
b 1 , . . . , Wn ) be the estimator of the leading
Test Construction. Let θb = θ(W
n
eigenvector by treating {Wi }i=1 as data. It is normalized to be a unit vector. We
reject H0G if
1
θb′ (
n

2n
X

i=n+1

1
Wi Wi′ )θb ≥ 1 + k ηN .
4

(33)

We now explain how the reduction achieves its goal. For simplicity, focus on the case where
p = 2n. Let ǫ = (ǫ1 , ..., ǫ2n ) ∈ {0, 1}2n where ǫi is the indicator of whether the i-th row of A0
(defined in Step 2 of Algorithm 1) belongs to the planted clique or not, and γ = (γ1 , ..., γ2n )
the indicators of the columns of A0 . In what follows, we discuss the distributions of W when
A ∼ H0G and H1G , respectively.
When A ∼ H0G , the ǫi ’s and γj ’s are all zeros. In this case, we can verify that the entries
of W are mutually independent and for each (i, j) the marginal distribution of Wij is close
to the N (0, 1) distribution (c.f., Lemma 7.1 in the supplement [19]). Hence, the rows of W
b 1 , . . . , Wn ) is
are close to i.i.d. random vectors from the Np (0, Ip ) distribution. Since θb = θ(W
2n
2
independent of {Wi }i=n+1 , the LHS of (40) is close in distribution to a χn random variable
scaled by p
n which concentrates around its expected value one. Indeed, it is upper bounded
by 1 + O( log(n)/n) with high probability.
If A ∼ H1G , then the (i, j)-th entry of A0 is an edge in the planted clique if and only if
ǫi = γj = 1. Moreover, the joint distribution of {ǫ1 , . . . , ǫ2n , γ1 , . . . , γ2n } is close to that of
4n i.i.d. Bernoulli random variables {e
ǫ1 , . . . , e
ǫ2n , γ
e1 , . . . , γ
e2n } with success probability δN =
k/N . For simplicity, suppose that these indicators are indeed i.i.d. Bernoulli(δN ) variables
14

{e
ǫ1 , . . . , e
ǫ2n , γ
e1 , . . . , γe2n }. Then, one can show that conditioning on γ
ej = 0, for any i ∈ [2n], the
conditional distribution of (Wij |e
γj = 0), after integrating over the conditional distribution
of e
ǫi , µi and (A0 )ij , is approximately N (0, 1). In contrast, conditioning on γ
ej = 1, for
any i ∈ [2n], the conditional distribution of (Wij |e
γj = 1) is approximately N (0, 1 + ηN ).
Therefore, conditioning on γ
e the distribution of the Wi ’s is close to that of 2n i.i.d. random
vectors sampled from
Np (0, τ θθ ′ + Ip ),

where θ = e
γ /ke
γ k and τ = ηN ke
γ k2 ,

i.e., a Gaussian spiked covariance
P model in (25). Here, the leading eigenvector θ has sparsity
level |supp(θ)| = |supp(e
γ )| = j γ
ej , which concentrates around its mean value nδN ≍ k if
b
N ≍ n. Thus, if θ estimates θ well, then the LHS of (40) approximately
follows a non-central
p
χn distribution scaled by n, which should exceed 1 + O( log(n)/n) with high probability
under the alternative hypothesis. Hence, Algorithm 1 is expected to yield a test with small
error for the Planted Clique problem (22) when θb is a good estimator.
The materialization of the foregoing discussion leads to the following result which demonstrates quantitatively that a decent estimator of the leading sparse eigenvector results in a
good test (by applying the reduction (30) – (33)) for the Planted Clique detection problem
(22).
5

N)
≤ c, cN ≤ n ≤
Theorem 5.2. For some sufficiently small constant c > 0, assume N (log
k4
b
N/12 and p ≥ 2n. Then, for any θ such that
n
1o
≤ β,
(34)
sup
Q kPθb − Pθ k2F >
3
Q∈Q(n,3k/2,p,kηN /2)

the test ψ defined by (30) – (33) satisfies
PH G ψ + PH G (1 − ψ) < β +
0

1

4n
′
+ C(n−1 + N −1 + e−C k ),
N

for sufficiently large n with some constants C, C ′ > 0.
If the estimator θb is uniformly consistent over Q(n, 3k/2, p, kηN /2), then β is close to zero.
Hence the conclusion of Theorem 5.2 implies that for appropriate growing sequences of n, N
and k, the testing error for (22) can be made smaller than any fixed nonzero probability.
Further invoking Hypothesis A, we obtain the following computational lower bounds for
sparse PCA.
Theorem 5.3. Suppose that Hypothesis A holds and that as n → ∞, 2n ≤ p ≤ na for some
s2
constant a > 1, n(log n)5 ≤ cs4 for some sufficiently small c > 0, and λ = 2430n(log(12n))
2 . If
for some δ ∈ (0, 2),
s2−δ log p
lim inf
> 0,
(35)
n→∞
nλ2
b
then for any randomized polynomial-time estimator θ,
lim inf
n→∞

n
1o 1
> .
Q kPθb − Pθ k2F >
3
4
Q∈Q(n,s,p,λ)
sup

15

(36)

In addition to estimation, we can also consider the following sparse PCA detection problem: Let Q denote the joint distribution of W1 , . . . , Wn , and we want to test
H0 : Q ∈ Q(n, s, p, 0),

v.s. H1 : Q ∈ Q(n, s, p, λ).

(37)

iid

Note that Q(n, s, p, 0) contains only one distribution Q0 where Wi ∼ Np (0, Ip ). Given any
testing procedure φ = φ(W1 , . . . , Wn ), we can obtain a solution to (22) by replacing the third
step in Algorithm 1 with the direct testing result of φ(W1 , . . . , Wn ). Following the lines of
the proof of Theorem 5.3, we have the following theorem.
Theorem 5.4. Under the same condition of Theorem 5.3, for any randomized polynomialtime test φ for testing (37),

 1
(38)
sup
Q(1 − φ) ≥ .
lim inf Q0 φ +
n→∞
4
Q∈Q(n,s,p,λ)
Remark 5.2. Theorems 5.3–5.4 are the first computational lower bounds for sparse PCA
that are valid in the setting of Gaussian single spiked covariance models (25).

5.2

Hardness of sparse CCA

In the second step, we show that computational hardness of sparse PCA under Gaussian
spiked covariance model implies the desired result in Theorem 5.1. To this end, we propose
the following reduction.
Algorithm 2: Reduction from Sparse PCA to Sparse CCA
Input:
1. Observations W1 , . . . , Wn ∈ Rp ;
2. Estimator u
b of the first leading canonical correlation coefficient u.
Output: An estimator θb of the leading eigenvector of L(W1 ).
1 Generate i.i.d. random vectors Z1 , . . . , Zn ∼ Np (0, Ip ). Set
1
Xi = √ (Wi + Zi ),
2

2

1
Yi = √ (Wi − Zi ),
2

i = 1, . . . , n.

Compute u
b=u
b(X1 , Y1 , . . . , Xn , Yn ). Set

b 1 , . . . , Wn ) = u
θb = θ(W
b/kb
uk.

(39)

(40)

iid

To see why Algorithm 2 is effective, one can verify that if Wi ∼ Np (0, τ θθ ′ + Ip ), then
iid

(Xi′ , Yi′ )′ ∼ Np+m (0, Σ) where

Σx = Σy =
with u = v = √

θ
,
τ /2+1

λ=

τ /2
τ /2+1 .

τ ′
θθ + Ip , Σxy = Σx (λuv ′ )Σy
2

(41)

This is a special case of the Gaussian canonical pair model

(2). Thus, the leading eigenvector of Wi aligns with the leading canonical coefficient vectors
of (Xi , Yi ). Exploiting this connection, we obtain the following theorem.
16

Theorem 5.5. Consider p = m, su = sv and λ ≤ 1. Then for any u
b such that
n
1 o
≤ β,
P L(b
u, u) >
300
P∈P(n,su ,sv ,p,m,1,λ/3;3)
sup

(42)

the estimator θb defined by Algorithm 2 satisfies

n
1o
≤ β.
Q kPθb − Pθ k2F >
3
Q∈Q(n,s,p,λ)
sup

If we start with an estimator u
b of the leading canonical coefficient vector, then we can
construct the reduction from Planted Clique to sparse CCA directly by essentially following
the steps in Algorithm 1 while using Algorithm 2 to construct θb from u
b in the third step.
Finally, the desired Theorem 5.1 is a direct consequence of Theorems 5.3 and 5.5.

6

Proofs

This section presents proofs of Theorems 4.1 and 4.2. The proofs of the other theoretical
results are given in the supplement [19].

6.1

Proof of Theorem 4.1

Before presenting the proof, we state some technical lemmas. The proofs of all the lemmas
are given in Section 9.3 in the supplement [19]. First, note that the estimator is normalized
b (0)
b (0)
with respect to Σ
x and Σy , while the truth U and V is normalized with respect to Σx and
b (0)
b (0)
Σy . To address this issue, we normalize the truth with respect to Σ
x and Σy to obtain
−1/2 and V
−1/2 . Also define Λ
1/2 Λ(V ′ Σ
1/2 .
e = U (U ′ Σ
b (0)
e = V (V ′ Σ
b (0)
e = (U ′ Σ
b (0)
b (0)
U
x U)
y V)
x U)
y V)
For notational convenience, define
r 
r 
ep 
em 
1
1
su + log
sv + log
, ǫn,v =
.
(43)
ǫn,u =
n
su
n
sv
The following lemma bounds the normalization effect.
Lemma 6.1. Assume ǫ2n,u + ǫ2n,v ≤ c for some sufficiently small constant c ∈ (0, 1). Then
there exist some constants C, C ′ > 0 only depending on c such that
1/2 e
e
e
kΣ1/2
x (U − U )kop ≤ Cǫn,u , kΣy (V − V )kop ≤ Cǫn,v , kΛ − Λkop ≤ C(ǫn,u + ǫn,v ),

with probability at least 1 − exp (−C ′ (su + log(ep/su ))) − exp (−C ′ (sv + log(em/sv ))).

e and Ve , let us state the following lemma, which asserts that the
Using the definitions of U
′
e
e
e
matrix A = U V is feasible to the optimization problem (18).

e=U
e Ve ′ . When A
e exists, we have
Lemma 6.2. Define A
1/2 e b (0) 1/2
b (0)
k(Σ
A(Σy ) k∗ = r
x )

and

17

1/2 e b (0) 1/2
b (0)
k(Σ
A(Σy ) kop = 1.
x )

As was argued in Section 4.1, the set Cr is the convex hull of Or . The following curvature
lemma shows that the relaxation Cr preserves the restricted strong convexity of the objective
function.
Lemma 6.3. Let F ∈ O(p, r), G ∈ O(m, r), K ∈ Rr×r and D = diag(d1 , ..., dr ) with
d1 ≥ ... ≥ dr > 0. If E satisfies kEkop ≤ 1 and kEk∗ ≤ r, then
hF KG′ , F G′ − Ei ≥
Define

dr
kF G′ − Ek2F − kK − DkF kF G′ − EkF .
2
′ b (0)
e xy = Σ
b (0)
Σ
x U ΛV Σy .

(44)

(45)

Lemma 6.4 is instrumental in determining the proper value of the tuning parameter required
in the program (18).
p
Lemma 6.4. Assume r [log(p + m)]/n ≤ c for some sufficiently small constant c ∈ (0, 1).
b (0)
Then there p
exist some constants C, C ′ > 0 only depending on M and c such that ||Σ
xy −
′
−C
e
Σxy ||∞ ≤ C [log(p + m)]/n, with probability at least 1 − (p + m)
.
We also need a lemma on restricted eigenvalue. For any p.s.d. matrix B, define
φB
max (k) =

u′ Bu
,
||u||0 ≤k,u6=0 u′ u

φB
min (k) =

max

u′ Bu
.
||u||0 ≤k,u6=0 u′ u
min

The following lemma is adapted from Lemma 12 in [18], and its proof is omitted.

Lemma 6.5. Assume n1 (ku ∧ p) log(ep/(ku ∧ p)) + (kv ∧ m) log(em/(kv ∧ m)) ≤ c for some
sufficiently small constant c > 0. Then
C ′ > 0 only depending
q there exist some constants C,q
on M and c such that for δu (ku ) =
we have

(ku ∧p) log(ep/(ku ∧p))
n
b (j)

b (j)

b (j)
Σ

b (j)
Σ

and δv (kv ) =

(kv ∧m) log(em/(kv ∧m))
,
n

Σx
x
M −1 − Cδu (ku ) ≤ φΣ
min (ku ) ≤ φmax (ku ) ≤ M + Cδu (ku ),
y
y
(kv ) ≤ φmax
(kv ) ≤ M + Cδv (kv ),
M −1 − Cδv (kv ) ≤ φmin

withprobability at least 1−exp −C ′ (ku ∧p) log(ep/(ku ∧p)) −exp −C ′ (kv ∧m) log(em/(kv ∧
m)) , for j = 0, 1, 2.

Finally, we need a result on subspace distance. Recall that for a matrix F , PF denotes
the projection matrix onto its column subspace.
Lemma 6.6. For any matrix F ∈ O(d, r) and any matrix G ∈ Rd×r , we have
inf kF − GW k2F =
W

1
kPF − PG k2F .
2

If further G ∈ O(d, r), then inf W ∈O(r,r) kF − GW k2F = 21 kPF − PG k2F .
Proofs of Lemma 6.1-6.6 are given in Section 9.3.1 of the supplement [19].

18

b (0)
b (0)
b (0)
b b
Proof of Theorem 4.1. In the rest of this proof, we denote Σ
x , Σy and Σxy by Σx , Σy and
b xy for notational convenience. We also let ∆ = A
b − A.
e The proof consists of two steps.
Σ
b x1/2 ∆Σ
b y1/2 kF . In the second
In the first step, we are going to derive an upper bound for kΣ
b 1/2
b 1/2
step, we derive a generalized cone condition and use it to lower bound kΣ
x ∆Σy kF by a
b 1/2
b 1/2
constant multiple of k∆kF and hence the upper bound on kΣ
x ∆Σy kF leads to an upper
bound on k∆kF .
e and Ve are well-defined with high probability. Thus, A
e is
Step 1.
By Lemma 6.1, U
well-defined with high probability, and we have
′
1/2
e
kΣ1/2
x (A − U V )Σy kop ≤ C(ǫn,u + ǫn,v ).

(46)

with probability at least 1 − exp (−C ′ (su + log(ep/su ))) − exp (−C ′ (sv + log(em/sv ))). Ace is feasible. Then, by the definition of A,
b we have
cording to Lemma 6.2, A
b xy , Ai
b − ρ||A||
b 1 ≥ hΣ
b xy , Ai
e − ρ||A||
e 1.
hΣ

After rearrangement, we have


b xy − Σ
e xy , ∆i,
e xy , ∆i ≤ ρ ||A||
e 1 − ||A
e + ∆||1 + hΣ
− hΣ

(47)

e xy is defined in (45). For the first term on the right hand side of (47), we have
where Σ
e 1 − ||A
e + ∆||1 = ||A
eSu Sv ||1 − ||A
eSu Sv + ∆Su Sv ||1 − ||∆(S S )c ||1
||A||
u v
≤ ||∆Su Sv ||1 − ||∆(Su Sv )c ||1 .

b xy − Σ
e xy , ∆i ≤ ||Σ
b xy −
For the second term on the right hand side of (47), we have hΣ
e xy ||∞ ||∆||1 . Thus when
Σ
b xy − Σ
e xy ||∞ ,
ρ ≥ 2||Σ
(48)

we have

e xy , ∆i ≤ 3ρ ||∆Su Sv ||1 − ρ ||∆(S S )c ||1 .
− hΣ
u v
2
2
Using Lemma 6.3, we can lower bound the left hand side of (49) as
e xy , ∆i = hΣ
b x1/2 U ΛV ′ Σ
b 1/2
b 1/2 e b b 1/2
− hΣ
y , Σx (A − A)Σy i
b x1/2 U
eΛ
e Ve ′ Σ
b 1/2
b 1/2 e b b 1/2
= hΣ
y , Σx (A − A)Σy i
≥

1
b 1/2
e b b 1/2 2
b 1/2 e b b 1/2
λr kΣ
x (A − A)Σy kF − δkΣx (A − A)Σy kF ,
2

e − ΛkF . Combining (49) and (50), we have
where δ = kΛ

(49)

(50)

b 1/2 ∆Σ
b 1/2 k2 ≤ 3ρ||∆Su Sv ||1 − ρ||∆(S S )c ||1 + 2δkΣ
b 1/2 ∆Σ
b 1/2 kF
λr kΣ
x
y
F
x
y
u v
1/2 b 1/2
b
≤ 3ρ||∆Su Sv ||1 + 2δkΣx ∆Σy kF .

(51)

2
2
b 1/2
b 1/2 2
kΣ
x ∆Σy kF ≤ 6ρ||∆Su Sv ||1 /λr + 4δ /λr .

(53)

(52)

Solving the quadratic equation (52) by Lemma 2 of [13], we have

19

Combining (51) and (53), we have
b 1/2
b 1/2 2
0 ≤ 3ρ||∆Su Sv ||1 − ρ||∆(Su Sv )c ||1 + δ2 /λr + λr kΣ
x ∆Σy kF
≤ 9ρ||∆Su Sv ||1 − ρ||∆(Su Sv )c ||1 + 5δ2 /λr ,

(54)

which gives rise to the generalized cone condition that we are going to use in Step 2. Finally,
√
by the bound ||∆Su Sv ||1 ≤ su sv ρk∆Su Sv kF and (53), we have
√
2
2
b 1/2
b 1/2 2
kΣ
x ∆Σy kF ≤ 6 su sv ρk∆Su Sv kF /λr + 4δ /λr ,

(55)

which completes the first step.
Step 2. By (54), we have obtained the following condition

||∆(Su Sv )c ||1 ≤ 9||∆Su Sv ||1 + 5δ2 /(ρλr ).

(56)

Due to the existence of the extra term 5δ2 /(ρλr ) on the RHS, we call it a generalized cone
b x1/2 ∆Σ
b 1/2
condition. In this step, we are going to lower bound kΣ
y kF by k∆kF on the generalized
cone. Motivated by the argument in [8], let the index set J1 = {(ik , jk )}tk=1 in (Su × Sv )c
correspond to the entries with the largest absolute values in ∆, and we define the set Je =
(Su × Sv ) ∪ J1 . Now we partition Jec into disjoint subsets J2 , ..., JK of size t (with |JK | ≤ t),
such that Jk is the set of
(double) indices corresponding to the entries of t largest absolute
Sk−1
values in ∆ outside Je ∪ j=2
Jj . By triangle inequality,
b 1/2
b 1/2
b 1/2
b 1/2
kΣ
x ∆Σy kF ≥ kΣx ∆JeΣy kF −
≥

q

bx
φΣ
min (su

+

by
Σ
(sv
t)φmin

K
X
k=2

+ t)k∆JekF −

b 1/2
b 1/2
kΣ
x ∆Jk Σy kF

q

K
X
by
bx
Σ
φΣ
(t)φ
(t)
k∆Jk kF .
max
max
k=2

By the construction of Jk , we have
K
X
k=2

k∆Jk kF ≤

K
K
X
√ X
||∆Jk−1 ||1 ≤ t−1/2 ||∆(Su Sv )c ||1
||∆Jk ||∞ ≤ t−1/2
t
k=2

k=2

r


su sv
5δ2
5δ2
−1/2
√ ,
≤t
9||∆Su Sv ||1 +
k∆JekF +
≤9
ρλr
t
ρλr t

(57)

where we have used the generalized cone condition (56). Hence, we have the lower bound
√
2
b 1/2
b 1/2
kΣ
x ∆Σy kF ≥ κ1 k∆JekF − κ2 δ /(ρλr t),
with

κ1 =

q

bx
φΣ
min (su

+

by
Σ
(sv
t)φmin

q
by
bx
Σ
κ2 = 5 φΣ
max (t)φmax (t).

r
q
by
su sv
bx
Σ
+ t) − 9
φΣ
max (t)φmax (t),
t

20

(58)

Taking t = c1 su sv for some sufficiently large constant c1 > 1, with high probability, κ1 can
be lower bounded by a positive constant κ0 only depending
p that by
p on M . To see this, note
−1
Lemma 6.5, (58) can be lower bounded by the difference of M − Cδu (2c1 su sv ) M −1 − Cδv (2c1 su sv )
p
−1/2 p
and 9c1
M + Cδu (c1 su sv ) M + Cδv (c1 su sv ), where δu and δv are defined as in Lemma
6.5. It is sufficient to show that δu (2c1 su sv ), δv (2c1 su sv ), δu (c1 su sv ) and δv (c1 su sv ) are sufficiently small to get a positive absolute constant κ0 . For the first term, when 2c1 su sv ≤ p,
it is bounded by 2c1 su svnlog(ep) and is sufficiently small under the assumption (13). When
2c1 su sv > p, it is bounded by 2c1 nsu sv and is also sufficiently small under (13). The same argument also holds for the other terms. Similarly, κ2 can be upper bounded by some constant.
Together with (55), this brings the inequality
√

√
k∆Jek2F ≤ C1 ( su sv ρ/λr )k∆JekF + C2 δ2 /λ2r + (δ2 /(ρλr t))2 .
Solving this quadratic equation, we have
k∆Jek2F ≤ C

 s s ρ2
δ2  δ2 2 
u v
√
.
+
+
λ2r
λ2r
ρλr t

(59)

By (57), we have
r
5δ2
su sv
√ .
k∆JekF +
k∆Jk kF ≤ 9
k∆Jec kF ≤
t
ρλ
t
r
k=2
K
X

(60)

Summing (59)
p and (60), we obtain a bound for k∆kF . According to Lemma 6.4, we may
choose ρ = γ [log(p
p + m)]/n for some large γ, so that
√ (48) holds with high probability. By
′
Lemma 6.1, δ ≤ C r(su + sv + log(p + m))/n ≤ C ρ t with high probability. Hence,
√
k∆kF ≤ C su sv ρ/λr ,

(61)

with high probability. This completes the second step. Finally, the triangle inequality leads
b − U V ′ kF ≤ k∆kF + kA
e − U V ′ kF . By (46) and (61), the proof is complete.
to kA

6.2

Proof of Theorem 4.2

b (1) − U ∗ .
Define U ∗ = U ΛV ′ Σy Vb (0) and ∆ = U

r+log p
≤ c for some sufficiently small constant c ∈
n
′
C, C > 0 only depending on M and c such that max

Lemma 6.7. Assume

(0, 1). Then there
b (1) b (0) −
exist some constants
1≤j≤p ||[Σxy V
p

∗
′
b (1)
Σ
x U ]j· || ≤ C (r + log p)/n, with probability at least 1 − exp − C (r + log p) .
The proof of Lemma 6.7 is given in Section 9.3.1 of the supplement [19].

b (1)
b (1)
b (1)
b b
Proof of Theorem 4.2. In the rest of this proof, we denote Σ
x , Σy and Σxy by Σx , Σy
b xy for simplicity of notation. Note that they depends on D1 , while the estimator Vb (0)
and Σ
depends on D0 . Hence, Vb (0) is independent of the sample covariance matrices occurring
in this proof. The proof consists of three steps. In the first step, we derive a bound for
b x ∆). In the second step, we derive a cone condition and use it to obtain a bound for
Tr(∆′ Σ
21

b x ∆) upper bounds k∆kF . In the last step, we derive the desired
k∆kF by arguing that Tr(∆′ Σ
b , U ).
bound for L(U
b (1) , we have Tr((U
b (1) )′ Σ
b xU
b (1) )−2 Tr((U
b (1) )′ Σ
b xy Vb (0) )+ρu Pp ||U
b (1) || ≤
Step 1. By definition of U
j=1
j·
b x U ∗ ) − 2 Tr((U ∗ )′ Σ
b xy Vb (0) ) + ρu Pp ||U ∗ ||. After rearrangement, we have
Tr((U ∗ )′ Σ
j·
j=1
′b

Tr(∆ Σx ∆) ≤ ρu

p h
X
j=1

i
h
i
b xy Vb (0) − Σ
b xU ∗ ) .
||Uj·∗ || − ||Uj·∗ + ∆j· || + 2 Tr ∆′ (Σ

(62)

For the first term on the right hand side of (62), we have
p
X
j=1

X
X
X

||Uj·∗ + ∆j· || −
||∆j· ||
||Uj·∗ || −
||Uj·∗ || − ||Uj·∗ + ∆j· || =
≤

X

j∈Su

j∈Suc

j∈Su

j∈Su

||∆j· || −

X

j∈Suc

||∆j· ||.

For the second term on the right hand side of (62), we have
p


 X
′ b
(0)
∗
b xy Vb (0) − Σ
b x U ∗ ]j· ||,
b
b
||∆j· || max ||[Σ
Tr ∆ (Σxy V − Σx U ) ≤
1≤j≤p

j=1

where [·]j· means the j-th row of the corresponding matrix. When
b xy Vb (0) − Σ
b x U ∗ ]j· ||,
ρu ≥ 4 max ||[Σ
1≤j≤p

we have

Since

P

j∈Su ||∆j· || ≤

b x ∆) ≤
Tr(∆′ Σ

ρu X
3ρu X
||∆j· || −
||∆j· ||.
2
2
c
j∈Su

(63)

(64)

j∈Su

√ qP
2
su
j∈Su ||∆j· || , (64) can be upper bounded by
√
s
3 s u ρu X
Tr(∆ Σx ∆) ≤
||∆j· ||2 .
2
′b

(65)

j∈Su

This completes the first step.
Step 2. The inequality (64) implies the cone condition
X
X
||∆j· ||.
||∆j· || ≤ 3
j∈Suc

(66)

j∈Su

Let the index set J1 = {j1 , ..., jt } in Suc correspond to the rows with the largest ℓ2 norm in ∆,
and we define the extended support Seu = Su ∪ J1 . Now we partition Seuc into disjoint subsets
J2 , ..., JK of size t (with |JK | ≤ t), such that Jk is the set of indices corresponding to the t

22

Sk−1
b x ∆) = kn−1/2 X∆k2 ,
rows with largest ℓ2 norm in ∆ outside Seu ∪ j=2
Jj . Note that Tr(∆′ Σ
F
where X = [X1 , ..., Xn ]′ ∈ Rn×p denotes the data matrix. By triangle inequality, we have
X
kn−1/2 X∆kF ≥ kn−1/2 X∆Seu ∗ kF −
kn−1/2 X∆Jk ∗ kF
≥

q

k≥2

b

x
φΣ
min (su + t)k∆Seu ∗ kF −

q

b

x
φΣ
max (t)

X
k≥2

k∆Jk ∗ kF ,

where for a subset B ⊂ [p], ∆B∗ = (∆ij 1{i∈B,j∈[r]} ), and
X
k≥2

√ X
√ X1 X
t
max ||∆j· || ≤ t
||∆j· ||
j∈Jk
t
k≥2
k≥2 j∈Jk−1
X
X
||∆j· ||
≤ t−1/2
||∆j· || ≤ 3t−1/2

k∆Jk ∗ kF ≤

j∈Suc

(67)

j∈Su

r s
r
su X
su
2
≤ 3
||∆j· || ≤ 3
k∆Seu ∗ kF .
t
t

(68)

j∈Su

In the above derivation, we have used the construction
of Jk and the
qcone condition (66).
q
p
bx
b
Σ
x
Hence, kn−1/2 X∆kF ≥ κk∆Seu ∗ kF with κ = φmin
(su + t) − 3 stu φΣ
max (t). In view of
Lemma 6.5, taking t = c1 su for some sufficiently large constant c1 , with high probability, κ
can be lower bounded by a positive constant κ0 only depending on M . Combining with (65),
we have
√
k∆Seu ∗ kF ≤ C su ρu /(2κ20 ).
(69)
By (67)-(68), we have

k∆(Seu )c ∗ kF ≤

X
p
−1/2
k∆Jk ∗ kF ≤ 3 su /tk∆Seu ∗ kF ≤ 3c1 k∆Seu ∗ kF .

(70)

k≥2

√
and
(70),
we
have
k∆k
≤
C
su ρ. By Lemma 6.7, we may choose ρu ≥
Summing
(69)
F
q
p
for some large γu so that (63) holds with high probability. Hence, k∆kF ≤
γu r+log
p n
C su (r + log p)/n with high probability. This completes the second step.
Step 3. Using the same argument in Step 2 of the proof of Theorem 3.1 (see supplementary
b , U ). The proof is complete.
material [19]), we obtain the desired bound for L(U

References

[1] T. W. Anderson. An Introduction to Multivariate Statistical Analysis. Wiley, 1958.
[2] S. Arora and B. Barak. Computational complexity: a modern approach. Cambridge
University Press, 2009.
[3] B. B. Avants, P. A. Cook, L. Ungar, J. C. Gee, and M. Grossman. Dementia induces
correlated reductions in white matter integrity and cortical thickness: a multivariate
neuroimaging study with sparse canonical correlation analysis. NeuroImage, 50(3):1004–
1016, 2010.
23

[4] Z. Bao, J. Hu, G. Pan, and W. Zhou. Canonical correlation coefficients of highdimensional normal vectors: finite rank case. arXiv preprint arXiv:1407.7194, 2014.
[5] S. R. Becker, E. J. Candès, and M. C. Grant. Templates for convex cone problems with
applications to sparse signal recovery. Mathematical Programming Computation, 3(3):
165–218, 2011.
[6] A. Belloni, V. Chernozhukov, and L. Wang. Square-root lasso: pivotal recovery of sparse
signals via conic programming. Biometrika, 98(4):791–806, 2011.
[7] Q. Berthet and P. Rigollet. Complexity theoretic lower bounds for sparse principal
component detection. In Conference on Learning Theory, pages 1046–1066, 2013.
[8] P. J. Bickel, Y. Ritov, and A. B. Tsybakov. Simultaneous analysis of Lasso and Dantzig
selector. The Annals of Statistics, 37(4):1705–1732, 2009.
[9] A. Birnbaum, I.M. Johnstone, B. Nadler, and D. Paul. Minimax bounds for sparse PCA
with noisy high-dimensional data. The Annals of Statistics, 41(3):1055–1084, 2013.
[10] L. Blum, F. Cucker, M. Shub, and S. Smale. Complexity and Real Computation. Springer,
2012.
[11] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and
statistical learning via the alternating direction method of multipliers. Foundations and
Trends in Machine Learning, 3(1):1–122, 2011.
[12] F. Bunea, J. Lederer, and Y. She. The group square-root lasso: Theoretical properties
and fast algorithms. Information Theory, IEEE Transactions on, 60(2):1313–1325, 2014.
[13] T. T. Cai, Z. Ma, and Y. Wu. Sparse PCA: Optimal rates and adaptive estimation. The
Annals of Statistics, 41(6):3074–3110, 2013.
[14] M. Chen, C. Gao, Z. Ren, and H. H. Zhou. Sparse CCA via precision adjusted iterative
thresholding. arXiv preprint arXiv:1311.6186, 2013.
[15] P. Diaconis and D. Freedman. Finite exchangeable sequences. The Annals of Probability,
pages 745–764, 1980.
[16] J. Douglas and H. Rachford. On the numerical solution of heat conduction problems in
two and three space variables. Transactions of the AMS, pages 421–439, 1956.
[17] V. Feldman, E. Grigorescu, L. Reyzin, S. Vempala, and Y. Xiao. Statistical algorithms
and a lower bound for detecting planted cliques. In Proceedings of the forty-fifth annual
ACM symposium on Theory of computing, pages 655–664. ACM, 2013.
[18] C. Gao, Z. Ma, Z. Ren, and H. H. Zhou. Minimax estimation in sparse canonical
correlation analysis. The Annals of Statistics, 43(5):2168–2197, 2015.
[19] C. Gao, Z. Ma, and H. H. Zhou. Supplement to “Sparse CCA: Adaptive estimation and
computational barriers”. 2015.

24

[20] G. H. Golub and C. F. Van Loan. Matrix Computations (3rd ed.). The Johns Hopkins
University Press, 1996.
[21] B. Hajek, Y. Wu, and J. Xu. Computational lower bounds for community detection on
random graphs. arXiv preprint arXiv:1406.6625, 2014.
[22] D. R. Hardoon and J. Shawe-Taylor. Sparse canonical correlation analysis. Machine
Learning, 83(3):331–353, 2011.
[23] H. Hotelling. Relations between two sets of variates. Biometrika, 28:321–377, 1936.
[24] I. M. Johnstone and A. Y. Lu. On consistency and sparsity for principal components
analysis in high dimensions. Journal of the American Statistical Association, 104(486):
682–693, 2009.
[25] I.M. Johnstone. Multivariate analysis and Jacobi ensembles: Largest eigenvalue, TracyWidom limits and rates of convergence. Ann. Stat., 36:2638–2716, 2008.
[26] L. Le Cam. Asymptotic Methods in Statistical Theory. Springer-Verlag, 1986.
[27] K.-A. Lê Cao, P. G. P. Martin, C. Robert-Granié, and P. Besse. Sparse canonical
methods for biological data integration: application to a cross-platform study. BMC
Bioinformatics, 10(1):34, 2009.
[28] K. Lounici, M. Pontil, S. van de Geer, and A. B. Tsybakov. Oracle inequalities and
optimal inference under group sparsity. The Annals of Statistics, 39(4):2164–2204, 2011.
[29] Z. Ma. Sparse principal component analysis and iterative thresholding. The Annals of
Statistics, 41(2):772–801, 2013.
[30] Z. Ma and Y. Wu. Computational barriers in minimax submatrix detection. arXiv
preprint arXiv:1309.5914, 2013.
[31] K. Mardia, J. Kent, and J. Bibby. Multivariate Analysis. Academic Press, 1979.
[32] Sahand N Negahban, Pradeep Ravikumar, Martin J Wainwright, and Bin Yu. A unified
framework for high-dimensional analysis of m-estimators with decomposable regularizers.
Statistical Science, 27(4):538–557, 2012.
[33] E. Parkhomenko, D. Tritchler, and J. Beyene. Sparse canonical correlation analysis
with application to genomic data integration. Statistical Applications in Genetics and
Molecular Biology, 8(1):1–34, 2009.
[34] B. Rossman. Average-case complexity of detecting cliques. PhD thesis, Massachusetts
Institute of Technology, 2010.
[35] T. Sun and C.-H. Zhang. Scaled sparse linear regression. Biometrika, 99(4):879–898,
2012.
[36] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. arXiv
preprint arXiv:1011.3027, 2010.
25

[37] V. Q. Vu, J. Cho, J. Lei, and K. Rohe. Fantope projection and selection: A near-optimal
convex relaxation of sparse pca. In Advances in Neural Information Processing Systems,
pages 2670–2678, 2013.
[38] S. Waaijenborg and A. H. Zwinderman. Sparse canonical correlation analysis for identifying, connecting and completing gene-expression networks. BMC Bioinformatics, 10
(1):315, 2009.
[39] T. Wang, Q. Berthet, and R. J. Samworth. Statistical and computational trade-offs in
estimation of sparse principal components. arXiv preprint arXiv:1408.5369, 2014.
[40] G. A. Watson. On matrix approximation problems with Ky Fan k norms. Numerical
Algorithms, 5(5):263–272, 1993.
[41] P.-Å. Wedin. Perturbation bounds in connection with singular value decomposition. BIT
Numerical Mathematics, 12(1):99–111, 1972.
[42] A. Wiesel, M. Kliger, and A. O. Hero III. A greedy approach to sparse canonical
correlation analysis. arXiv preprint arXiv:0801.2748, 2008.
[43] D. M. Witten, R. Tibshirani, and T. Hastie. A penalized matrix decomposition, with
applications to sparse principal components and canonical correlation analysis. Biostatistics, 10(3):515–534, 2009.
[44] B. Yu. Assouad, Fano, and Le Cam. In Festschrift for Lucien Le Cam, pages 423–435.
Springer, 1997.
[45] M. Yuan and Y. Lin. Model selection and estimation in regression with grouped variables.
Journal of the Royal Statistical Society: Series B, 68(1):49–67, 2006.
[46] Y. Zhang, M. J. Wainwright, and M. I. Jordan. Lower bounds on the performance of
polynomial-time algorithms for sparse linear regression. arXiv preprint arXiv:1402.1918,
2014.

26

Supplement to “Sparse CCA: Adaptive Estimation and Computational
Barriers”
7

Proofs of Results in Section 5

In this section, we present the proofs of Theorems 5.1–5.5. Here we do not consider the
issue of discretization. The main purpose is to help the readers get the intuition behind the
problem without worrying about rigor at the theoretical computer science level. A rigorous
treatment of the computational lower bounds is deferred to Section 8 where the asymptotic
equivalent discretization and the statement of rigorous results for the discretized models will
be presented.

7.1

Proof of Theorem 5.2

The reason for introducing the two distributions (28)–(29) is to match specific mixtures of
them to φ0 and φ̄µ respectively as summarized in the following lemma. For two probability
distributions P and Q, the total variation distance is defined as TV(P, Q) = supB |P(B) −
Q(B)|. We also write TV(p, q) if p and q are the densities of P and Q, respectively.
Lemma 7.1. There exists
√ an absolute constant C > 0, such that for all integers N ≥ 12,
k ≤ N/12 and all |µ| ≤ 3 ηN log N ,
TV(hµ,0 , φ0 ) ≤ CN −3

and

TV(hµ,1 , φ̄µ ) ≤ CN −3 ,

where hµ,0 = 12 (fµ,0 + fµ,1 ) and hµ,1 = δN fµ,1 + (1 − δN ) 21 (fµ,0 + fµ,1 ).
The proof of the above lemma will be given in Section 7.4. To facilitate the proof of Theorem 5.2, we need to state and prove another two lemmas which characterize the distributions
of the Wi ’s under H0G and H1G respectively. Let L({Wi }2n
i=1 ) denote the joint distribution of
2n
′
{Wi }i=1 . In addition, denote Np (0, τ θθ + Ip ) by Qθ,τ . When τ = 0, Np (0, Ip ) is denoted by
Q0 . The first lemma shows that under H0G , the joint distribution of {Wi }2n
i=1 is close in total
variation to that of a random sample of size 2n from Q0 .
Lemma 7.2. Suppose A ∼ G(N, 1/2). There exists an absolute constant C > 0 such that
2n
−1
TV(L({Wi }2n
.
i=1 ), Q0 ) ≤ CN

Proof. Recall ηN defined in (26) and hµ,0 in Lemma√7.1. Let ν be
√ N (0, ηN ), and ν̄ be the
distribution obtained by restricting ν on the set [−3 ηN log N , 3 ηN log N ]. Then the µi ’s
in (30) are i.i.d. r.v.’s following the distribution ν̄.
For each i ∈ [2n] and each j ∈ [2n], define i.i.d. random variables W ij ∼ N (0, 1). For each
i ∈ [2n] and 2n < j ≤ p, define W ij = Wij . Let W i = (W i1 , ..., W ip )′ . It is straightforward
2n
to verify that L({W i }2n
i=1 ) = Q0 . By the data-processing inequality, we have
2n
n
n
TV(L({Wi }i=1
), Q2n
0 ) ≤ TV(L({Wi }i=1 ), L({W i }i=1 )).

27

Hence, it is sufficient to bound TV(L({Wi }ni=1 ), L({W i }ni=1 )). Conditioning on µi , Wij follows
hµi ,0 when A ∼ G(N, k). Therefore,
Z
TV(Wij , W ij ) = TV( hµi ,0 dν̄(µi ), φ0 )
≤

sup
√

|µi |≤3 ηN log N

TV(hµi ,0 , φ0 ) ≤ CN −3 .

Here the
last inequality is due to Lemma 7.1. Applying Lemma 7 of [30], we obtain TV(L({Wi }ni=1 ), L({W i }ni=1 ))
P
n Pn
−1 . This completes the proof.
i=1
j=1 TV(Wij , W ij ) ≤ CN
The next lemma shows that the joint distribution {Wi }2n
i=1 is close in total variation to a
mixture of the joint distribution of a random sample of size 2n from Qθ,τ . Here, the mixture
is defined by a prior distribution π on the (θ, τ ) pair, which is supported on a region where θ
is sparse and τ is bounded away from zero. For notational convenience, for any distribution
R
Pβ indexed by parameter β ∈ B and any probability
R measure π on B, we let Pβ dν(β)
denote the probability measure P defined by P(E) = Pβ (E)dν(β) for any event E. When
β ∼ ν Ris a random variable and Pβ = L(W |β) is the conditional distribution of W |β, we also
write L(W |β)dν(β) to represent the marginal distribution of W after integrating out β.
Lemma 7.3. Suppose A ∼ G(N, 1/2, k). There exists a distribution π supported on the set

(θ, τ ) : θ ∈ S p−1 , |supp(θ)| ≤ 3k/2, τ ∈ [kηN /2, 3kηN /2] ,
(71)
such that for some absolute constants C1 , C2 > 0,


Z
4n
1
−C2 k
2n
2n
.
+
+
TV(L({Wi }i=1 ), Qθ,τ dπ(θ, τ )) ≤ C1 e
N
N

Proof. Recall ηN defined in (26) and hµ,0 and hµ,1 defined in Lemma 7.1. As in the proof of
Lemma
√ 7.2, let ν√ be N (0, ηN ), and ν̄ the distribution obtained by restricting ν on the set
[−3 ηN logRN , 3 ηN log N ]. Then the µi ’s in (30) are i.i.d. r.v.’s followingRν̄. Simple calculus
shows that φ0 (x)dν(µ) = φ0 (x) is the density function of N (0, 1), and φ̄µ (x)dν(µ) gives
the density function of N (0, 1 + ηN ).
We first focus on the case p = 2n. The case of p ≥ 2n will be treated at the end of the
proof. Recall that (ǫ1 , ..., ǫ2n ) are the indicators of the rows of A0 whether the corresponding
vertices belong to the planted clique, and (γ1 , ..., γp ) are the corresponding indicators of the
columns of A0 . Let (e
ǫ1 , ..., e
ǫ2n ) and (e
γ1 , ..., γep ) be i.i.d. Bernoulli random variables with
e
e0 )ij = 1 if e
mean δN = k/N . Define a matrix A0 , where an entry (A
ǫi = γ
ej = 1 and is an
f
independent instantiation of the Bernoulli(1/2) distribution otherwise. Then, we define W
with entries
fij = (B0 )ij (1 − (A
e0 )ij ) + (B1 )ij (A
e0 )ij .
W
Then, by Theorem 4 of [15] and the data-processing inequality, we have
f ), L(W )) ≤ TV(L(e
TV(L(W
ǫ, γ
e), L(ǫ, γ)) ≤

4n
.
N

(72)

f , conditioning on µi and
Recall hµ,0 and hµ,1 defined in Lemma 7.1. By the definition of W
fij ∼ hµ ,0 , while conditioning on µi and γ
fij ∼ hµ ,1 .
γ
ej = 0, W
ej = 1, W
i
i
28

Further define W ij by setting
W ij |(e
γj = 0, µi ) ∼ φ0 ,

W ij |(e
γj = 1, µi ) ∼ φ̄µi ,

where φ̄µi is √
defined according to (27). By Lemma 7.1 and Lemma 7 of [30], uniformly over
maxi |µi | ≤ 3 ηN log N, we have
p
2n X


 X

f |e
fij |e
TV L(W
γ , µ), L(W |e
TV L(W
γj , µi ), L(W ij |e
γ , µ) ≤
γj , µi ) ≤ CN −1
i=1 j=1

for some constant C > 0.
Next, we integrate the above bound over µ. To this end, first note that
Z
Z
dν(µ) =
φ0 (x)dx ≤ CN −4 .
TV(ν, ν̄) =
√
√
|x|>3 log N

|µ|>3 ηN log N

R

R
f |e
f |e
With slight abuse of notation, let L(W
γ , µ)dν̄(µ) (resp.
L(W
γ , µ)dν(µ)) denote the
f |e
conditional distribution
of
W
γ
if
the
coordinates
of
µ
=
(µ
,
.
.
.
,
µ
1
2n ) were i.i.d. following
R
R
ν̄ (resp. ν), and let L(W |e
γ , µ)dν̄(µ) and L(W |e
γ , µ)dν(µ) be analogously defined. Then,
conditioning on γ
e, we obtain
Z
Z
f
TV( L(W |e
γ , µ)dν̄(µ), L(W |e
γ , µ)dν(µ))
Z
Z
f |e
≤ TV( L(W
γ , µ)dν̄(µ), L(W |e
γ , µ)dν̄(µ))
Z
Z
+ TV( L(W |e
γ , µ)dν̄(µ), L(W |e
γ , µ)dν(µ))
≤

sup
√

maxi |µi |≤3 ηN log N

f |e
TV(L(W
γ , µ), L(W |e
γ , µ)) + CnTV(ν̄, ν) ≤ CN −1 .

Here, the first inequality comes from the triangle inequality, the second
Pnfrom thePdefinition
of total variation distance. For each given γ
e = (e
γ1 , ..., γen ), define s = j=1 γ
ej = nj=1 γ
ej2 =
ke
γj k2 , θ = s−1/2 γ
e and τ = sηN . Note that both θ and τ are functions of γ
e. Then observe
R
R
f |e
f |e
,
which
implies
for
L(
W
γ
)
=
L(
W
γ
,
µ)dν̄(µ),
γ , µ)dν(µ) = Q2n
that L(W |e
θ,τ


−1
f |e
.
TV L(W
γ ), Q2n
θ,τ ≤ CN

Define the event Q = {e
γ : |s − k| ≤ k/2}. Then, by Bernstein’s inequality, P(Qc ) ≤ e−Ck . Let
π
e be the joint distribution of (θ, τ ), and π be the distribution obtained from renormalizing
the restriction of π
e on {(θ(e
γ ), τ (e
γ )) : γ
e ∈ Q} which is exactly the set in (71). Then we have
c
−Ck
f |e
f |θ, τ ) since there
TV(π, π
e) ≤ CP(Q ) ≤ Ce
. In addition, we note that L(W
γ ) = L(W
exists one-to-one identification between the pair (θ, τ ) and e
γ . Therefore, we have
Z
Z
f ), L(W
f |θ, τ )dπ(θ, τ ))
f ), Q2n dπ(θ, τ )) ≤ TV(L(W
TV(L(W
θ,τ
Z
Z
f |θ, τ )dπ(θ, τ ), Q2n dπ(θ, τ ))
+TV( L(W
θ,τ
f |θ, τ ), Q2n )
≤ TV(e
π , π) + sup TV(L(W
θ,τ


θ,τ

−Ck

≤ C e

29

+N

−1



.

R
f ) = L(W
f |θ, τ )de
Here, the Rsecond inequality holds since L(W
π(θ, τ ). Hence, by (72),

4n
−Ck
−1
2n
+N
+ N . Note that on the support of π, the
TV(L(W ), Qθ,τ dπ(θ, τ )) ≤ C e
parameter (θ, τ ) belongs to the set (71). An application of data-processing inequality leads
to the conclusion. When p ≥ 2n, we may first analyze the distribution of the first 2n coordinates using the above arguments. The remaining 2n − p coordinates are exact, and the total
variation bound is zero.
P
′
b
Proof of Theorem 5.2. Abbreviate n1 2n
i=n+1 Wi Wi by Σ. We can rewrite the testing function
ψ as
n
o
b θb ≥ 1 + kηN /4 .
ψ(W ) = ψ(A, µ, B0 , B1 ) = 1 θb′ Σ

Here, µ = (µ1 , . . . , µ2n ) collects the random variables in (30). Thus, it is clear that ψ is a
randomized test for the Planted Clique detection problem (22). Note that for any (θ, τ ) in
the support of π, we have
 3k
kηN 
.
(73)
Qnθ,τ ∈ Q n, , p,
2
2
We now bound the testing errors. For Type-I error, Lemma 7.2 implies
PH G ψ ≤ Qn0 ψ + CN −1 .
0

b are independent. Conditioning on θb and using Bernstein’s
Note that under Qn0 , θb and Σ
inequality, we have
2n 

X
b 2 > 1 + kηN ,
b θb = 1 + 1
|θb′ Wi |2 − ||θ||
θb′ Σ
n
4
i=n+1



4
b
with probability at most exp − N 2Cnk
(log N )4 . Integrating over θ, we have

PH G ψ ≤ exp −
0

Cnk 4 
+ CN −1 ≤ C(n−1 + N −1 ),
N 2 (log N )4

where the last inequality holds under the assumptions
some sufficiently small constant c > 0.
Turn to the Type-II error. Lemma 7.3 implies

N (log N )5
k4

(74)

≤ c and cN ≤ n ≤ N/12 for

 4n

,
(75)
PH G (1 − ψ) ≤ Qπ (1 − ψ) + C e−Ck + N −1 +
1
N
R
where we have used the notation Qπ = Qnθ,τ dπ. For each Qnθ,τ in the support of π, Wn+i has
√
representation Wn+i = τ gi θ + ǫi , where the gi ’s and the ǫi ’s are independently distributed
according to N (0, 1) and Np (0, Ip ), and are independent across i = 1, ..., n, and τ ≥ kηN /2.
Therefore,
√
n
n
n
 1X
1 X
2
τ b′ X ′ b
2
′
2
′
′
2
b θb = τ |θb θ|
gi ǫi θ.
gi +
|θb ǫi | +
θθ
θb Σ
n
n
n
i=1

i=1

30

i=1

After rearrangement, we have
n

b θb − (1 + τ )
θb′ Σ

≤

1X 2
(gi − 1) + τ min{|(θb − θ)′ θ|2 , |(θb + θ)′ θ|2 }
n

+

i=1
n
X

1
n

i=1

n

(|θb′ ǫi |2 − 1) +

2X
b ,
gi (ǫ′i θ)
n
i=1


where min{|(θb − θ)′ θ|2 , |(θb + θ)′ θ|2 } is bounded by min{|(θb − θ)′ θ|2 , |(θb + θ)′ θ|2 } ≤ min ||θb −
θ||2 , ||θb + θ||2 ≤ kPθb − Pθ k2F . Together with (34), the above bound implies that for each
(θ, τ ) pair in the support of π,
n
1o
≤ β.
Qnθ,τ min{|(θb − θ)′ θ|2 , |(θb + θ)′ θ|2 } >
3

(76)

By Bernstein’s inequality, we have
Qnθ,τ

r
n
n
n
n 1X
1 X b′ 2
2X
′
log n o
2
′b
≤ n−C .
(gi − 1) +
(|θ ǫi | − 1) +
gi (ǫi θ) > C
n
n
n
n
i=1

i=1

i=1

Combining the above analysis and using the assumptions that
N/12, we have

N (log N )5
k4

≤ c and cN ≤ n ≤

′

Qnθ,τ (1 − ψ) ≤ β + n−C .

(77)

Integrating over (θ, τ ) according to the prior π and applying (75), we obtain

 4n
′
.
PH G (1 − ψ) ≤ β + n−C + C e−Ck + N −1 +
1
N

Summing up the Type-I and Type-II errors, we have
PH G ψ + PH G (1 − ψ) ≤ β +
0

1

4n
′
+ C(n−1 + N −1 + e−C k ).
N

(78)

Thus, the proof is complete.

7.2

Proofs of Theorems 5.3, 5.4 and 5.1

Consider a Planted Clique detection problem with N = 12n and 
k = ⌊2s/3⌋. Then the
assumptions of Theorem 5.2 are satisfied. Thus, supQ∈Q(n,s,p,λ) Q kPθb − Pθ k2F > 13 ≤ 41
′
implies a testing error bounded by 41 + 31 + C(n−1 + N −1 + e−C k ), which is smaller than 2/3
as n → ∞. This contradicts the condition (35) that implies the condition of Hypothesis A.
Hence, we must have (36) and the proof of Theorem 5.3 is complete. To prove Theorem 5.4,
note that according to the proof of Theorem 5.2, a testing error bound 1/4 for the sparse PCA
′
problem implies a testing error bound 14 + 31 + C(n−1 + N −1 + e−C k ) for the Planted Clique
detection problem. The same argument we have just used leads to the desired conclusion.
Finally, Theorem 5.1 can be derived from Theorem 5.5 and Theorem 5.2 by applying a similar
argument.
31

7.3

Proof of Theorem 5.5

Let Wi ∼ Np (0, τ θθ ′ + Ip ), then (Xi′ , Yi′ )′ ∼ Np+m (0, Σ) with Σ given in (41). We complete
the proof by noting
o 16 min kb
n
u − uk2 , kb
u + uk2
kPθb − Pθ k2F ≤ 4 min kθb − θk2 , kθb + θk2 ≤
kuk2


3 2
σ 2 (Σx )
L(b
u, u) ≤ 16 1 +
≤ 16 max
L(b
u, u).
2
2
σmin (Σx )

7.4

Proof of Lemma 7.1

√
We first verify that (28)–(29) are proper density functions when |µ| ≤ 3 ηN log N , which is
a corollary of the following lemma.
√
√
Lemma 7.4. If k ≤ N/12, |µ| ≤ 3 ηN log N and |x| ≤ 3 log N , then
4
−1
φ̄µ (x) − φ0 (x) ≤ φ0 (x).
δN
5
Proof. By definition,


µ2 
µ2 
−1
φ̄µ (x) − φ0 (x) = (2δN )−1 φ0 (x) exp µx −
δN
+ exp − µx −
−2 .
2
2
Under the conditions of the lemma, we have |µx| +

µ2
2

≤ 12 , and so



µ2 
4
µ4
µ2 
+ exp − µx −
− 2 ≤ µ2 + |µx|2 +
≤ 8µ2 log N.
exp µx −
2
2
3
3

We complete the proof by combining the last two displays.

The following lemma controls the rescaling constants in (28) and (29).
Lemma 7.5. There exists an absolute constant C > 0 such that for any |µ| ≤ 1, |Mi − 1| ≤
CN −4 for i = 0, 1.
Proof. Note that
Z
Z
1 = fµ,0 (x)dx = M0
−M0

Z


−1
(φ̄µ (x) − φ0 (x)) dx
φ0 (x) − δN

√
|x|>3 log N

= M0 − M0

Z


−1
(φ̄µ (x) − φ0 (x)) dx
φ0 (x) − δN

√
|x|>3 log N

The integral on the RHS is upper bounded by
Z
Z

−1
−1
φ
(x)dx
+
δ
1 + δN
0
N
√


−1
(φ̄µ (x) − φ0 (x)) dx.
φ0 (x) − δN

√
|x|>3 log N

|x|>3 log N

φ̄µ (x)dx ≤ CN −4 ,

where the last inequality comes from standard Gaussian tail bounds. This readily implies
|M0 − 1| ≤ CN −4 The desired bound on M1 follows from similar arguments.
32

Proof of Lemma 7.1. Define
−1
(φ̄µ (x) − φ0 (x)),
gi (x) = φ0 (x) − (−1)i δN

for i = 0, 1.

Then we have for i = 0 and 1,
fµ,i (x) = gi (x) − (1 − Mi 1{|x|≤3√log N } )gi (x),
By Lemma 7.4 and Lemma 7.5,
Z
Z
(1 − Mi 1{|x|≤3√log N } )gi (x) dx
|fµ,i (x) − gi (x)|dx ≤
Z
Z
≤ |1 − Mi | |gi (x)|dx + Mi
√

|x|>3 log N

|gi (x)|dx ≤ CN −3 .

(79)

Therefore, we have


Z
1
1
1
TV
(fµ,0 + fµ,1 ), φ0 =
φ0 (x) − (fµ,0 (x) + fµ,1 (x)) dx
2
2
2
Z
Z
1
1
1 X
≤
|fµ,i (x) − gi (x)|dx ≤ CN −3 ,
φ0 (x) − (g0 (x) + g1 (x)) dx +
2
2
4
i=0,1

where the last inequality is due to the identity φ0 = 12 (g0 + g1 ) and (79). In addition, we
have


1
TV δN fµ,1 + (1 − δN ) (fµ,0 + fµ,1 ) , φ̄µ
2
Z
1 − δN
1
(fµ,0 (x) + fµ,1 (x)) − φ̄µ (x) dx
δN fµ,1 (x) +
=
2
2
Z
1
1 − δN
≤
(g0 (x) + g1 (x)) − φ̄µ (x) dx
δN g1 (x) +
2
2
X 1 − (−1)i δN Z
+
|fµ,i (x) − gi (x)|dx
4
≤ CN

i=0,1
−3

.

Here, the last inequality is due to the identity δN g1 +
completes the proof.

8

1−δN
2

(g0 + g1 ) = φ̄µ and (79). This

Discretization and Computational Lower bounds

To formally address the computational complexity issue in a continuous statistical model,
we adopt the framework in [30]. After introducing the asymptotically equivalent discretized
models, we state the computational lower bounds for sparse PCA and sparse CCA under
the discretized models in Section 8.1. The necessary modifications to Algorithms 1 and 2
are spelled out in Section 8.2 to ensure that they are truly of randomized polynomial time
complexity.
33

For any t ∈ N, define the function [·]t : R → 2−t Z by
 
[x]t = 2−t 2t x .

(80)

For any matrix, the function is defined component-wise. Let
(p,n)

EM


iid
= L(X1 , . . . , Xn ) : Xi ∼ Np (µ, Σ), 1/M ≤ σmin (Σ) ≤ σmax (Σ) ≤ M

be the class of joint distributions of n i.i.d. samples from all multivariate Gaussian distributions with spectrum contained in [1/M, M ], and
(p,n,t)

EM

(p,n)

= {L([X1 ]t , . . . , [Xn ]t ) : L(X1 , . . . , Xn ) ∈ EM

}

be its discretized counterpart. The following lemma bounds the Le Cam distance [26] between
the two classes of distributions. Its proof is given below in Section 8.3.
(p,n)

Lemma 8.1. When 2t t−1/2 ≥ 2(pM )3/2 , the Le Cam distance between EM
(p,n) (p,n,t)
) ≤ n(pM )3/2 t1/2 2−t .
satisfies ∆(EM , EM

(p,n,t)

and EM

For any t ∈ N, we define discretized sparse PCA parameter space as
Qt (n, s, p, λ) = {L([W ]t ) : L(W ) ∈ Q(n, s, p, λ)} ,
and discretized sparse CCA probability space as
P t (n, su , sv , p, m, r, λ; M ) = {L([X]t , [Y ]t ) : L(X, Y ) ∈ P(n, su , sv , p, m, r, λ; M )} .
In view of Theorem 5.1, we are primarily interested in P(n, su , sv , p, m, 1, λ; 3) and its discretized counterpart.
For the sparse PCA parameter spaces, with the choice of n, s, p and λ in Theorem
(p,n)
(p,n,t)
5.3, under condition (35), Q(n, s, p, λ) ⊂ E4
and Qt (n, s, p, λ) ⊂ E4
. Thus, if we
set the discretization level at t = ⌈4 log2 (p + n)⌉, then Q(n, s, p, λ) and Qt (n, s, p, λ) are
asymptotically equivalent. Similarly, with the choice of n, su , sv , p, m, λ in Theorem 5.1, un(p+m,n)
and
der condition (23), when t = ⌈4 log2 (p + m + n)⌉, P(n, su , sv , p, m, 1, λ; 3) ⊂ E5
(p+m,n,t)
t
P (n, su , sv , p, m, 1, λ; 3) ⊂ E5
are also asymptotically equivalent. Therefore, with the
foregoing discretization levels, the statistical difficulties of the original sparse PCA and CCA
problems are asymptotically equivalent to those of the discretized problems. In particular, the
conditions for any procedure to be consistent are the same for the original and the discretized
parameter spaces.

8.1

Computational lower bounds for discretized models

We now state computational lower bounds for the discretized sparse PCA and sparse CCA
problems. The meaning of “randomized polynomial-time estimators” is now based on the
probabilistic Turing machine computation model rather than the computation model mentioned in Remark 5.1.

34

Theorem 8.1. Let t = ⌈4 log2 (p + n)⌉. Under the condition of Theorem 5.3, for any ranb
domized polynomial-time estimator θ,
n
1o 1
> .
(81)
lim inf
sup
Q kPθb − Pθ k2F >
n→∞ Q∈Qt (n,s,p,λ)
3
4
Theorem 8.2. Let t = ⌈4 log2 (p + m + n)⌉. Under the condition of Theorem 5.1, for any
randomized polynomial-time estimator u
b,
n
1 o 1
lim inf
sup
P L(b
u, u) >
> .
(82)
n→∞ P∈P t (n,s ,s ,p,m,1,λ;3)
300
4
u v
To prove these theorems, we need to modify Algorithms 1 and 2 which are not compatible
with the Turing machine computation model. The details are spelled out in the next subsection. After these modifications, the proofs can be obtained by essentially following the lines
of the proofs of their continuous counterparts while controlling some additional negligible
terms in total variation bounds due to additional truncation. The details are omitted.

8.2

Randomized polynomial-time reduction for discretized models

We first introduce a way to approximately sample with polynomial time complexity from a
distribution obtained from discretizing a continuous distribution with density [30, Section
4.2]. The modifications to Algorithms 1 and 2 then follow.
For any w, K ∈ N and K + w + 1 < b ∈ N, define the discrete distribution Aw,b,K [F] with
probability mass function Aw,b,K [f ] as

 R K −w
−2 +i2
f
(x)dx
i − 1
−2K +(i−1)2−w
 ,
(83)
Aw,b,K [f ] − 2K + w = 
R 2K
2
K f (x)dx
−2

b

for i ∈ [2K+w+1 − 1], and let
K

−w

Aw,b,K [f ](2 − 2

) =1−

2K+w+1
X −1
i=1

Aw,b,K [f ](−2K + (i − 1)2−w ).

(84)

In (83), [·]b is the quantization defined previously in (80), and (84) ensures that Aw,b,K [F] is
a proper probability distribution. By the definition of total variation distance, it is straightforward to verify that the approximation error in total variation distance by Aw,b,K [F] to the
distribution of [U 1{U ∈[−2K ,2K ]} ]w with U ∼ F is upper bounded by 2K+w+1−b . As discussed
in Section 4.2 of [30], regardless of the original distribution F, the computational complexity
of drawing a random number from Aw,b,K (F) is O(b2K+w ). This fact is crucial in ensuring
the modified reduction below is of randomized polynomial-time.
Randomized polynomial-time reduction For any t ∈ N, we let
p
w = t + ⌈4 log2 p⌉, K = ⌈log2 (3 log(N + p))⌉, b = w + K + 1 + ⌈4 log2 p⌉.

(85)

The reduction in Algorithm 1 is modified to Algorithm 3 and the reduction in Algorithm 2 is
modified to Algorithm 4. As in the continuous case, a direct reduction from Planted Clique
35

Algorithm 3: Reduction from Planted Clique to Sparse PCA (in Discretized Gaussian
Single Spiked Model)
Input:
1. Graph adjacency matrix A ∈ {0, 1}N ×N ;
2. Estimator θb for the leading eigenvector θ.
Output: A solution to the hypothesis testing problem (22).
e 0 ]. Set
1 Initialization. Generate i.i.d. random variables ξ1 , . . . , ξ2n ∼ Aw,b,K [Φ
1/2

µi = [ηN ]w ξi ,

2

i = 1, . . . , 2n.

Gaussianization. Generate two matrices B0 , B1 ∈ R2n×2n where conditioning on the
µi ’s, all the entries are mutually independent satisfying
L((B0 )ij |µi ) = Aw,b,K [Fµi ,0 ] and

L((B1 )ij |µi ) = Aw,b,K [Fµi ,1 ].

Let A0 ∈ {0, 1}2n×2n be the lower–left 2n × 2n submatrix of the matrix A. Generate a
′ ]′ ∈ R2n×p where for each i ∈ [2n], if j ∈ [2n], we set
matrix W = [W1′ , . . . , W2n
Wij = (B0 )ij (1 − (A0 )ij ) + (B1 )ij (A0 )ij .
3

If 2n < j ≤ p, we let Wij be an independent draw from Aw,b,K [N (0, 1)].
b 1 ]t , . . . , [Wn ]t ) be the estimator of the leading
Test Construction. Let θb = θ([W
eigenvector by treating {[Wi ]t }ni=1 as data. It is normalized to be a unit vector. We
reject H0G if
θb′

2n
 1 X

n

[Wi ]t [Wi ]′t

i=n+1

 
1
θb ≥ 1 + [ k ηN ]w .
w
4

to discretized sparse CCA can be obtained by constructing the estimator θb in the third step
of Algorithm 3 from Algorithm 4.
By (85) and the discussion following (84), the complexity for sampling any random variable in the above reduction is O(p8 (log p)3/2 ), and in total, we need to generate no more than
O(n(p + n)) random variables. Hence, the total complexity for random number generation
is O(p10 (log p)3/2 ) in view of the condition p ≥ 2n. On the other hand, it is straightforward
b have complexity
to verify that all the other computations (except for the estimator u
b or θ)
10
3/2
no more than O(p (log p) ). Since the conditions of Theorems 8.2 and 8.1 ensure that for
some constant a > 1, 2n ≤ p ≤ na and n ≤ N/12, we obtain that the additional computational complexity induced by the proposed reductions is O(N 10a (log N )3/2 ). Therefore, they
are of randomized polynomial-time complexity.

8.3

Proof of Lemma 8.1

We need the following lemma for the proof.

36

Lemma 8.2. For X ∼ Np (µ, Σ) with M −1 ≤ σmin (Σ) ≤ σmax (Σ) ≤ M and U = (U1 , . . . , Up )′
iid

where Ui ∼ Unif[0, 1], we have for any t−1/2 2t ≥ 2(pM )3/2 ,

TV(X, [X]t + 2−t U ) ≤ (pM )3/2 t1/2 2−t .
(p,n,t)

comes from discretizing a corresponding distribution in
Since each distribution in EM
(p,n,t)
(p,n)
(p,n)
) = 0. On the other hand,
EM on a grid with equal spacing 2−t , we have δ(EM , EM
(p,n)
(p,n,t)
, EM ) ≤ n(pM )3/2 t1/2 2−t . This completes
Lemma 8.2 and Lemma 7 of [30] lead to δ(EM
the proof.
−t
Proof of Lemma 8.2. Let f and g denote the density functions of X and [X]
Qpt + 2−t U , respectively. Then g is a piecewise constant function. For any (x1 , ..., xp ) ∈ B = i=1 [2 ij , 2−t (ij +
1)), where ij ∈ Z, we have
Z
1
f (x1 , ..., xp )dx1 ...dxp ,
g(x1 , ..., xp ) =
ν(B) B

where ν is the Lebesgue measure. Hence,
sup
||x−µ||∞ ≤K

≤
≤

g(x)
−1 ≤
f (x)
′

e|(x−µ) Σ

sup

sup
||x−µ||∞ ≤K
||x−y||∞ ≤2−t

f (x)
−1
f (y)

−1 (x−µ)−(y−µ)′ Σ−1 (y−µ)|/2

||x−µ||∞ ≤K
||x−y||∞ ≤2−t

ekΣ

sup

−1

−1 k k(x−µ)(x−µ)′ −(y−µ)(y−µ)′ k /2
F
F

||x−µ||∞ ≤K
||x−y||∞ ≤2−t

−1



≤ exp p3/2 M K2−t − 1
≤

(87)

3 3/2
p M K2−t ,
2

Algorithm 4: Reduction from Discretized Sparse PCA to Discretized Sparse CCA
Input:
1. Observations W1 , . . . , Wn ∈ (2−t Z)p ;
2. Estimator u
b of the first leading canonical correlation coefficient u.
Output: An estimator θb of the leading eigenvector of L(W1 ).
1 Generate i.i.d. random vectors Zi = (Zi1 , . . . , Zip )′ for i ∈ [n] with
iid

Zij ∼ Aw,b,K [N (0, 1)]. Set

h 1 i
(Wi + Zi ),
Xi = √
2 w

2

(86)

h 1 i
Yi = √
(Wi − Zi ),
2 w

Compute u
b=u
b([X1 ]t , [Y1 ]t , . . . , [Xn ]t , [Yn ]t ). Set

b 1 , . . . , Wn ) = [b
θb = θ(W
u/kb
uk]w .
37

i = 1, . . . , n.

whenever p3/2 M K2−t ≤ 21 . The inequality (86) holds since
|(x − µ)′ Σ−1 (x − µ) − (y − µ)′ Σ−1 (y − µ)|


= Tr Σ−1 [(x − µ)(x − µ)′ − (y − µ)(y − µ)′ ]

≤ kΣ−1 kF k(x − µ)(x − µ)′ − (y − µ)(y − µ)′ kF

√
by Cauchy-Schwarz inequality. The inequality (87) holds because kΣ−1 kF ≤ pkΣ−1 kop ≤
√
pM and k(x − µ)(x − µ)′ − (y − µ)(y − µ)′kF ≤ p||x − y||∞ (||x − µ||∞ + ||y − µ||∞) ≤ 2pK2−t .
Note that
Z
Z
Z
g(x)
− 1 dx.
f (x)
|f (x) − g(x)|dx +
|f − g| ≤
f (x)
||x−µ||∞ ≤K
||x−µ||∞ >K
According to Gaussian tail probability, the first term can be bounded by 2p

q

√
2
2 M − (K−1)
2M
.
π K−1 e

The second term is bounded by 23 p3/2 M K2−t according to our previous analysis. Choosing
√
K = 2M t log 2 + 1, we obtain the bound 2(pM )3/2 t1/2 2−t for Rall t−1/2 2t ≥ 2(pM )3/2 . The
conclusion follows the simple fact that TV(X, [X]t + 2−t U ) = 21 |f − g|.

9
9.1

Additional Proofs
Proof of Theorem 3.1

We first present a bound for the estimator defined by (9) under the joint loss.
Theorem 9.1. Assume (13) for some sufficiently small c > 0. Then there exist constants
C, C ′ > 0 only depending on c such that




b (0) (Vb (0) )′ − U V ′ Σ1/2 k2 ≤ C r(su + sv ) + su log ep + sv log em ,
kΣ1/2
U
y
F
x
nλ2
su
sv

with P-probability at least 1 − exp (−C ′ (su + log(ep/su ))) − exp (−C ′ (sv + log(em/sv ))) uniformly over P ∈ P(n, su , sv , p, m, r, λ; M ).
Theorem 9.1 is similar to Theorem 1 of [18], except that the loss function depends on the
marginal covariances so that the error bound is independent of M . Its proof is omitted given
the similarity with that of Theorem 1 of [18].
∗
b (1)
b (1)
b (1)
From now on, we omit the superscript in Σ
x , Σy and Σxy for simplicity. Let U =
b (1) − U ∗ . The following lemmas are needed in the proof of Theorem
U ΛV ′ Σy Vb (0) and ∆ = U
3.1. Their proofs are given in Section 9.3.2.
u)
Lemma 9.1. Assume su log(ep/s
≤ c for some sufficiently small constant c ∈ (0, 1). Then,
n
there exist some constants C, C ′ > 0 only depending on c, such that with probability at least
1 − exp(−C ′ su log(ep/su )),
r
r
s
log(ep/s
)
su log(ep/su )
u
u
b (0) kop ≤ 1 + C
, k(Vb (0) )′ Σy Vb (0) − Ikop ≤ C
.
kΣ1/2
y V
n
n

38

u)
Lemma 9.2. Assume su log(ep/s
≤ c for some sufficiently small constant c > 0. Then, there
n
′
exist some constants C, C > 0 only depending on c, such that

′
2
′
1/2
2
b 1/2 2
(1 − δC
)kΣ1/2
x ∆kF ≤ kΣx ∆kF ≤ (1 + δC )kΣx ∆kF ,
q
u)
′
′
.
with probability at least 1 − exp(−C su log(ep/su )), with δC = C su log(ep/s
n

Lemma 9.3. Assume n1 (sv log(em/sv ) + su log(ep/su ) + rsu ) ≤ c for some sufficiently small
constant c > 0. Then, there exist some constants C, C ′ > 0 only depending on c, such that
r


rsu + su log(ep/su ) 1/2
′ b
(0)
kΣx ∆kF ,
Tr ∆ (Σxy − Σxy )Vb
≤C
n
with probability at least 1 − exp (−C ′ (su log(ep/su ) + rsu )) − exp(−C ′ sv log(em/sv )).

Lemma 9.4. Assume n1 (sv log(em/sv ) + su log(ep/su ) + rsu ) ≤ c for some sufficiently small
constant c > 0. Then, there exist some constants C, C ′ > 0 only depending on c, such that
r


b x − Σx )U ∗ ≤ C rsu + su log(ep/su ) kΣ1/2 ∆kF ,
Tr ∆′ (Σ
x
n
with probability at least 1 − exp (−C ′ (su log(ep/su ) + rsu )) − exp(−C ′ sv log(em/sv )).
1/2

Proof. The proof consists of two steps. First, we derive a bound for kΣx ∆kF . Next, we
b , U ).
derive the desired bound for L(U
Step 1. By the definition of the estimator, we have

b (1) )′ Σ
b xU
b (1) ) − 2 Tr((U
b (1) )′ Σ
b xy Vb (0) ) ≤ Tr((U ∗ )′ Σ
b x U ∗ ) − 2 Tr((U ∗ )′ Σ
b xy Vb (0) ).
Tr((U

After rearrangement, we have


b x ∆) ≤ 2 Tr ∆′ (Σ
b xy Vb (0) − Σ
b xU ∗ )
Tr(∆′ Σ


b xy − Σxy )Vb (0) + 2 Tr ∆′ (Σ
b x − Σx )U ∗ .
≤ 2 Tr ∆′ (Σ

Using Lemma 9.2, Lemma 9.3 and Lemma 9.4, we have
r
1 1/2 2
rsu + su log(ep/su ) 1/2
kΣx ∆kF ≤ 4C
kΣx ∆kF ,
2
n
1/2

with high probability, which immediately implies a bound for kΣx ∆k2F . This completes
Step 1.
Step 2. We claim that

−1
′
b (0) ≤ C ,
σmin
Σ1/2
x U ΛV Σy V
λr
b − Σ1/2 U
b (1) ((U
b (1) )′ Σx U
b (1) )−1/2 kF ≤ C rsu + su log(ep/su ) ,
kΣx1/2 U
x
n
39

(88)
(89)

with high probability. The two claims (88) and (89) will be proved in the end. We bound
b , U ) by
L(U
q
b , U ) = inf kΣ1/2
b
L(U
x (U W − U )kF
W ∈O(r,r)

1/2 b (1) b (1) ′
b
b (1) )−1/2 kF
≤ kΣ1/2
((U ) Σx U
x U − Σx U
b (1) ((U
b (1) )′ Σx U
b (1) )−1/2 W − Σ1/2 U kF .
+ inf kΣ1/2 U
W ∈O(r,r)

x

x

With high probability, we could further bound the rightmost side by
r
1
rsu + su log(ep/su )
C
+ √ kPΣ1/2 Ub − PΣ1/2 U kF
x
x
n
2
r


rsu + su log(ep/su )
′
−1
b (0) kΣ1/2
≤ C
Σ1/2
+ Cσmin
x ∆kF
x U ΛV Σy V
n
r
rsu + su log(ep/su )
≤ C
.
nλ2

(90)
(91)

The bound (90) is due to the claim (89), Lemma 6.6 and the fact that PΣ1/2 Ub = PΣ1/2 Ub (1) .
The inequality (91) is derived from the sin-theta theorem [41]. Thus, we have obtained
b , U ). To finish the proof, we need to prove (88) and (89). Since
the desired bound for L(U
1/2
Σx U ∈ O(p, r), we have

−1
′
b (0) ≤ λ−1 k(V ′ Σy Vb (0) )−1 kop .
σmin
Σ1/2
x U ΛV Σy V
Thus, it is sufficient to bound k(V ′ Σy Vb (0) )−1 kop . By Theorem 9.1 and sin-theta theorem [41],
kPΣ1/2 Vb (0) − PΣ1/2 V kF is sufficiently small. In view of Lemma 6.6, there exists W ∈ O(r, r),
y
y
such that
kΣy1/2 Vb (0) (Vb (0) Σy Vb (0) )−1/2 − Σ1/2
y V W kF
is sufficiently small. Therefore, together with Lemma 9.1,

kV ′ Σy Vb (0) − W kop
≤ kV ′ Σy Vb (0) − V ′ Σy V W (Vb (0) Σy Vb (0) )1/2 kop + kW kop k(Vb (0) Σy Vb (0) )1/2 − Ikop
b (0) Σy Vb (0) )1/2 kop + k(Vb (0) Σy Vb (0) )1/2 − Ikop
≤ kΣy1/2 Vb (0) (Vb (0) Σy Vb (0) )−1/2 − Σ1/2
y V W kF k(V

is also sufficiently small. By Weyl’s inequality [20, p.449], |σmin (V ′ Σy Vb (0) )−1| ≤ kV ′ Σy Vb (0) −
W kop is sufficiently small. Hence, k(V ′ Σy Vb (0) )−1 kop ≤ 2 with high probability, which implies
the desired bound in (88). Finally, we need to prove (89). We have
b − Σ1/2
b (1) ((U
b (1) )′ Σx U
b (1) )−1/2 kF
kΣx1/2 U
x U
b (1) kF k((U
b (1) )′ Σ
b (2) U
b (1) )−1/2 − ((U
b (1) )′ Σx U
b (1) )−1/2 kop
≤ kΣx1/2 U
x


b (1) )′ (Σ
b (2)
b (1) kop .
≤ C kΣx1/2 U ΛV ′ Σy Vb (0) kF + kΣ1/2
x − Σx )U
x ∆kF k(U
1/2

1/2

We have already shown that kΣx ∆kF is sufficiently small. The term kΣx U ΛV ′ Σy Vb (0) kF
√
√
√
is bounded by rkV ′ Σy Vb (0) kop ≤ r(1 + kV ′ Σy Vb (0) − W kop ) ≤ C r by using the bound
40

b (1) )′ (Σ
b (2)
b (1) kop , note that Σ
b (2)
derived for kV ′ Σy Vb (0) − W kop . To bound k(U
x − Σx )U
x only
(1)
b
depends on D2 and is independent of U . Using union bound and an ǫ-net argument (see,
1/2
for example, [36]) and the fact q
that r ≤ su (which is implied by Σx U ∈ O(p, r)), we have
b (1) )′ (Σ
b (2)
b (1) kop ≤ C rsu +su log(ep/su ) with high probability. Hence, the proof is
k(U
x − Σx )U
n
complete.

9.2

Proof of Theorem 3.2

For any probability
measures P, Q, define the Kullback-Leibler(KL) divergence by D(P||Q) =

R
dP
log dQ dP. The following result is Lemma 14 in [18]. It gives explicit formula for the KL
divergence between distributions generated by a special kind of covariance matrices.
"
#
′
Ip
λU(i) V(i)
Lemma 9.5. For i = 1, 2, let Σ(i) =
with λ ∈ (0, 1), U(i) ∈ O(p, r)
′
λV(i) U(i)
Im
and V(i) ∈ O(m, r). Let P(i) denote the distribution of a random i.i.d. sample of size n from
the Np+m (0, Σ(i) ) distribution. Then
D(P(1) ||P(2) ) =

nλ2
′
kU V ′ − U(2) V(2)
k2F .
2(1 − λ2 ) (1) (1)

The main tool for our proof is the following Fano’s lemma [44, Lemma 3].
Proposition 9.1. Let (Θ, ρ) be a metric space and {Pθ : θ ∈ Θ} a collection of probability
measures. For any totally bounded T ⊂ Θ, denote by M(T, ρ, ǫ) the ǫ-packing number of T
with respect to ρ, i.e., the maximal number of points in T whose pairwise minimum distance
in ρ is at least ǫ. Define the KL diameter of T by dKL (T ) , supθ,θ′ ∈T D(Pθ ||Pθ′ ). Then


 ǫ2
dKL (T ) + log 2
2 b
.
inf sup Pθ ρ θ(X), θ ≥
≥1−
4
log M(T, ρ, ǫ)
b
θ θ∈Θ

(92)

Finally, we lower bound the prediction loss by the squared subspace distance. Its proof
is given in Section 9.3.2.
Proposition 9.2. Suppose the eigenvalues of Σx lie in the interval [M1 , M2 ]. Then, we have
kPUb − PU k2F ≤
A similar inequality holds for L(Vb , V ).

√ M2
b , U ).
2
L(U
M1

Proof of Theorem 3.2. Let us first give an outline of the proof. By Proposition 9.2, we have


b , U ) ≥ Cǫ2 ≥ inf sup P kP b − PU k2F ≥ C1 ǫ2 ,
inf sup P L(U
U
b P∈P
U

b P∈P
U

for any rate ǫ2 . Therefore, it is sufficient to derive a lower bound for the loss kPUb − PU k2F .
Without loss of generality, we assume su /3 is an integer and su ≤ 3p/4. The case su >

41

3p/4 is harder and thus it shares the same lower bound. The subset of covariance class
F(p, m, su , sv , r, λ; M ) we consider is





e 0
Ip
λU V0′
U
e ∈ B,
:U =
T = Σ=
,U
λV0 U ′
Im
0 ur

p−2su /3
ur ∈ R
, ||ur || = 1, |supp(ur )| ≤ su /3 ,


′
where V0 = Ir 0′ ∈ O(m, r) and B is a subset of O(2su /3, r − 1) to be specified later.
e and the vector ur . As U
e and ur vary, we
From the construction, U depends on the matrix U
∗
always have U ∈ O(p, r). We use T (ur ) to denote a subset of T where ur = u∗r is fixed, and
e ∗ ) to denote a subset of T where U
e =U
e ∗ is fixed.
use T (U
rsu
∗
The proof has three steps. In the first step, we derive the part nλ
2 using the subset T (ur )
s
log(ep/s
)
u
u
for some particular u∗r . In the second step, we derive the other part
using the
nλ2
∗
∗
e
e
subset T (U ) for some fixed U . Finally, we combine the two results in the third step.


e0 = Ir−1 0′ ′ ∈
Step 1. Let u∗r = (1, 0, ..., 0)′ , and we consider the subset T (u∗r ). Let U
√
O(2su /3, r − 1) and ǫ0 ∈ (0, r] to be specified later. Define
By Lemma 9.5,

e ∈ O(2su /3, r − 1) : kU
e −U
e0 kF ≤ ǫ0 }.
B = B(ǫ0 ) = {U

dKL (T (u∗r )) =

2 2
nλ2
e(1) − U
e(2) k2 ≤ 2nλ ǫ0 .
k
U
F
2(1 − λ2 )
1 − λ2
∈B(ǫ0 )

sup
e(i)
U

(93)

Here, the equality is due to the definition of V0 and the inequality due to the definition of
B(ǫ0 ). We now establish a lower bound for the packing number of T (u∗r ). For some α ∈ (0, 1)
e(1) , . . . , U
e(N ) } ⊂ O(2su /3, r − 1) be a maximal set such that for
to be specified later, let {U
any i 6= j ∈ [N ],
√
e(i) U
e′ − U
e0 U
e0′ kF ≤ ǫ0 ,
e(i) U
e′ − U
e(j) U
e ′ kF ≥ 2αǫ0 .
(94)
kU
kU
(i)
(i)
(j)
Then by [13, Lemma 1], for some absolute constant C > 1,
N≥



1
Cα

(r−1)(2su /3−r+1)

.

e(i) U
e′ −
It is easy to see that the loss function kPU(i) − PU(j) k2F on the subset T (u∗r ) equals kU
(i)
√
e(j) U
e ′ k2 . Thus, for ǫ = 2αǫ0 with sufficiently small α, log M(T (u∗ ), ρ, ǫ) ≥ (r−1)(2su /3−
U
r

(j) F

1
1
1
1
r+1) log Cα
≥ (r−1)( 16 su −1) log Cα
≥ 12
rsu log Cα
when r is sufficiently large and r ≤ su /2.
rsu
2
Taking ǫ0 = c1 nλ2 for sufficiently small c1 , we have



ǫ2
2
inf sup P kPUb − PU kF ≥
≥1−
4
b T (u∗ )
U
r

2c1 rsu
1−λ2 + log 2
1 .
1
12 rsu log Cα

(95)

Since λ is bounded away from 1, we may choose sufficiently small c1 and α, so that the right
hand side of (95) can be lower bounded by 0.9. This completes the first step.
42

u)
can be obtained from the rank-one argument spelled out in
Step 2. The part su log(ep/s
nλ2


e ∗ ) with U
e ∗ = Ir−1 0′ ′ ∈ O(2su /3, r − 1).
[14]. To be rigorous, consider the subset T (U
e ∗ ), the loss function is
Restricting on the set T (U

kPU(i) − PU(j) k2F = kur,(i) u′r,(i) − ur,(j) u′r,(j) k2F .

Let X = [X1 X2 ] with X1 ∈ Rn×(r−1) and X2 ∈ Rn×(p−r+1) , and Y = [Y1 Y2 ] with
Y1 ∈ Rn×(r−1) and Y2 ∈ Rn×(m−r+1) . Then it is further equivalent to estimating u1 under projection loss based on the observation (X2 , Y2 ), because (X2 , Y2 ) is a sufficient statistic
for ur . Applying the argument in [14, Appendix G] and choosing the appropriate constant,
we have


su log(ep/su )
∧
c
inf sup P kPUb − PU k2F ≥ C
(96)
0 ≥ 0.9,
nλ2
b
U
e ∗)
T (U
for some constant C > 0. This completes the second step.
Step 3. For any P ∈ P, by union bound, we have

P kPUb − PU k2F ≥ ǫ21 ∨ ǫ22


≥ 1 − P kPUb − PU k2F < ǫ21 − P kPUb − PU k2F < ǫ22


= P kPUb − PU k2F ≥ ǫ21 + P kPUb − PU k2F ≥ ǫ22 − 1.

rsu
Taking supT (u∗ )∪T (Ue ∗ ) on both sides of the inequality, and letting ǫ21 = C1 nλ
2 in (95) and
r

u)
ǫ22 = C2 su log(ep/s
∧ c0 in (96), we have
nλ2


sup P kPUb − PU k2F ≥ ǫ21 ∨ ǫ22 ≥ 0.9 + 0.9 − 1 = 0.8,

P∈P

e
where we have used the identity supU∈T
e
e ∗ ) (f (ur ) + g(U )) = supur ∈T (U
e ∗ ) f (ur ) +
(u∗r ),ur ∈T (U
e
supUe ∈T (u∗ ) g(U ). Careful readers may notice that we have assume sufficiently large r in Step
r
1. For r which is not sufficiently large, a similar rank-one argument as in Step 2 gives the
desired lower bound. Thus, the proof is complete.

9.3

Proofs of technical lemmas

This section gathers the proofs of all technical results used in the above sections. The proofs
are organized according to the order of their first appearance. To simplify notation, we denote
b x(j) , Σ
b (j)
b (j)
b b
b
Σ
y and Σxy by Σx , Σy and Σxy for j ∈ {0, 1, 2} whenever there is no confusion from
the context.
9.3.1

Proofs of lemmas in Section 6

In order to prove Lemma 6.1, we need an auxiliary result.

43

Lemma 9.6. Assume n1 (su + sv + log(ep/su ) + log(em/sv )) ≤ c for some sufficiently small
constant c ∈ (0, 1). Then there exist some constants C, C ′ > 0 only depending on c such that
s 

ep
1
′b
′b
1/2
kU Σx U − Ikop ∨ k(U Σx U ) − Ikop ≤ C
,
su + log
n
su
s 

em
1
′b
′b
1/2
,
kV Σy V − Ikop ∨ k(V Σy V ) − Ikop ≤ C
sv + log
n
sv
with probability at least 1 − exp(−C ′ (su + log(ep/su ))) − exp(C ′ (sv + log(em/sv ))).
Proof. Using the definition of operator norm and the sparsity of U , we have
b x U − Ir kop = kU ′ (Σ
b x − Σx )U kop = k(USu ∗ )′ (Σ
b xSu Su − ΣxSu Su )USu ∗ kop
kU ′ Σ
−1/2
b xSu Su − ΣxSu Su )(USu ∗ v) ≤ kΣ1/2 USu ∗ k2 kΣ−1/2 Σ
b
= sup (USu ∗ v)′ (Σ
op
xSu Su xSu Su ΣxSu Su − Ikop ,
xSu Su
||v||=1

1/2

−1/2

−1/2

b xSu Su Σ
where kΣxSu Su USu ∗ k2op ≤ 1 and kΣxSu Su Σ
xSu Su − Ikop is bounded by the desired rate
b x U )1/2 −
with high probability according to Lemma 16 in [18]. Lemma 15 in [18] implies k(U ′ Σ
b x U − Ikop , and thus k(U ′ Σ
b x U )1/2 − Ikop also shares same upper bound. The
Ikop ≤ CkU ′ Σ
′
′
b y V − Ikop ∨ k(V Σ
b y V )1/2 − Ikop can be derived by the same argument.
upper bound for kV Σ
Hence, the proof is complete.
Proof of Lemma 6.1. According to the definition, we have

1/2
′b
1/2
e
b x U )−1/2 kop ,
kΣ1/2
− Ikop k(U ′ Σ
x (U − U )kop ≤ kΣx U kop k(U Σx U )
1/2
′b
1/2
e
b y V )−1/2 kop ,
kΣ1/2
− Ikop k(V ′ Σ
y (V − V )kop ≤ kΣy V kop k(V Σy V )
e − Λkop ≤ k(U ′ Σ
b x U )1/2 − Ikop kΛ(V ′ Σ
b y V )1/2 kop
kΛ

b y V )1/2 − Ikop .
+kΛkop k(V ′ Σ

Applying Lemma 9.6, the proof is complete.

e , we have U
e ′Σ
b xU
e = I, and thus Σ
b 1/2
e
Proof of Lemma 6.2. By the definition of U
x U ∈ O(p, r).
b 1/2
e
Similarly Σ
y V ∈ O(m, r). Thus,
b 1/2
e b 1/2
b 1/2 e
b 1/2 e
kΣ
x AΣy kop ≤ kΣx U kop kΣy V kop ≤ 1.

(97)

b y V )−1 (V ′ Σ
b y V )) = r.
Tr(Q′ Q) = Tr((V ′ Σ

(98)

′
b 1/2
e b 1/2
e
Now let us use the notation Q = Σ
x AΣy . Then, by the definition of A, we have Q Q =
′b
−1 ′ b 1/2
b 1/2
Σ
y V (V Σy V ) V Σy , and

Combining (97) and (98), it is easy to see that all eigenvalues of Q′ Q are 1. Thus, we have
kQk∗ = r and kQkop = 1. The proof is complete.

44

Proof of Lemma 6.3. Denote F = [f1 , ..., fr ], G = [g1 , ..., gr ] and cj = fj′ Egj . By kEkop ≤ 1,
we have |cj | ≤ 1. The left hand side of (44) is lower bounded by hF KG′ , F G′ − Ei ≥
hF DG′ , F G′ − Ei − P
kK − DkF kF G′ − EkF , where hF DG′ , F G′ − Ei = hD, I − F ′ EGi =
P
r
r
l=1 dl (1 − cl ) ≥ dr
l=1 (1 − cl ). The first term on the right hand side of (44) is

dr 
dr
kF G′ k2F + kEk2F − 2 Tr(F ′ EG)
kF G′ − Ek2F =
2
2
r

X
dr 
≤
cj
Tr(F ′ F G′ G) + kEkop kEk∗ − 2
2
j=1

≤ dr

r
X
j=1

(1 − cj ).

This completes the proof.
b xy − Σ
e xy ||∞ can be upper bounded by the
Proof of Lemma 6.4. Using triangle inequality, ||Σ
following sum,
b xy − Σxy ||∞ + ||(Σ
b x − Σx )U ΛV ′ Σy ||∞
||Σ
b y − Σy )||∞ + ||(Σ
b x − Σx )U ΛV ′ (Σ
b y − Σy )||∞ .
+||Σx U ΛV ′ (Σ

The first term can be bounded by the desired rate by union bound and Bernstein’s inequality
[36, Prop. 5.16]. For the second term, it can be written as
n

1X
max
(Xij [Xi′ U ΛV ′ Σy ]k − EXij [Xi′ U ΛV ′ Σy ]k ) ,
j,k n
i=1

where Xij is the j-th element of Xi and the notation [·]k means the k-th element of the
referred vector. Thus, it is a maximum over average of centered sub-exponential random
variables. Then, by Bernstein’s inequality and union bound, it is also bounded by the desired
rate.
we can bound the third term. For the last term, it can be bounded by
Pr Similarly,
′ b
′ b
b
b
λ
||(
Σ
−
Σ
x
x )ul vl (Σy − Σy )||∞ , where for each l, ||(Σx − Σx )ul vl (Σy − Σy )||∞ can be
l=1 l
written as
n
n
1 X
 1 X

max
(Xij Xi′ ul − EXij Xi′ ul )
(Yik Yi′ vl − EYik Yi′ vl ) .
j,k
n
n
i=1

i=1

with the desired probability using union bound and
It can be bounded by the rate log(p+m)
n
. Under the
Bernstein’s inequality. Hence, the last term can be bounded by λ1 r log(p+m)
n
q
log(p+m)
assumption that r
is bounded by a constant, it can further be bounded by the
n
q
rate log(p+m)
with high probability. Combining the bounds of the four terms, the proof is
n
complete.
Proof of Lemma 6.6. By the property of least squares, we have
inf kF − GW k2F = kF − G(G′ G)− G′ F k2F = kF − PG F k2F = r − Tr(PF PG ).
W

Since kPF − PG k2F = 2r − 2 Tr(PF PG ), the proof is complete.
45

Proof of Lemma 6.7. By the definition of U ∗ , we have Σxy Vb (0) = Σx U ∗ . Thus,

b xy Vb (0) − Σ
b x U ∗ ]j· || ≤ max ||[(Σ
b xy − Σxy )Vb (0) ]j· || + max ||[(Σ
b x − Σx )U ∗ ]j· ||.
max ||[Σ

1≤j≤p

1≤j≤p

1≤j≤p

b x − Σx )U ∗ ]j· ||. Note that the sample covariance can be
Let us first bound max1≤j≤p ||[(Σ
written as
n

 X
1/2 1
b
Zi Zi′ Σ1/2
Σx = Σx
x ,
n
i=1

where
1/2

{Zi }ni=1

are i.i.d. Gaussian vectors distributed as N (0, Ip ). Let Tj′ be the j-th row of

Σx , and then we have
n

b x − Σx )U ∗ ]j· =
[(Σ

1X ′
∗
′ 1/2 ∗
(Tj Zi Zi′ Σ1/2
x U − Tj Σx U ).
n
i=1

For each i and j, define vector

(j)
Wi

1/2

#
Tj′ Zi
.
=
1/2
(U ∗ )′ Σx Zi
"

(j)

(j)

b x −Σx )U ∗ ]j· || ≤ k 1
Since Tj′ Zi Zi′ Σx U ∗ is a submatrix of Wi (Wi )′ , we have ||[(Σ
n
(j)

(j)

EWi (Wi )′ )kop . Hence, for any t > 0, we have
n
o
b x − Σx )U ∗ ]j· || > t
P max ||[(Σ

Pn

(j)
(j) ′
i=1 (Wi (Wi ) −

1≤j≤p

p
n
o
n 1X
X
(j)
(j)
(j)
(j)
(Wi (Wi )′ − EWi (Wi )′ )kop > t
P k
≤
n

≤
(j)

j=1
p
X
j=1

i=1


n
exp C1 r − C2 n min

t

t2

,
kW (j) kop kW (j) k2op

o
,

(99)

(j)

where W (j) = EWi (Wi )′ , and we have used concentration inequality for sample covariance
[36, Thm. 5.39]. Since kW (j) kop ≤ C3 for some constant C3 only depending on M , (99) can
be bounded by


exp C1′ (r + log p) − C2′ n(t ∧ t2 ) .

p
Take t2 = C4 r+log
for some sufficiently large C4 , and under the assumption n−1 (r + log p) ≤
n
q
b x −Σx )U ∗ ]j· || ≤ C r+log p with probability at least 1−exp(−C ′ (r+log p)).
C1 , max1≤j≤p ||[(Σ
n

b xy − Σxy )Vb (0) ]j· ||. Let us sketch the
Similar arguments lead to the bound of max1≤j≤p ||[(Σ
proof. Note that we may write
n

1X ′
(0)
b
b
[(Σxy − Σxy )V ]j =
Tj Zi Yi′ Vb (0) − E(Tj′ Zi Yi′ Vb (0) ) .
n
i=1

46

(j)

Then, define Hi

=




Tj′ Zi
, and we have
(Vb (0) )′ Yi

b xy − Σxy )Vb (0) ]j || ≤ max k 1
max ||[(Σ
1≤j≤p n
1≤j≤p

n
X
i=1

(j)

Using the same argument, we can bound this term by C
1 − exp(−C ′ (r + log p)). Thus, the proof is complete.
9.3.2

(j)

(j)

(j)

(Hi (Hi )′ − EHi (Hi )′ )kop .
q

r+log p
n

with probability at least

Proofs of lemmas in Section 9
1/2

(0)

Proof of Lemma 9.1. Let Tv = Sbv ∪Sv , where Sbv = supp(Vb (0) ). First, let us bound kΣyTv Tv VbTv ∗ kop .
b y Vb (0) = Ir , we have
Since (Vb (0) )′ Σ
1/2
(0)
1/2
b −1/2
b 1/2 Vb (0) kop ≤ kΣ1/2 Σ
b −1/2 kop kΣ
kΣyTv Tv VbTv ∗ kop ≤ kΣyTv Tv Σ
yTv Tv yTv Tv kop
yTv Tv
yTv Tv

−1/2 −1/2
1/2
b
b −1 Σ1/2 k1/2 = σmin (Σ−1/2 Σ
= kΣyTv Tv Σ
yTv Tv yTv Tv ΣyTv Tv )
yTv Tv yTv Tv op
r
−1/2

su log(ep/su )
−1/2 b
−1/2
≤ 1 − kΣyTv Tv ΣyTv Tv ΣyTv Tv − Ikop
≤ 1+C
,
n
with probability at least 1 − exp(−C ′ su log(ep/su )), where the last inequality is by Lemma
12 of [18]. Hence,

b yTv Tv )Vb (0) kop
b y )Vb (0) kop = k(Vb (0) )′ (ΣyTv Tv − Σ
k(Vb (0) )′ Σy Vb (0) − Ikop = k(Vb (0) )′ (Σy − Σ
Tv ∗
Tv ∗
r
su log(ep/su )
−1/2 b
−1/2
(0)
1/2
,
≤ kΣyTv Tv VbTv ∗ k2op kΣyTv Tv Σ
yTv Tv ΣyTv Tv − Ikop ≤ 4C
n
with probability at least 1 − exp(−C ′ su log(ep/su )). The proof is completed by realizing
1/2
(0)
1/2
kΣyTv Tv VbTv ∗ kop = kΣy Vb (0) kop .

b ). Using the definition of FrobeProof of Lemma 9.2. Let Tu = Sbu ∪ Su , where Sbu = supp(U
nius norm, we have
2
′ b
′ b
b 1/2 2
kΣ1/2
x ∆kF − kΣx ∆kF = Tr(∆ (Σx − Σx )∆) = Tr((∆Tu ∗ ) (ΣxTu Tu − ΣxTu Tu )∆Tu ∗ )
r
su log(ep/su ) 1/2 2
−1/2 b
−1/2
≤ kΣxTu Tu ∆Tu ∗ k2F kΣxTu Tu Σ
kΣx ∆kF ,
xTu Tu ΣxTu Tu − Ikop ≤ C
n
1/2

with high probability, where we have used kΣxTu Tu ∆Tu ∗ k2F = kΣx ∆k2F and Lemma 12 in
[18] in the last inequality. After rearrangement, the proof is complete.
b x is constructed from D0 and Σ
b y is constructed from D1 .
Proof of Lemma 9.3. In this proof, Σ
b ) and Sbv = supp(Vb (0) ).
We use the notation Tu = Su ∪Sbu and Tv = Sv ∪Sbv , where Sbu = supp(U
Note that Tu depends on D1 and Tv depends on D0 . We first condition on D0 , and then we
have

b xy − Σxy )Vb (0) = hΣ
b xyTu Tv − ΣxyTu Tv , ∆Tu ∗ (Vb (0) )′ i
Tr ∆′ (Σ
Tv ∗
1/2

(0)

1/2

−1/2

−1/2

b xyTu Tv − ΣxyTu Tv )Σ
≤ kΣyTv Tv VbTv ∗ kop kΣxTu Tu ∆Tu ∗ kF hΣxTu Tu (Σ
yTv Tv , KTu i

1/2
(0)
1/2
−1/2 b
−1/2
≤ kΣyTv Tv VbTv ∗ kop kΣxTu Tu ∆Tu ∗ kF sup hΣxT T (Σ
xyT Tv − ΣxyT Tv )ΣyTv Tv , KT i
T

47

where T ranges over all subsets with cardinality bounded by 2su , and for each such T ,
1/2
1/2
(0)
1/2
b (0) ′ 1/2
KT = kΣxT T ∆T ∗ (VbTv ∗ )′ ΣyTv Tv k−1
F ΣxT T ∆T ∗ (VTv ∗ ) ΣyTv Tv satisfying kKT kF = 1. We do not
put Tv in the subscript of K because conditioning on D0 , Tv is fixed. For each T , we can
−1/2 b
−1/2
use Lemma 7 in [18] to bound hΣxT T (Σ
xyT Tv − ΣxyT Tv )ΣyTv Tv , KT i . A direct union bound
argument leads to
r
rsu + su log(ep/su )
−1/2 b
−1/2
sup hΣxT T (ΣxyT Tv − ΣxyT Tv )ΣyTv Tv , KT i ≤ C
,
n
T

(0)
1/2
with probability at least 1−exp (−C ′ (su log(ep/su ) + rsu )). By Lemma 9.1, we have kΣyTv Tv VbTv ∗ kop =
1/2
1/2
1/2
kΣy Vb (0) kop ≤ 2 with high probability. Finally, observing that kΣxTu Tu ∆Tu ∗ kF = kΣx ∆kF ,
we have completed the proof.

Proof of Lemma 9.4. It is omitted due to similarity to that of Lemma 9.3.

Proof of Proposition 9.2. Let the singular value decomposition of U be U = ΘRH ′ . Then
we have HRΘ′ Σx ΘRH ′ = U ′ Σx U = I, from which we derive Θ′ Σx Θ = R−2 . Using Lemma
6.6, we have
√
b W − ΘkF
kPUb − PU kF = 2 inf kU
W
√
√
b W HR−1 − ΘRH ′ HR−1 kF ≤ 2 inf kU
b W − U kF kR−1 kop
≤ 2 inf kU
W
W
√
−1/2
′
1/2
b
inf kΣ1/2
≤ 2M1
x (U W − U )kF kΘ Σx Θkop
W
√
b
≤ 2(M2 /M1 )1/2 inf kΣ1/2
x (U W − U )kF
W
√
b
≤ 2(M2 /M1 )1/2 inf kΣ1/2
x (U W − U )kF .
W ∈O(r,r)

1/2

b W − U )k2 = Tr((U
b W − U )′ Σx (U
b W − U )), the proof is complete.
Finally, by kΣx (U
F

10

Implementation of (18)

To implement the convex programming (18), we turn to the Alternating Direction Method
b x and Σ
b y for Σ
b (0)
of Multipliers (ADMM) [16, 11]. In the rest of this section, we write Σ
x and
(0)
b
Σy for notational convenience.
First, note that (18) can be rewritten as
minimize
subject to
where

f (F ) + g(G),
b 1/2 F Σ
b 1/2 − G = 0,
Σ
x

b xy , F i + ρkF k1 ,
f (F ) = −hΣ

g(G) = ∞1{kGk∗ >r} + ∞1{kGkop >1} .
48

(100)

y

(101)
(102)

Thus, the augmented Lagrangian form of the problem is
b 1/2 F Σ
b 1/2 − Gi + η kΣ
b 1/2 F Σ
b 1/2 − Gk2 .
Lη (F, G, H) = f (F ) + g(G) + hH, Σ
x
y
y
F
2 x

(103)

Following the generic algorithm spelled out in Section 3 of [11], suppose after the kth iteration, the matrices are (F k , Gk , H k ), then we update the matrices in the (k + 1)-th
iteration as follows:
F k+1 = argmin Lη (F, Gk , H k ),

(104)

Gk+1 = argmin Lη (F k+1 , G, H k ),

(105)

k+1 b 1/2
b 1/2
H k+1 = H + η(Σ
Σy − Gk+1 ).
x F

(106)

F

G
k

The algorithm iterates over (104) – (106) till some convergence criterion is met. It is clear
that the update (106) for the dual variable H is easy to calculate. Moreover the updates
(104) and (105) can be solved easily and have explicit meaning in giving solution to sparse
CCA. We are going to show that (104) can be viewed as a Lasso problem. Thus, this step
targets at the sparsity of the matrix U V ′ . The update (105) turns out to be equivalent to
a singular value capped soft thresholding problem, and it targets at the low-rankness of the
1/2
1/2
matrix Σx U V ′ Σy . In what follows, we study in more details the updates for F and G.
First, we note that (104) is equivalent to
b 1/2 F Σ
b 1/2 − Gk k2
b 1/2 F Σ
b 1/2 i + η kΣ
F k+1 = argmin f (F ) + hH k , Σ
x
y
F
x
y
2
F
1 k 1 b −1/2 b b −1/2 2
η b 1/2 b 1/2
k
Σxy Σy )kF + ρkF k1 .
= argmin kΣ
x F Σy − (G − H + Σx
2
η
η
F

(107)

Thus, it is clear that the update of F in (104) can be viewed as a Lasso problem as summarized
in the following proposition. Here and after, for any positive semi-definite matrix A, A−1/2
denotes the principal square root of its pseudo-inverse.
Proposition 10.1. Let vec be the vectorization operation of any matrix and ⊗ the Kronecker
product. Then vec(F k+1 ) is the solution to the following standard Lasso problem
min kΓx − bk2F +
x

2ρ
kxk1
η

b −1/2 Σ
b xy Σ
b −1/2
b y1/2 ⊗ Σ
b x1/2 and b = vec(Gk − 1 H k + 1 Σ
).
where Γ = Σ
y
η
η x

Remark 10.1. It is worth mentioning that the vectorized formulation in Proposition 10.1
is for illustration only. In practice, we solve the problem in (107) directly, since the vectorized version, especially the Kronecker product, would great increase the computation cost.
The solver to (107) can be easily implemented in standard software packages for convex
programming, such as TFOCS [5].
Since each update of F is the solution of some Lasso problem, it should be sparse in the
sense that its vector ℓ1 norm is well controlled.
49

Turning to the update for G, we note that (105) is equivalent to
η b 1/2 k+1 b 1/2
Gk+1 = argmin g(G) − hH k , Gi + kΣ
Σy − Gk2F
x F
2
G
1
η
k+1 b 1/2 2
b 1/2
Σy )kF
= argmin kG − ( H k + Σ
x F
2
η
G
+ ∞1{kGk∗ >r} + ∞1{kGkop >1}

1
b 1/2 F k+1 Σ
b 1/2 )k2
= argmin kG − ( H k + Σ
x
y
F
η
G

+ ∞1{kGk∗ >r} + ∞1{kGkop >1} .

(108)

The solution to the last display has a closed form according to the following result.
Proposition 10.2. Let G∗ be the solution to the optimization problem:
minimize kG − W kF

subject to kGk∗ ≤ r, kGkop ≤ 1.
P
′
Let the SVDPof W be W = m
i=1 ωi ai bi with ω1 ≥ · · · ≥ ωm ≥ 0 the ordered singular values.
m
Then G∗ = i=1 gi ai b′i where for any i, gi = 1 ∧ (ωi − γ ∗ )+ for some γ which is the solution
to
m
X
1 ∧ (ωi − γ)+ ≤ r.
minimize γ,
subject to γ > 0,
i=1

Proof. The proof essentially follows that of Lemma 4.1 in [37]. In addition to the fact that
the current problem deals
P with asymmetric matrix, the only difference that we now have an
inequality constraint i gi ≤ r rather than an equality constraint as in [37]. The asymmetry
of the current problem does not matter since it is orthogonally invariant.
Here and after, we call the operation in Proposition 10.2 singular value capped soft thresholding (SVCST) and write G∗ = SVCST(W ). Thus, any update for G results from the SVCST
operation of some matrix, and so it has well controlled singular values.
In summary, the convex program (18) is implemented as Algorithm 5.

11

Numerical Studies

This section presents numerical results demonstrating the competitive finite sample performance of the proposed adaptive estimation procedure CoLaR on simulated datasets.
Simulation settings We consider three simulation settings. In all these settings, we set
p = m, Σx = Σy = Σ, and r = 2 with λ1 = 0.9 and λ2 = 0.8. Moreover, the nonzero rows of
both U and V are set at {1, 6, 11, 16, 21}. The values at the nonzero coordinates are obtained
from normalizing (with respect to Σ) random numbers drawn from the uniform distribution
on the finite set {−2, 1, 0, 1, 2}. The choices of Σ in the three settings are as follows:
1. Identity: Σ = Ip .

50

Algorithm 5: An ADMM algorithm for SCCA
Input:
b x, Σ
b y and Σ
b xy ,
1. Sample covariance matrices Σ
2. Penalty parameter ρ,
3. Rank r,
4. ADMM parameter η and tolerance level ǫ.
1

2
3
4
5

6

b
Output: Estimated sparse canonical correlation signal A.
b xy ), G0 = 0, H 0 = 0.
Initialize: k = 0, F 0 = SVCST(Σ
repeat
Update F k+1 as in (104)
(Lasso) ;
k+1
−1
k+1 Σ
b 1/2
b 1/2
Update G
← SVCST(η H k + Σ
(SVCST) ;
x F
y )
1/2 k+1 b 1/2
k+1
k
k+1
b
Update H
← H + η(Σx F
Σy − G ) ;
k ←k+1 ;
until max{kF k+1 − F k kF , ρkGk+1 − Gk kF } ≤ ǫ;
b = F k.
Return A
2. Toeplitz: Σ = (σij ) where σij = 0.3|i−j| for all i, j ∈ [p]. In other words, Σx and Σy
are Toeplitz matrices.
q
0 / σ 0 σ 0 ). We set Σ0 = (σ 0 ) = Ω−1 where Ω = (ω ) with
3. SparseInv: Σ = (σij
ij
ij
ii jj
ωij = 1{i=j} + 0.5 × 1{|i−j|=1} + 0.4 × 1{|i−j|=2} ,

i, j ∈ [p].

In other words, Σx and Σy have sparse inverse matrices.
In all three settings, we normalize the variance of each coordinate to be one.
Implementation details The proposed CoLaR estimator in Section 4.1 has two stages.
The convex program (18) in the first stage can be solved via an ADMM algorithm [11]. The
details of the ADMM approach are presented in Section 10. The optimization problem (19)
in the second stage can be solved by a standard group-Lasso algorithm [45].
In
p all numerical results reported in this section, we used the same penalty level ρ =
0.55 log(p + m)/n in (18) and we used η = 2 in (106). Inp(19), we used five-fold cross
validation to select a common penalty parameter ρu = ρv = b (r + log p)/n. In particular,
test , Y test ) and the other four
for l = 1, . . . , 5, we use one fold of the data as the test set (X(l)
(l)
train , Y train ). For any choice of b, we solved (19) on (X train , Y train )
folds as the training set (X(l)
(l)
(l)
(l)
b(l) , Vb(l) ). Then we computed the sum of canonical correlations between
to obtain estimates (U
b(l) ∈ Rn×r and Y test Vb(l) ∈ Rn×r to obtain CVl (b). Finally, CV(b) = P5 CVl (b).
X test U
(l)

l=1

(l)

Among all the candidate penalty parameters, we select the b value such that CV(b) is maximized. The candidate b values used in the simulation below are {0.5, 1, 1.5, 2}. Throughout
the simulation, we used all the sample {(Xi , Yi )}ni=1 to form the sample covariance matrices
used in (18) – (20).

51

In addition to the performance of CoLaR, we also report that of the method proposed
in [43] (denoted by PMA here and on). The PMA seeks the solution to the optimization
problem
b xy v, subject to ||u|| ≤ 1, ||v|| ≤ 1, ||u||1 ≤ c1 , ||v||1 ≤ c2 .
max u′ Σ
u,v

The solution is used to estimate the first canonical pair (b
u1 , vb1 ). Then the same procedure is
′
′
b
b
b
repeated after Σxy is replaced by Σxy − (b
u1 Σxy vb1 )b
u1 vb1 , and the solution gives the estimator
of the second canonical pair (b
u2 , vb2 ). This process is repeated until u
br , vbr is obtained. Note
that the normalization constraint kuk ≤ 1 and kvk ≤ 1 implicitly assumes that the marginal
covariance matrices Σx and Σy are identity matrices. We used the R implementation of the
method (function CCA in the PMA package in R) by the authors of [43]. To remove undesired
amplification of error caused by normalization, we renormalized each individual u
bj with
b x and each individual vbj with respect to Σ
b y before calculating the error under
respect to Σ
the loss (7). For each simulated dataset, we set the sparsity penalty parameters penaltyx
and penaltyz of the function CCA at each of the eleven different values {0.6l : l = 0, 1, . . . , 10}
and only the smallest estimation error out of all eleven trials was used to compute the error
reported in the tables below.
Results Tables 1 – 3 report, in each of the three settings, the medians of the prediction
errors of CoLaR and PMA out of 100 repetitions for four different configurations of (p, m, n)
values.
In each table, the columns U -PMA and V -PMA report the medians of the smallest estimation errors out of the eleven trials on each simulated dataset. The columns U -init and
V -init report the median estimation errors of the renormalized r left singular vectors and
right singular vectors of the solutions to the initialization step (18), where the renormalization is the same as in (20) and in both (18) and renormalization we used all the n pairs
of observations. Last but not least, the columns U -CoLaR and V -ColaR report the median
estimation errors of the CoLaR estimators where both stages were carried out.
In all simulation settings, both the renormalized initial estimators and the CoLaR estimators consistently outperform PMA. Comparing the last four columns within each table,
we also find that the CoLaR estimators with both stages carried out significantly improve
over the renormalized initial estimators, which is in accordance with our theoretical results
in Section 4.
In summary, the proposed method delivers consistent and competitive performance in all
three covariance settings across all dimension and sample size configurations, and its behavior
agrees well with the theory.
(p, m, n)
(300, 300, 200)
(600, 600, 200)
(300, 300, 500)
(600, 600, 500)

U -PMA
2.1316
3.4154
0.2683
2.0335

V -PMA
2.1297
3.3584
0.2701
2.0368

U -init
0.2653
0.3167
0.1207
0.1448

V -init
0.1712
0.2087
0.0665
0.0817

U -CoLaR
0.0498
0.0671
0.0135
0.0166

V -CoLaR
0.0646
0.0776
0.0159
0.0203

Table 1: Prediction errors (Identity): Median in 100 repetitions.

52

(p, m, n)
(300, 300, 200)
(600, 600, 200)
(300, 300, 500)
(600, 600, 500)

U -PMA
2.1853
3.4247
0.2358
2.1214

V -PMA
2.1840
3.4852
0.2191
2.0889

U -init
0.2885
0.3236
0.1202
0.1408

V -init
0.1706
0.2004
0.0664
0.0811

U -CoLaR
0.0511
0.0638
0.0135
0.0176

V -CoLaR
0.0601
0.0764
0.0166
0.0209

Table 2: Prediction errors (Toeplitz): Median in 100 repetitions.
(p, m, n)
(300, 300, 200)
(600, 600, 200)
(300, 300, 500)
(600, 600, 500)

U -PMA
2.9697
4.6908
2.3967
2.8707

V -PMA
2.9619
4.3339
2.0620
2.8609

U -init
0.5552
0.5596
0.2695
0.3068

V -init
0.5718
0.6133
0.1917
0.2368

U -CoLaR
0.1568
0.2123
0.0242
0.0338

V -CoLaR
0.1194
0.1572
0.0219
0.0271

Table 3: Prediction errors (SparseInv): Median in 100 repetitions.
Model misspecification We now examine the performance of our estimator when the
model is misspecified. To this end, we consider the case where there are three pairs of nontrivial canonical correlations present in the data but we set r = 2 in our algorithm. As before,
we consider three different types of marginal covariance matrices: Identity, Toeplitz and
SparseInv. In addition, we generate the first two pairs of canonical correlation vectors
in the same way as before. For the generation of the third pair of canonical directions,
we consider two different scenarios. In the first scenario, the support of the third pair of
canonical direction vectors are set at {1, 6, 11, 16, 21} and so they are the same as those of
the first two pairs. In the second scenario, we put no constraint on the support of these
vectors. For both scenarios, we set λ3 = 0.3. Table 4 reports the prediction errors of the first
two pairs of canonical correlations in both scenarios when (p, m, n) = (300, 300, 500). The
implementation details are exactly the same as before. The first two columns contain results
in the first scenario, and the third and the fourth columns the second scenario. Comparing
these results with their counterparts in correctly specified models (the last two cells in the
second last rows of Tables 1–3), we have found that the performance of our estimator was
robust to model misspecification in both scenarios.
Covariance
Identity
Toeplitz
SparseInv

Scenario 1
U -CoLaR V -CoLaR
0.0190
0.0195
0.0197
0.0197
0.0348
0.0263

Scenario 2
U -CoLaR V -CoLaR
0.0144
0.0160
0.0138
0.0147
0.0221
0.0328

Table 4: Prediction errors with misspecified models: Median in 100 repetitions. (p, m, n) =
(300, 300, 500).

53

