Sparse CCA: Adaptive Estimation and Computational Barriers

Introduction

where $\Sigma_{x}=\text{Cov}(X),\Sigma_{y}=\text{Cov}(Y),\Sigma_{xy}=\text{Cov}(X,Y)$, ${u}_{0}=0$, and ${v}_{0}=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}^{\prime}X,v_{j}^{\prime}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 population covariance, and $\Sigma_{x}$, $\Sigma_{y}$, and $\Sigma_{xy}$ are replaced by their sample versions $\widehat{\Sigma}_{x}$, $\widehat{\Sigma}_{y}$, and $\widehat{\Sigma}_{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 . 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, .

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 . 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

Here $U=[u_{1},...,u_{r}]$ and $V=[v_{1},...,v_{r}]$ collect the canonical coefficient vectors and $\Lambda=\mathop{\text{diag}}(\lambda_{1},\dots,\lambda_{r})$ with $1>\lambda_{1}\geq\cdots\geq\lambda_{r}>0$ are the ordered canonical correlations. Let $S_{u}={\rm supp}(U)$ and $S_{v}={\rm 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 $S_{u}$ and $S_{v}$ to be small, namely $|S_{u}|\leq s_{u}$ and $|S_{v}|\leq s_{v}$ for some $s_{u}\leq p$ and $s_{v}\leq m$. Under this model, Gao et al. showed that the minimax rate for estimating $U$ and $V$ under the joint loss function $\|\widehat{U}\widehat{V}^{\prime}-UV^{\prime}\|_{\rm F}^{2}$ is

However, to achieve the rate, Gao et al. used a computationally infeasible and nonadaptive procedure, which requires exhaustive search of all possible subsets with the given cardinality and the knowledge of $s_{u}$ and $s_{v}$. 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?

We now introduce the main contributions of the present paper from three different viewpoints as suggested by the three questions we have raised.

The joint loss $\|\widehat{U}\widehat{V}^{\prime}-UV^{\prime}\|_{\rm F}^{2}$ studied by characterizes the 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,\lambda_{r},p$ and $s_{u}$, but not on either $m$ or $s_{v}$. 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$.

As pointed out in and , 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 . The additional difficulty in sparse CCA mainly comes from the presence of the nuisance parameters $\Sigma_{x}$ and $\Sigma_{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 by considering the smallest convex set containing all matrices of the form $AB^{\prime}$ with both $A$ and $B$ being rank-$r$ orthogonal matrices. The convex program can be efficiently solved by the Alternating Direction Method with Multipliers (ADMM) . Based on the output of the first stage, we formulate a sparse linear regression problem in the second stage to improve estimation accuracy, and the final estimator $\widehat{U}$ and $\widehat{V}$ can be obtained via a group-Lasso algorithm . Under the sample size condition that

for some sufficiently large constant $C>0$, we show $\widehat{U}$ and $\widehat{V}$ recover the true canonical coefficient matrices $U$ and $V$ within optimal error rates adaptively with high probability.

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 , namely $n\geq C{s^{2}\log p}/{\lambda^{2}}$ where $s$ is the sparsity of the leading eigenvector and $\lambda$ the gap between the leading eigenvalue and the rest of the spectrum. Berthet and Rigollet 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 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 . The same remark also applies to the subsequent work on sparse PCA estimation . 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.

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 .

3 Notation

Problem Formulation

Consider a canonical pair model where the observed pairs of measurement vectors $(X_{i}^{\prime},Y_{i}^{\prime})^{\prime}$, $i=1,\dots,n$ are i.i.d. from a multivariate Gaussian distribution $N_{p+m}(0,\Sigma)$ where

with the cross-covariance matrix $\Sigma_{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 and sparse CCA problems when one seeks multiple sparse vectors simultaneously.

Recall that for any matrix $A$, ${\rm supp}(A)$ collects the indices of nonzero rows in $A$. Adopting the above notion of sparsity, we define $\mathcal{F}(s_{u},s_{v},p,m,r,\lambda;M)$ to be the collection of all covariance matrices $\Sigma$ with the structure (2) satisfying

where $n$ is the sample size. We shall allow $s_{u},s_{v},p,m,r,\lambda$ to vary with $n$, while $M>1$ is restricted to be an absolute constant.

2 Prediction loss

By symmetry, we can define $L(\widehat{V},V)$ by simply replacing $U$, $\widehat{U}$, $X^{\star}$ and $\Sigma_{x}$ in (7) and (8) with $V$, $\widehat{V}$, $Y^{\star}$ and $\Sigma_{y}$.

A related loss function is $\|P_{\widehat{U}}-P_{U}\|_{\rm F}^{2}$ measuring the difference between two subspaces. By Proposition 9.2 in the supplementary material , the prediction loss $L(\widehat{U},U)$ is a stronger loss function. That is, $\|P_{\widehat{U}}-P_{U}\|_{\rm F}^{2}\leq CL(\widehat{U},U)$ for some constant $C>0$ only depending on $M$. Actually, $L(\widehat{U},U)$ is strictly stronger. To see this, let $\Sigma_{x}=I_{p}$, $U\in O(p,r)$ and $\widehat{U}=2U$. Then, $\|P_{\widehat{U}}-P_{U}\|_{\rm F}^{2}=0$, while $L(\widehat{U},U)=\inf_{W\in O(r)}\|\Sigma_{x}^{1/2}(\widehat{U}W-U)\|_{\rm F}^{2}=\inf_{W\in O(r)}(5r-\mathop{\sf Tr}(W))=r>0$. In this paper, we will focus on the stronger loss $L(\widehat{U},U)$, and provide brief remarks on results for $\|P_{\widehat{U}}-P_{U}\|_{\rm F}^{2}$.

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.

To obtain minimax upper bound, we propose a two-stage combinatorial optimization procedure. We split the data into three equal size batches $\mathcal{D}_{0}=\{(X_{i}^{\prime},Y_{i}^{\prime})^{\prime}\}_{i=1}^{n_{0}}$, $\mathcal{D}_{1}=\{(X_{i}^{\prime},Y_{i}^{\prime})^{\prime}\}_{i=n_{0}+1}^{2n_{0}}$ and $\mathcal{D}_{2}=\{(X_{i}^{\prime},Y_{i}^{\prime})^{\prime}\}_{i=2n_{0}+1}^{n}$, and denote the sample covariance matrices computed on each batch by $\widehat{\Sigma}_{x}^{(j)},\widehat{\Sigma}_{y}^{(j)}$ and $\widehat{\Sigma}_{xy}^{(j)}$ for $j\in\{0,1,2\}$.

In the first stage, we find $(\widehat{U}^{(0)},\widehat{V}^{(0)})$ which solves the following program:

In the second stage, we further refine the estimator for $U$ by finding $\widehat{U}^{(1)}$ solving

The final estimator is a normalized version of $\widehat{U}^{(1)}$, defined as

The purpose of sample splitting employed in the above procedure is to facilitate the proof.

where the expectation is with respect to the distribution $(X^{\prime},Y^{\prime})^{\prime}\sim N_{p+m}(0,\Sigma)$. 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 $\mathop{\sf Tr}(V^{\prime}\Sigma_{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^{\prime}Y=\Lambda U^{\prime}X+E$, where $E\sim N(0,I_{r}-\Lambda^{2})$. Then, (10) is a least square formulation of the regression interpretation. However, CCA is different from regression because the response $V^{\prime}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 $\Sigma$ are replaced by the initial estimator $\widehat{V}^{(0)}$ and sample covariance matrix from an independent sample. Therefore, $\widehat{U}^{(1)}$ can be viewed as an estimator of $U\Lambda$. Hence, 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).

for some sufficiently small constant $c>0$. Then there exist constants $C,C^{\prime}>0$ only depending on $c$ such that

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 $\Sigma_{x}$ are bounded in the interval $[M_{1},M_{2}]$. The convergence rate of $L(\widehat{U},U)$ would still be \frac{1}{n\lambda^{2}}s_{u}\big{(}r+\log\frac{ep}{s_{u}}\big{)}, because the dependence on $M_{1},M_{2}$ has been implicitly built into the prediction loss. On the other hand, a convergence rate for the loss $\|P_{\widehat{U}}-P_{U}\|_{\rm F}^{2}$ would be \big{(}\frac{M_{2}}{M_{1}}\big{)}\frac{1}{n\lambda^{2}}s_{u}\big{(}r+\log\frac{ep}{s_{u}}\big{)}, with an extra factor of the condition number of $\Sigma_{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 $s_{v}$. However, it can be shown that $\exp\left(-C^{\prime}(s_{v}+\log(em/s_{v}))\right)\leq m^{-C^{\prime}/2}$, and so the corresponding term in the tail probability goes to as long as $m\rightarrow\infty$. The optimality of this upper bound can be justified by the following minimax lower bound.

Assume that $r\leq\frac{s_{u}\wedge s_{v}}{2}$. Then there exists some constant $C>0$ only depending on $M$ and an absolute constant $c_{0}>0$, such that

where $\mathcal{P}=\mathcal{P}(n,s_{u},s_{v},p,m,r,\lambda;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.

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 $s_{u}$ and $s_{v}$. 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.

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 split into three batches $\mathcal{D}_{0},\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ of equal sizes and for $j=0,1,2$, let $\widehat{\Sigma}_{x}^{(j)},\widehat{\Sigma}_{y}^{(j)}$ and $\widehat{\Sigma}_{xy}^{(j)}$ be defined as before.

Naturally, we relax it to $(\widehat{\Sigma}_{x}^{(0)})^{1/2}F(\widehat{\Sigma}_{y}^{(0)})^{1/2}\in\mathcal{C}_{r}$ where

is the smallest convex set containing $\mathcal{O}_{r}$. The relation (17) is stated in the proof of Theorem 3 of . Combining (15) – (17), we use the following convex program for the first stage in our adaptive estimation scheme:

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

A related but different convex relaxation was proposed in 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 $\left\{P:\mathop{\sf Tr}(P)=r,0\preceq P\preceq I_{p}\right\}$. Such an idea is not directly applicable in the current setting due to the asymmetric nature of the matrices included in the set $\mathcal{O}_{r}$ in (16).

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

2 Theoretical guarantees

We first state the upper bound for the solution $\widehat{A}$ to the convex program (18).

for some sufficiently large constant $C_{1}>0$. Then there exist positive constants $\gamma_{1},\gamma_{2}$ and $C,C^{\prime}$ only depending on $M$ and $C_{1}$, such that when $\rho=\gamma\sqrt{{\log(p+m)}/{n}}$ for $\gamma\in[\gamma_{1},\gamma_{2}]$,

Note that the error bound in Theorem 4.1 can be much larger than the optimal rate for joint estimation of $UV^{\prime}$ established in . Nonetheless, under the sample size condition (21), it still ensures that $\widehat{A}$ is close to $UV^{\prime}$ in Frobenius norm distance. This fact, together with the proposed refinement scheme (19) – (20), guarantees the optimal rates of convergence for the estimator (20) as stated in the following theorem.

Assume (21) holds for some sufficiently large $C_{1}\geq 0$. Then there exist constants $\gamma$ and $\gamma_{u}$ only depending on $C_{1}$ and $M$ such that if we set $\rho=\gamma^{\prime}\sqrt{{[\log(p+m)]}/{n}}$ and $\rho_{u}=\gamma^{\prime}_{u}\sqrt{({r+\log p})/{n}}$ for any $\gamma^{\prime}\in[\gamma,C_{2}\gamma]$ and $\gamma_{u}^{\prime}\in[\gamma_{u},C_{2}\gamma_{u}]$ for some absolute constant $C_{2}>0$, there exist a constants $C,C^{\prime}>0$ only depending on $C_{1},C_{2}$ and $M$, such that

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 the eigenvalues of $\Sigma_{x}$ all lie in the interval $[M_{1},M_{2}]$, then the convergence rate of $L(\widehat{U},U)$ would be \big{(}\frac{M_{2}}{M_{1}}\big{)}^{2}\frac{s_{u}\left(r+\log p\right)}{n\lambda^{2}} and a convergence rate of $\|P_{\widehat{U}}-P_{U}\|_{\rm F}^{2}$ would be \big{(}\frac{M_{2}}{M_{1}}\big{)}^{3}\frac{s_{u}\left(r+\log p\right)}{n\lambda^{2}}. Compared with Remark 3.1, there is an extra factor \big{(}\frac{M_{2}}{M_{1}}\big{)}^{2}, which is also present for the Lasso error bounds . Evidence has been given in the literature that such an extra factor can be intrinsic to all polynomial-time algorithms .

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 $\rho$ and $\rho_{u}$ in (18) and (19) do not depend on $s_{u}$ or $s_{v}$. 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 . A careful examination of the proofs shows that the dependence of $\rho$ and $\rho_{u}$ on $M$ is through $\|\Sigma_{x}\|_{\rm op}^{1/2}\|\Sigma_{y}\|_{\rm op}^{1/2}$ and $\|\Sigma_{x}\|_{\rm op}$, 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.

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.

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) – (LABEL:eq:para-space), 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.

Let $N$ be a positive integer and $k\in[N]$. We denote by $\mathcal{G}(N,1/2)$ the Erdős-Rényi graph on $N$ vertices where each edge is drawn independently with probability $1/2$, and by $\mathcal{G}(N,1/2,k)$ the random graph generated by first sampling from $\mathcal{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\in\{0,1\}^{N\times N}$ of an instance from either $\mathcal{G}(N,1/2)$ or $\mathcal{G}(N,1/2,k)$, the Planted Clique detection problem of parameter $(N,k)$ refers to testing the following hypotheses

It is widely believed that when $k=O(N^{1/2-\delta})$, 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.

For any sequence $k=k(N)$ such that $\limsup_{N\to\infty}\frac{\log k}{\log N}<\frac{1}{2}$ and any randomized polynomial-time test $\psi$,

Evidence supporting this hypothesis has been provided in . Computational lower bounds in several statistical problems have been established by assuming the above hypothesis and its close variants, including sparse PCA detection and estimation in classes defined by a restricted covariance concentration condition, submatrix detection and community detection .

Under Hypothesis A, the necessity of condition (21) is supported by the following theorem.

Suppose that Hypothesis A holds and that as $n\to\infty$, $p=m$ satisfying $2n\leq p\leq n^{a}$ for some constant $a>1$, $s_{u}=s_{v}$, $n(\log n)^{5}\leq cs_{u}^{4}$ for some sufficiently small $c>0$, and $\lambda=\frac{s_{u}s_{v}}{7290n(\log(12n))^{2}}$. If for some $\delta\in(0,1)$,

then for any randomized polynomial-time estimator $\widehat{u}$,

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.

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 $\Theta(1)$ . To be rigorous under the probabilistic Turing machine model , 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 .

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.

1 Hardness of sparse PCA under Gaussian spiked covariance model

Gaussian single spiked model refers to the distribution $N_{p}(0,\Sigma)$ where $\Sigma=\tau\theta\theta^{\prime}+I_{p}$. Here, $\theta$ is the eigenvector of unit length and $\tau>0$ is the eigenvalue. Define the following Gaussian single spiked model parameter space for sparse PCA

The minimax estimation rate for $\theta$ under the loss $\|P_{\widehat{\theta}}-P_{\theta}\|_{\rm F}^{2}$ is $\frac{\lambda+1}{n\lambda^{2}}s\log\frac{ep}{s}$. See, for instance, . However, to achieve the above minimax rate via computationally efficient methods such as those proposed in , researchers have required the sample size to satisfy $n\geq C\frac{s^{2}\log p}{\lambda^{2}}$ for some sufficiently large constant $C>0$. Moreover, no computationally efficient 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 of sparse PCA under Gaussian spiked covariance model (25) when $\liminf_{n\to\infty}\frac{s^{2-\delta}\log p}{n\lambda^{2}}>0$ for some $\delta>0$.

We note that previous computational lower bounds for sparse PCA in 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 , 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 $\Sigma=\tau\theta\theta^{\prime}+I_{p}$. 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.

Suppose we have an estimator $\widehat{\theta}=\widehat{\theta}(W_{1},\dots,W_{n})$ of the leading sparse eigenvector, we 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

denote the density function of the Gaussian mixture $\frac{1}{2}N(\mu,1)+\frac{1}{2}N(-\mu,1)$. Next, let $\widetilde{\Phi}_{0}$ be the restriction of the $N(0,1)$ distribution on the interval $[-3\sqrt{\log N},3\sqrt{\log N}]$. For any $|\mu|\leq 3\sqrt{\eta_{N}\log N}$, define two probability distributions $\mathcal{F}_{\mu,0}$ and $\mathcal{F}_{\mu,1}$ with densities

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\geq 12n$ and $p\geq 2n$.

We now explain how the reduction achieves its goal. For simplicity, focus on the case where $p=2n$. Let ${\epsilon}=({\epsilon}_{1},...,{\epsilon}_{2n})\in\{0,1\}^{2n}$ where $\epsilon_{i}$ is the indicator of whether the $i$-th row of $A_{0}$ (defined in Step 2 of Algorithm 1) belongs to the planted clique or not, and ${\gamma}=({\gamma}_{1},...,{\gamma}_{2n})$ the indicators of the columns of $A_{0}$. In what follows, we discuss the distributions of $W$ when $A\sim H_{0}^{G}$ and $H_{1}^{G}$, respectively.

When $A\sim H_{0}^{G}$, the $\epsilon_{i}$’s and $\gamma_{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 $W_{ij}$ is close to the $N(0,1)$ distribution (c.f., Lemma 7.1 in the supplement ). Hence, the rows of $W$ are close to i.i.d. random vectors from the $N_{p}(0,I_{p})$ distribution. Since $\widehat{\theta}=\widehat{\theta}(W_{1},\dots,W_{n})$ is independent of $\{W_{i}\}_{i=n+1}^{2n}$, the LHS of (40) is close in distribution to a $\chi^{2}_{n}$ random variable scaled by $n$ which concentrates around its expected value one. Indeed, it is upper bounded by $1+O(\sqrt{\log(n)/n})$ with high probability.

If $A\sim H_{1}^{G}$, then the $(i,j)$-th entry of $A_{0}$ is an edge in the planted clique if and only if $\epsilon_{i}=\gamma_{j}=1$. Moreover, the joint distribution of $\{\epsilon_{1},\dots,\epsilon_{2n},\gamma_{1},\dots,\gamma_{2n}\}$ is close to that of $4n$ i.i.d. Bernoulli random variables $\{\widetilde{\epsilon}_{1},\dots,\widetilde{\epsilon}_{2n},\widetilde{\gamma}_{1},\dots,\widetilde{\gamma}_{2n}\}$ with success probability $\delta_{N}=k/N$. For simplicity, suppose that these indicators are indeed i.i.d. Bernoulli($\delta_{N}$) variables $\{\widetilde{\epsilon}_{1},\dots,\widetilde{\epsilon}_{2n},\widetilde{\gamma}_{1},\dots,\widetilde{\gamma}_{2n}\}$. Then, one can show that conditioning on $\widetilde{\gamma}_{j}=0$, for any $i\in[2n]$, the conditional distribution of $(W_{ij}|\widetilde{\gamma}_{j}=0)$, after integrating over the conditional distribution of $\widetilde{\epsilon}_{i}$, $\mu_{i}$ and $(A_{0})_{ij}$, is approximately $N(0,1)$. In contrast, conditioning on $\widetilde{\gamma}_{j}=1$, for any $i\in[2n]$, the conditional distribution of $(W_{ij}|\widetilde{\gamma}_{j}=1)$ is approximately $N(0,1+\eta_{N})$. Therefore, conditioning on $\widetilde{\gamma}$ the distribution of the $W_{i}$’s is close to that of $2n$ i.i.d. random vectors sampled from

i.e., a Gaussian spiked covariance model in (25). Here, the leading eigenvector $\theta$ has sparsity level $|{\rm supp}(\theta)|=|{\rm supp}(\widetilde{\gamma})|=\sum_{j}\widetilde{\gamma}_{j}$, which concentrates around its mean value $n\delta_{N}\asymp k$ if $N\asymp n$. Thus, if $\widehat{\theta}$ estimates $\theta$ well, then the LHS of (40) approximately follows a non-central $\chi_{n}$ distribution scaled by $n$, which should exceed $1+O(\sqrt{\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 $\widehat{\theta}$ 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).

For some sufficiently small constant $c>0$, assume $\frac{N(\log N)^{5}}{k^{4}}\leq c$, $cN\leq n\leq N/12$ and $p\geq 2n$. Then, for any $\widehat{\theta}$ such that

the test $\psi$ defined by (30) – (33) satisfies

for sufficiently large $n$ with some constants $C,C^{\prime}>0$.

If the estimator $\widehat{\theta}$ is uniformly consistent over $\mathcal{Q}(n,3k/2,p,k\eta_{N}/2)$, then $\beta$ 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.

Suppose that Hypothesis A holds and that as $n\rightarrow\infty$, $2n\leq p\leq n^{a}$ for some constant $a>1$, $n(\log n)^{5}\leq cs^{4}$ for some sufficiently small $c>0$, and $\lambda=\frac{s^{2}}{2430n(\log(12n))^{2}}$. If for some $\delta\in(0,2)$,

then for any randomized polynomial-time estimator $\widehat{\theta}$,

Under the same condition of Theorem 5.3, for any randomized polynomial-time test $\phi$ for testing (37),

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).

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.

To see why Algorithm 2 is effective, one can verify that if $W_{i}\stackrel{{\scriptstyle iid}}{{\sim}}N_{p}(0,\tau\theta\theta^{\prime}+I_{p})$, then $(X_{i}^{\prime},Y_{i}^{\prime})^{\prime}\stackrel{{\scriptstyle iid}}{{\sim}}N_{p+m}(0,\Sigma)$ where

with $u=v=\frac{\theta}{\sqrt{\tau/2+1}},\,\lambda=\frac{\tau/2}{\tau/2+1}$. This is a special case of the Gaussian canonical pair model (2). Thus, the leading eigenvector of $W_{i}$ aligns with the leading canonical coefficient vectors of $(X_{i},Y_{i})$. Exploiting this connection, we obtain the following theorem.

Consider $p=m$, $s_{u}=s_{v}$ and $\lambda\leq 1$. Then for any $\widehat{u}$ such that

the estimator $\widehat{\theta}$ defined by Algorithm 2 satisfies

If we start with an estimator $\widehat{u}$ 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 $\widehat{\theta}$ from $\widehat{u}$ in the third step. Finally, the desired Theorem 5.1 is a direct consequence of Theorems 5.3 and 5.5.

Proofs

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

Before presenting the proof, we state some technical lemmas. The proofs of all the lemmas are given in Section 9.3 in the supplement . First, note that the estimator is normalized with respect to $\widehat{\Sigma}^{(0)}_{x}$ and $\widehat{\Sigma}^{(0)}_{y}$, while the truth $U$ and $V$ is normalized with respect to $\Sigma_{x}$ and $\Sigma_{y}$. To address this issue, we normalize the truth with respect to $\widehat{\Sigma}_{x}^{(0)}$ and $\widehat{\Sigma}_{y}^{(0)}$ to obtain $\widetilde{U}=U(U^{\prime}\widehat{\Sigma}_{x}^{(0)}U)^{-1/2}$ and $\widetilde{V}=V(V^{\prime}\widehat{\Sigma}_{y}^{(0)}V)^{-1/2}$. Also define $\widetilde{\Lambda}=(U^{\prime}\widehat{\Sigma}_{x}^{(0)}U)^{1/2}\Lambda(V^{\prime}\widehat{\Sigma}_{y}^{(0)}V)^{1/2}$. For notational convenience, define

The following lemma bounds the normalization effect.

Assume $\epsilon_{n,u}^{2}+\epsilon_{n,v}^{2}\leq c$ for some sufficiently small constant $c\in(0,1)$. Then there exist some constants $C,C^{\prime}>0$ only depending on $c$ such that

with probability at least $1-\exp\left(-C^{\prime}(s_{u}+\log(ep/s_{u}))\right)-\exp\left(-C^{\prime}(s_{v}+\log(em/s_{v}))\right)$.

Using the definitions of $\widetilde{U}$ and $\widetilde{V}$, let us state the following lemma, which asserts that the matrix $\widetilde{A}=\widetilde{U}\widetilde{V}^{\prime}$ is feasible to the optimization problem (18).

Define $\widetilde{A}=\widetilde{U}\widetilde{V}^{\prime}$. When $\widetilde{A}$ exists, we have

As was argued in Section 4.1, the set $\mathcal{C}_{r}$ is the convex hull of $\mathcal{O}_{r}$. The following curvature lemma shows that the relaxation $\mathcal{C}_{r}$ preserves the restricted strong convexity of the objective function.

Lemma 6.4 is instrumental in determining the proper value of the tuning parameter required in the program (18).

Assume $r\sqrt{{[\log(p+m)]}/{n}}\leq c$ for some sufficiently small constant $c\in(0,1)$. Then there exist some constants $C,C^{\prime}>0$ only depending on $M$ and $c$ such that $||\widehat{\Sigma}_{xy}^{(0)}-\widetilde{\Sigma}_{xy}||_{\infty}\leq C\sqrt{{[\log(p+m)]}/{n}}$, with probability at least $1-(p+m)^{-C^{\prime}}$.

We also need a lemma on restricted eigenvalue. For any p.s.d. matrix $B$, define

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

Assume \frac{1}{n}\big{(}(k_{u}\wedge p)\log(ep/(k_{u}\wedge p))+(k_{v}\wedge m)\log(em/(k_{v}\wedge m))\big{)}\leq c for some sufficiently small constant $c>0$. Then there exist some constants $C,C^{\prime}>0$ only depending on $M$ and $c$ such that for $\delta_{u}(k_{u})=\sqrt{\frac{(k_{u}\wedge p)\log(ep/(k_{u}\wedge p))}{n}}$ and $\delta_{v}(k_{v})=\sqrt{\frac{(k_{v}\wedge m)\log(em/(k_{v}\wedge m))}{n}}$, we have

with probability at least 1-\exp\big{(}-C^{\prime}(k_{u}\wedge p)\log(ep/(k_{u}\wedge p))\big{)}-\exp\big{(}-C^{\prime}(k_{v}\wedge m)\log(em/(k_{v}\wedge m))\big{)}, for $j=0,1,2$.

Finally, we need a result on subspace distance. Recall that for a matrix $F$, $P_{F}$ denotes the projection matrix onto its column subspace.

If further $G\in O(d,r)$, then $\inf_{W\in O(r,r)}\|F-GW\|_{\rm F}^{2}=\frac{1}{2}\|P_{F}-P_{G}\|_{\rm F}^{2}.$

Proofs of Lemma 6.1-6.6 are given in Section 9.3.1 of the supplement .

In the rest of this proof, we denote $\widehat{\Sigma}_{x}^{(0)}$, $\widehat{\Sigma}_{y}^{(0)}$ and $\widehat{\Sigma}_{xy}^{(0)}$ by $\widehat{\Sigma}_{x}$, $\widehat{\Sigma}_{y}$ and $\widehat{\Sigma}_{xy}$ for notational convenience. We also let $\Delta=\widehat{A}-\widetilde{A}$. The proof consists of two steps. In the first step, we are going to derive an upper bound for $\|\widehat{\Sigma}_{x}^{1/2}\Delta\widehat{\Sigma}_{y}^{1/2}\|_{\rm F}$. In the second step, we derive a generalized cone condition and use it to lower bound $\|\widehat{\Sigma}_{x}^{1/2}\Delta\widehat{\Sigma}_{y}^{1/2}\|_{\rm F}$ by a constant multiple of $\|\Delta\|_{\rm F}$ and hence the upper bound on $\|\widehat{\Sigma}_{x}^{1/2}\Delta\widehat{\Sigma}_{y}^{1/2}\|_{\rm F}$ leads to an upper bound on $\|\Delta\|_{\rm F}$.

Step 1. By Lemma 6.1, $\widetilde{U}$ and $\widetilde{V}$ are well-defined with high probability. Thus, $\widetilde{A}$ is well-defined with high probability, and we have

with probability at least $1-\exp\left(-C^{\prime}(s_{u}+\log(ep/s_{u}))\right)-\exp\left(-C^{\prime}(s_{v}+\log(em/s_{v}))\right)$. According to Lemma 6.2, $\widetilde{A}$ is feasible. Then, by the definition of $\widehat{A}$, we have

where $\widetilde{\Sigma}_{xy}$ is defined in (45). For the first term on the right hand side of (47), we have

For the second term on the right hand side of (47), we have $\langle\widehat{\Sigma}_{xy}-\widetilde{\Sigma}_{xy},\Delta\rangle\leq||\widehat{\Sigma}_{xy}-\widetilde{\Sigma}_{xy}||_{\infty}||\Delta||_{1}$. Thus when

Using Lemma 6.3, we can lower bound the left hand side of (49) as

where $\delta=\|\widetilde{\Lambda}-\Lambda\|_{\rm F}$. Combining (49) and (50), we have

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

which gives rise to the generalized cone condition that we are going to use in Step 2. Finally, by the bound $||\Delta_{S_{u}S_{v}}||_{1}\leq\sqrt{s_{u}s_{v}}\rho\|\Delta_{S_{u}S_{v}}\|_{\rm F}$ and (53), we have

Step 2. By (54), we have obtained the following condition

Due to the existence of the extra term $5\delta^{2}/(\rho\lambda_{r})$ on the RHS, we call it a generalized cone condition. In this step, we are going to lower bound $\|\widehat{\Sigma}_{x}^{1/2}\Delta\widehat{\Sigma}_{y}^{1/2}\|_{\rm F}$ by $\|\Delta\|_{\rm F}$ on the generalized cone. Motivated by the argument in , let the index set $J_{1}=\{(i_{k},j_{k})\}_{k=1}^{t}$ in $(S_{u}\times S_{v})^{c}$ correspond to the entries with the largest absolute values in $\Delta$, and we define the set $\widetilde{J}=(S_{u}\times S_{v})\cup J_{1}$. Now we partition $\widetilde{J}^{c}$ into disjoint subsets $J_{2},...,J_{K}$ of size $t$ (with $|J_{K}|\leq t$), such that $J_{k}$ is the set of (double) indices corresponding to the entries of $t$ largest absolute values in $\Delta$ outside $\widetilde{J}\cup\bigcup_{j=2}^{k-1}J_{j}$. By triangle inequality,

where we have used the generalized cone condition (56). Hence, we have the lower bound

Taking $t=c_{1}s_{u}s_{v}$ for some sufficiently large constant $c_{1}>1$, with high probability, $\kappa_{1}$ can be lower bounded by a positive constant $\kappa_{0}$ only depending on $M$. To see this, note that by Lemma 6.5, (58) can be lower bounded by the difference of $\sqrt{M^{-1}-C\delta_{u}(2c_{1}s_{u}s_{v})}\sqrt{M^{-1}-C\delta_{v}(2c_{1}s_{u}s_{v})}$ and $9c_{1}^{-1/2}\sqrt{M+C\delta_{u}(c_{1}s_{u}s_{v})}\sqrt{M+C\delta_{v}(c_{1}s_{u}s_{v})}$, where $\delta_{u}$ and $\delta_{v}$ are defined as in Lemma 6.5. It is sufficient to show that $\delta_{u}(2c_{1}s_{u}s_{v})$, $\delta_{v}(2c_{1}s_{u}s_{v})$, $\delta_{u}(c_{1}s_{u}s_{v})$ and $\delta_{v}(c_{1}s_{u}s_{v})$ are sufficiently small to get a positive absolute constant $\kappa_{0}$. For the first term, when $2c_{1}s_{u}s_{v}\leq p$, it is bounded by $\frac{2c_{1}s_{u}s_{v}\log(ep)}{n}$ and is sufficiently small under the assumption (13). When $2c_{1}s_{u}s_{v}>p$, it is bounded by $\frac{2c_{1}s_{u}s_{v}}{n}$ and is also sufficiently small under (13). The same argument also holds for the other terms. Similarly, $\kappa_{2}$ can be upper bounded by some constant.

Together with (55), this brings the inequality

Summing (59) and (60), we obtain a bound for $\|\Delta\|_{\rm F}$. According to Lemma 6.4, we may choose $\rho=\gamma\sqrt{[{\log(p+m)}]/{n}}$ for some large $\gamma$, so that (48) holds with high probability. By Lemma 6.1, $\delta\leq C\sqrt{{r(s_{u}+s_{v}+\log(p+m))}/{n}}\leq C^{\prime}\rho\sqrt{t}$ with high probability. Hence,

with high probability. This completes the second step. Finally, the triangle inequality leads to $\|\widehat{A}-UV^{\prime}\|_{\rm F}\leq\|\Delta\|_{\rm F}+\|\widetilde{A}-UV^{\prime}\|_{\rm F}$. By (46) and (61), the proof is complete. ∎

2 Proof of Theorem 4.2

Define $U^{*}=U\Lambda V^{\prime}\Sigma_{y}\widehat{V}^{(0)}$ and $\Delta=\widehat{U}^{(1)}-U^{*}$.

Assume $\frac{r+\log p}{n}\leq c$ for some sufficiently small constant $c\in(0,1)$. Then there exist some constants $C,C^{\prime}>0$ only depending on $M$ and $c$ such that $\max_{1\leq j\leq p}||[\widehat{\Sigma}^{(1)}_{xy}\widehat{V}^{(0)}-\widehat{\Sigma}^{(1)}_{x}U^{*}]_{j\cdot}||\leq C\sqrt{{(r+\log p)}/{n}}$, with probability at least 1-\exp\big{(}-C^{\prime}(r+\log p)\big{)}.

The proof of Lemma 6.7 is given in Section 9.3.1 of the supplement .

In the rest of this proof, we denote $\widehat{\Sigma}_{x}^{(1)}$, $\widehat{\Sigma}_{y}^{(1)}$ and $\widehat{\Sigma}_{xy}^{(1)}$ by $\widehat{\Sigma}_{x}$, $\widehat{\Sigma}_{y}$ and $\widehat{\Sigma}_{xy}$ for simplicity of notation. Note that they depends on $\mathcal{D}_{1}$, while the estimator $\widehat{V}^{(0)}$ depends on $\mathcal{D}_{0}$. Hence, $\widehat{V}^{(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 $\mathop{\sf Tr}(\Delta^{\prime}\widehat{\Sigma}_{x}\Delta)$. In the second step, we derive a cone condition and use it to obtain a bound for $\|\Delta\|_{\rm F}$ by arguing that $\mathop{\sf Tr}(\Delta^{\prime}\widehat{\Sigma}_{x}\Delta)$ upper bounds $\|\Delta\|_{\rm F}$. In the last step, we derive the desired bound for $L(\widehat{U},U)$.

Step 1. By definition of $\widehat{U}^{(1)}$, we have $\mathop{\sf Tr}((\widehat{U}^{(1)})^{\prime}\widehat{\Sigma}_{x}\widehat{U}^{(1)})-2\mathop{\sf Tr}((\widehat{U}^{(1)})^{\prime}\widehat{\Sigma}_{xy}\widehat{V}^{(0)})+\rho_{u}{\sum_{j=1}^{p}||\widehat{U}^{(1)}_{j\cdot}||}\leq\mathop{\sf Tr}((U^{*})^{\prime}\widehat{\Sigma}_{x}U^{*})-2\mathop{\sf Tr}((U^{*})^{\prime}\widehat{\Sigma}_{xy}\widehat{V}^{(0)})+\rho_{u}{\sum_{j=1}^{p}||U^{*}_{j\cdot}||}$. After rearrangement, we have

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

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

where $[\cdot]_{j\cdot}$ means the $j$-th row of the corresponding matrix. When

Since $\sum_{j\in S_{u}}||\Delta_{j\cdot}||\leq\sqrt{s_{u}}\sqrt{\sum_{j\in S_{u}}||\Delta_{j\cdot}||^{2}}$, (64) can be upper bounded by

Step 2. The inequality (64) implies the cone condition

where for a subset $B\subset[p]$, $\Delta_{B*}=(\Delta_{ij}{\mathbf{1}_{\left\{{i\in B,j\in[r]}\right\}}})$, and

In the above derivation, we have used the construction of $J_{k}$ and the cone condition (66). Hence, $\|n^{-1/2}X\Delta\|_{\rm F}\geq\kappa\|\Delta_{\widetilde{S}_{u}*}\|_{\rm F}$ with $\kappa=\sqrt{\phi_{\min}^{\widehat{\Sigma}_{x}}(s_{u}+t)}-3\sqrt{\frac{s_{u}}{t}}\sqrt{\phi_{\max}^{\widehat{\Sigma}_{x}}(t)}$. In view of Lemma 6.5, taking $t=c_{1}s_{u}$ for some sufficiently large constant $c_{1}$, with high probability, $\kappa$ can be lower bounded by a positive constant $\kappa_{0}$ only depending on $M$. Combining with (65), we have

Summing (69) and (70), we have $\|\Delta\|_{\rm F}\leq C\sqrt{s_{u}}\rho$. By Lemma 6.7, we may choose $\rho_{u}\geq\gamma_{u}\sqrt{\frac{r+\log p}{n}}$ for some large $\gamma_{u}$ so that (63) holds with high probability. Hence, $\|\Delta\|_{\rm F}\leq C\sqrt{{s_{u}(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 material ), we obtain the desired bound for $L(\widehat{U},U)$. The proof is complete. ∎

References

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.

There exists an absolute constant $C>0$, such that for all integers $N\geq 12$, $k\leq N/12$ and all $|\mu|\leq 3\sqrt{\eta_{N}\log N}$,

where $h_{\mu,0}=\frac{1}{2}(f_{\mu,0}+f_{\mu,1})$ and $h_{\mu,1}=\delta_{N}f_{\mu,1}+\left(1-\delta_{N}\right)\frac{1}{2}(f_{\mu,0}+f_{\mu,1})$.

Suppose $A\sim\mathcal{G}(N,1/2)$. There exists an absolute constant $C>0$ such that

Recall $\eta_{N}$ defined in (26) and $h_{\mu,0}$ in Lemma 7.1. Let $\nu$ be $N\left(0,\eta_{N}\right)$, and $\bar{\nu}$ be the distribution obtained by restricting $\nu$ on the set $[-3\sqrt{\eta_{N}\log N},3\sqrt{\eta_{N}\log N}]$. Then the $\mu_{i}$’s in (30) are i.i.d. r.v.’s following the distribution $\bar{\nu}$.

Hence, it is sufficient to bound ${\sf TV}(\mathcal{L}(\{W_{i}\}_{i=1}^{n}),\mathcal{L}(\{\overline{W}_{i}\}_{i=1}^{n}))$. Conditioning on $\mu_{i}$, $W_{ij}$ follows $h_{\mu_{i},0}$ when $A\sim\mathcal{G}(N,k)$. Therefore,

Here the last inequality is due to Lemma 7.1. Applying Lemma 7 of , we obtain ${\sf TV}(\mathcal{L}(\{W_{i}\}_{i=1}^{n}),\mathcal{L}(\{\overline{W}_{i}\}_{i=1}^{n}))\leq\sum_{i=1}^{n}\sum_{j=1}^{n}{\sf TV}(W_{ij},\overline{W}_{ij})\leq CN^{-1}$. This completes the proof. ∎