Computational and Statistical Boundaries for Submatrix Localization in a Large Noisy Matrix

Introduction

where $M$ is the signal of interest and $Z$ is noise, is ubiquitous in statistics and is used in a wide range of applications. When $M$ and $Z$ are matrices, many interesting problems arise under a variety of structural assumptions on $M$ and the distribution of $Z$. Examples include sparse principal component analysis (PCA) (Vu and Lei,, 2012; Berthet and Rigollet, 2013b, ; Birnbaum et al.,, 2013; Cai et al.,, 2013, 2015), non-negative matrix factorization (Lee and Seung,, 2001), non-negative PCA (Zass and Shashua,, 2006; Montanari and Richard,, 2014). Under the conventional statistical framework, one is looking for optimal statistical procedures for recovering the signal or detecting its presence.

As the dimensionality of the data becomes large, the computational concerns associated with statistical procedures come to the forefront. In particular, problems with a combinatorial structure or non-convex constraints pose a significant computational challenge because naive methods based on exhaustive search are typically not computationally efficient. Trade-off between computational efficiency and statistical accuracy in high-dimensional inference has drawn increasing attention in the literature. In particular, Chandrasekaran et al., (2012) and Wainwright, (2014) considered a general class of linear inverse problems, with different emphasis on convex geometry and decomposition of statistical and computational errors. Chandrasekaran and Jordan, (2013) studied an approach for trading off computational demands with statistical accuracy via relaxation hierarchies. Berthet and Rigollet, 2013a ; Ma and Wu, (2013); Zhang et al., (2014) focused on computational requirements for various statistical problems, such as detection and regression.

In the present paper, we study the interplay between computational efficiency and statistical accuracy in submatrix localization based on a noisy observation of a large matrix. The problem considered in this paper is formalized as follows.

This model can be further extended to the case of multiple submatrices as

where $|R_{s}|=k_{s}^{(m)}$ and $|C_{s}|=k_{s}^{(n)}$ denote the support set of the $s$-th submatrix. For simplicity, we first focus on the single submatrix and then extend the analysis to the model (3) in Section 2.5.

There are two fundamental questions associated with the submatrix model (2). One is the detection problem: given one observation of the $X$ matrix, decide whether it is generated from a distribution in the submatrix model or from the pure noise model. Precisely, the detection problem considers testing of the hypotheses

The other is the localization problem, where the goal is to exactly recover the signal index sets $R_{m}$ and $C_{n}$ (the support of the mean matrix $M$). It is clear that the localization problem is at least as hard (both computationally and statistically) as the detection problem. As we show in this paper, the localization problem requires larger signal to noise ratio $\lambda/\sigma$, as well as a more detailed exploitation of the submatrix structure.

If the signal to noise ratio is sufficiently large, it is computationally easy to localize the submatrix. On the other hand, if this ratio is small, the localization problem is statistically impossible. To quantify this phenomenon, we identify two distinct thresholds ($\sf SNR_{s}$ and $\sf SNR_{c}$) for $\lambda/\sigma$ in terms of parameters $m,n,k_{m},k_{n}$. The first threshold, $\sf SNR_{s}$, captures the statistical boundary, below which no method (possibly exponential time) can succeed with probability going to one in the minimax sense. The exhaustive search method successfully finds the submatrix above this threshold. The second threshold, $\sf SNR_{c}$, corresponds to the computational boundary, above which an adaptive (with respect to the parameters) linear time spectral algorithm finds the signal. Below this threshold, no polynomial time algorithm can succeed, under the hidden clique hypothesis, described later.

2 Prior Work

There is a growing body of work in statistical literature on submatrix problems. Shabalin et al., (2009) provided a fast iterative maximization algorithm to solve the submatrix localization problem. However, as with many EM type algorithms, the theoretical result is very sensitive to initialization. Arias-Castro et al., (2011) studied the detection problem for a cluster inside a large matrix. Butucea and Ingster, (2013); Butucea et al., (2013) formulated the submatrix detection and localization problems under Gaussian noise and determined sharp statistical transition boundaries. For the detection problem, Ma and Wu, (2013) provided a computational lower bound result under the assumption that hidden clique detection is computationally difficult.

Balakrishnan et al., (2011); Kolar et al., (2011) focused on statistical and computational trade-offs for the submatrix localization problem. They provided a computationally feasible entry-wise thresholding algorithm, a row/column averaging algorithm, and a convex relaxation for sparse SVD to investigate the minimum signal to noise ratio that is required in order to localize the submatrix. Under the sparse regime $k_{m}\precsim m^{1/2}$ and $k_{n}\precsim n^{1/2}$, the entry-wise thresholding turns out to be the “near optimal” polynomial-time algorithm (which we will show a de-noised spectral algorithm that perform slightly better in Section 2.4). However, for the dense regime when $k_{m}\succsim m^{1/2}$ and $k_{n}\succsim n^{1/2}$, the algorithms provided in Kolar et al., (2011) are not optimal in the sense that there are other polynomial-time algorithm that can succeed in finding the submatrix with smaller SNR. Concurrently with our work, Chen and Xu, (2014) provided a convex relaxation algorithm that improves the SNR boundary of Kolar et al., (2011) in the dense regime. On the downside, the implementation of the method requires a full SVD on each iteration, and therefore does not scale well with the dimensionality of the problem. Furthermore, there is no computational lower bound in the literature to guarantee the optimality of the SNR boundary achieved in Chen and Xu, (2014).

A problem similar to submatrix localization is that of clique finding. Deshpande and Montanari, (2013) presented an iterative approximate message passing algorithm to solve the latter problem with sharp boundaries on SNR. However, in contrast to submatrix localization, where the signal submatrix can be located anywhere within the matrix, the clique finding problem requires the signal to be centered on the diagonal.

We would like to emphasize the difference between detection and localization problems. When $M$ is a vector, Donoho and Jin, (2004) proposed the “higher criticism” approach to solve the detection problem under the Gaussian sequence model. Combining the results in (Donoho and Jin,, 2004; Ma and Wu,, 2013), in the computationally efficient region, there is no loss in treating $M$ in model (2) as a vector and applying the higher criticism method to the vectorized matrix for the problem of submatrix detection. In fact, the procedure achieves sharper constants in the Gaussian setting. However, in contrast to the detection problem, we will show that for localization, it is crucial to utilize the matrix structure, even in the computationally efficient region.

3 Notation

Denote the asymptotic notation $a(n)=\Theta(b(n))$ if there exist two universal constants $c_{l},c_{u}$ such that $c_{l}\leq\varliminf\limits_{n\rightarrow\infty}a(n)/b(n)\leq\varlimsup\limits_{n\rightarrow\infty}a(n)/b(n)\leq c_{u}$. $\Theta^{*}$ is asymptotic equivalence hiding logarithmic factors in the following sense: $a(n)=\Theta^{*}(b(n))$ iff there exists $c>0$ such that $a(n)=\Theta(b(n)\log^{c}n)$. Additionally, we use the notation $a(n)\asymp b(n)$ as equivalent to $a(n)=\Theta(b(n))$, $a(n)\succsim b(n)$ iff $\lim_{n\rightarrow\infty}a(n)/b(n)=\infty$ and $a(n)\precsim b(n)$ iff $\lim_{n\rightarrow\infty}a(n)/b(n)=0$.

We define the zero-mean sub-Gaussian random variable $\bf z$ with sub-Gaussian parameter $\sigma$ in terms of its Laplacian. If there exists a universal constant $c>0$,

Clearly, Gaussian and Bernoulli measures, and more general product measures of zero-mean sub-Gaussian random variables satisfy this isotropic definition up to a constant scalar factor.

4 Our Contributions

To state our main results, let us first define a hierarchy of algorithms in terms of their worst-case running time on instances of the submatrix localization problem:

The set ${\sf LinAlg}$ contains algorithms $\mathcal{A}$ that produce an answer (in our case, the localization subset $\hat{R}^{\mathcal{A}}_{m},\hat{C}^{\mathcal{A}}_{n}$) in time linear in $m\times n$ (the minimal computation required to read the matrix). The classes ${\sf PolyAlg}$ and ${\sf ExpoAlg}$ of algorithms, respectively, terminate in polynomial and exponential time, while ${\sf AllAlg}$ has no restriction.

Combining Theorem 3 and 4 in Section 2 and Theorem 5 in Section 3, the statistical and computational boundaries for submatrix localization can be summarized as follows.

Consider the submatrix localization problem under the model (2). The computational boundary ${\sf SNR_{c}}$ for the dense case when $\min\{k_{m},k_{n}\}\succsim\max\{m^{1/2},n^{1/2}\}$ is

where (7) holds under the Hidden Clique hypothesis $\sf HC_{l}$ (see Section 2.1). For the sparse case when $\max\{k_{m},k_{n}\}\precsim\min\{m^{1/2},n^{1/2}\}$, the computational boundary is ${\sf SNR_{c}}=\Theta^{*}(1)$, more precisely

The statistical boundary ${\sf SNR_{s}}$ is

under the minimal assumption $\max\{k_{m},k_{n}\}\precsim\min\{m,n\}$.

If we parametrize the submatrix model as $m=n,k_{m}\asymp k_{n}\asymp k=\Theta^{*}(n^{\alpha}),\lambda/\sigma=\Theta^{*}(n^{-\beta})$, for some $0<\alpha,\beta<1$, we can summarize the results of Theorem 1 in a phase diagram, as illustrated in Figure 1.

To explain the diagram, consider the following cases. First, the statistical boundary is

which gives the line separating the red and the blue regions. For the dense regime $\alpha\geq 1/2$, the computational boundary given by Theorem 1 is

which corresponds to the line separating the blue and the green regions. For the sparse regime $\alpha<1/2$, the computational boundary is $\Theta(1)\precsim{\sf SNR_{c}}\precsim\Theta(\sqrt{\log\frac{m\vee n}{k_{m}k_{n}}})$, which is the horizontal line connecting $(\alpha=0,\beta=0)$ to $(\alpha=1/2,\beta=0)$.

As a key part of Theorem 1, we provide various linear time spectral algorithms that will succeed in localizing the submatrix with high probability in the regime above the computational threshold. Furthermore, the method is adaptive: it does not require the prior knowledge of the size of the submatrix. This should be contrasted with the method of Chen and Xu, (2014) which requires the prior knowledge of $k_{m},k_{n}$; furthermore, the running time of their SDP-based method is superlinear in $nm$. Under the hidden clique hypothesis, we prove that below the computational threshold there is no polynomial time algorithm that can succeed in localizing the submatrix. This is a new result that has not been established in the literature. We remark that the computational lower bound for localization requires a technique different from the lower bound for detection; the latter has been resolved in Ma and Wu, (2013).

Beyond localization of one single submatrix, we generalize both the computational and statistical story to a growing number of submatrices in Section 2.5. As mentioned earlier, the statistical boundary for one single submatrix localization has been investigated by Butucea et al., (2013) in the Gaussian case. Our result focuses on the computational intrinsic difficulty of localization for a growing number of submatrices, at the expense of not providing the exact constants for the thresholds.

The phase transition diagram in Figure 1 for localization should be contrasted with the corresponding result for detection, as shown in (Butucea and Ingster,, 2013; Ma and Wu,, 2013). For a large enough submatrix size (as quantified by $\alpha>2/3$), the computationally-intractable-but-statistically-possible region collapses for the detection problem, but not for localization. In plain words, detecting the presence of a large submatrix becomes both computationally and statistically easy beyond a certain size, while for localization there is always a gap between statistically possible and computationally feasible regions. This phenomenon also appears to be distinct to that of other problems like estimation of sparse principal components (Cai et al.,, 2013), where computational and statistical easiness coincide with each other over a large region of the parameter spaces.

5 Organization of the Paper

The paper is organized as follows. Section 2 establishes the computational boundary, with the computational lower bounds given in Section 2.1 and upper bound results in Sections 2.2-2.4. An extension to the case of multiple submatrices is presented in Section 2.5. The upper and lower bounds for statistical boundary for multiple submatrices are discussed in Section 3. A discussion is given in Section 4. Technical proofs are deferred to Section 5. In addition to the spectral method given in Section 2.2 and 2.4, Appendix A contains a new analysis of a known method that is based on a convex relaxation (Chen and Xu,, 2014). Comparison of computational lower bounds for localization and detection is included in Appendix B.

Computational Boundary

We characterize in this section the computational boundaries for the submatrix localization problem. Sections 2.1 and 2.2 consider respectively the computational lower bound and upper bound. The computational lower bound given in Theorem 2 is based on the hidden clique hypothesis.

Theoretical Computer Science identifies a range problems which are believed to be “hard,” in the sense that in the worst-case the required computation grows exponentially with the size of the problem. Faced with a new computational problem, one might try to reduce any of the “hard” problems to the new problem, and therefore claim that the new problem is as hard as the rest in this family. Since statistical procedures typically deal with a random (rather than worst-case) input, it is natural to seek token problems that are believed to be computationally difficult on average with respect to some distribution on instances. The hidden clique problem is one such example (for recent results on this problem, see Feldman et al., (2013); Deshpande and Montanari, (2013)). While there exists a quasi-polynomial algorithm, no polynomial-time method (for the appropriate regime, described below) is known. Following several other works on reductions for statistical problems, we work under the hypothesis that no polynomial-time method exists.

Let us make the discussion more precise. Consider the hidden clique model $\mathcal{G}(N,\kappa)$ where $N$ is the total number of nodes and $\kappa$ is the number of clique nodes. In the hidden clique model, a random graph instance is generated in the following way. Choose $\kappa$ clique nodes uniformly at random from all the possible choices, and connect all the edges within the clique. For all the other edges, connect with probability $1/2$.

Consider the random instance of hidden clique model $\mathcal{G}(N,\kappa)$. For any sequence $\kappa(N)$ such that $\kappa(N)\leq N^{\beta}$ for some $0<\beta<1/2$, there is no randomized polynomial time algorithm that can find the planted clique with probability tending to $1$ as $N\rightarrow\infty$. Mathematically, define the randomized polynomial time algorithm class $\sf PolyAlg$ as the class of algorithms $\mathcal{A}$ that satisfies

Consider the hidden clique model $\mathcal{G}(N,\kappa)$. For any sequence of $\kappa(N)$ such that $\kappa(N)\leq N^{\beta}$ for some $0<\beta<1/2$, there is no randomized polynomial time algorithm that can distinguish between

with probability going to $1$ as $N\rightarrow\infty$. Here $\mathcal{P}_{\sf ER}$ is the Erdős-Rényi model, while $\mathcal{P}_{\sf HC}$ is the hidden clique model with uniform distribution on all the possible locations of the clique. More precisely,

The hidden clique hypothesis has been used recently by several authors to claim computational intractability of certain statistical problems. In particular, Berthet and Rigollet, 2013a ; Ma and Wu, (2013) assumed the hypothesis ${\sf HC_{d}}$ and Wang et al., (2014) used $\sf HC_{l}$. Localization is harder than detection, in the sense that if an algorithm $\mathcal{A}$ solves the localization problem with high probability, it also correctly solves the detection problem. Assuming that no polynomial time algorithm can solve the detection problem implies impossibility results in localization as well. In plain language, $\sf HC_{l}$ is a milder hypothesis than $\sf HC_{d}$.

We will provide two computational lower bound results, one for localization and the other for detection, in Theorems 2 and 6. The latter one will be deferred to Appendix B to contrast the difference of constructions between localization and detection. The detection computational lower bound was first proved in Ma and Wu, (2013). For the localization computational lower bound, to the best of our knowledge, there is no proof in the literature. Theorem 2 ensures the upper bound in Lemma 1 being sharp.

Consider the submatrix model (2) with parameter tuple $(m=n,k_{m}\asymp k_{n}\asymp n^{\alpha},\lambda/\sigma=n^{-\beta})$, where $\frac{1}{2}<\alpha<1,~{}\beta>0$. Under the computational assumption $\sf HC_{l}$, if

it is not possible to localize the true support of the submatrix with probability going to $1$ within polynomial time.

Our algorithmic reduction for localization relies on a bootstrapping idea based on the matrix structure and a cleaning-up procedure introduced in Lemma 13 given in Section 5. These two key ideas offer new insights in addition to the usual computational lower bound arguments. Bootstrapping introduces an additional randomness on top of the randomness in the hidden clique. Careful examination of these two $\sigma$-fields allows us to write the resulting object into mixture of submatrix models. For submatrix localization we need to transform back the submatrix support to the original hidden clique support exactly, with high probability. In plain language, even though we lose track of the exact location of the support when reducing the hidden clique to submatrix model, we can still recover the exact location of the hidden clique with high probability. For technical details of the proof, please refer to Section 5.

2 Adaptive Spectral Algorithm and Computational Upper Bound

In this section, we introduce linear time algorithm that solves the submatrix localization problem above the computational boundary $\sf SNR_{c}$. Our proposed localization Algorithms 1 and 2 is motivated by the spectral algorithm in random graphs (McSherry,, 2001; Ng et al.,, 2002).

The proposed algorithm has several advantages over the localization algorithms that appeared in literature. First, it is a linear time algorithm (that is, $\Theta(mn)$ time complexity). The top singular vectors can be evaluated using fast iterative power methods, which is efficient both in terms of space and time. Secondly, this algorithm does not require the prior knowledge of $k_{m}$ and $k_{n}$ and automatically adapts to the true submatrix size.

Lemma 1 below justifies the effectiveness of the spectral algorithm.

Consider the submatrix model (2), Algorithm 1 and assume $\min\{k_{m},k_{n}\}\succsim\max\{m^{1/2},n^{1/2}\}$. There exist a universal $C>0$ such that when

the spectral method succeeds in the sense that $\hat{R}_{m}=R_{m},\hat{C}_{n}=C_{n}$ with probability at least $1-m^{-c}-n^{-c}-2\exp\left(-c(m+n)\right)$.

3 Dense Regime

We are now ready to state the SNR boundary for polynomial-time algorithms (under an appropriate computational assumption), thus excluding the exhaustive search procedure. The results hold under the dense regime when $k\succsim n^{1/2}$.

Consider the submatrix model (2) and assume $\min\{k_{m},k_{n}\}\succsim\max\{m^{1/2},n^{1/2}\}$. There exists a critical rate

for the signal to noise ratio $\sf SNR_{c}$ such that for $\lambda/\sigma\succsim{\sf SNR_{c}}$, both the adaptive linear time Algorithm 1 and the robust polynomial time Algorithm 5 will succeed in submatrix localization, i.e., $\hat{R}_{m}=R_{m},\hat{C}_{n}=C_{n}$, with high probability. For $\lambda/\sigma\precsim{\sf SNR_{c}}$, there is no polynomial time algorithm that will work under the hidden clique hypothesis $\sf HC_{l}$.

The proof of the above theorem is based on the theoretical justification of the spectral Algorithm 1 and convex relaxation Algorithm 5, and the new computational lower bound result for localization in Theorem 2. We remark that the analyses can be extended to multiple, even growing number of submatrices case. We postpone a proof of this fact to Section 2.5 for simplicity and focus on the case of a single submatrix.

4 Sparse Regime

Under the sparse regime when $k\precsim n^{1/2}$, a naive plug-in of Lemma 1 requires the ${\sf SNR_{c}}$ to be larger than $\Theta(n^{1/2}/k)\succsim\sqrt{\log n}$, which implies the vanilla spectral Algorithm 1 is outperformed by simple entrywise thresholding. However, a modified version with entrywise soft-thresholding as a preprocessing de-noising step turns out to provide near optimal performance in the sparse regime. Before we introduce the formal algorithm, let us define the soft-thresholding function at level $t$ to be

Soft-thresholding as a de-noising step achieving optimal bias-and-variance trade-off has been widely understood in the wavelet literature, for example, see Donoho and Johnstone, (1998).

Now we are ready to state the following de-noised spectral Algorithm 2 to localize the submatrix under the sparse regime when $k\precsim n^{1/2}$.

Lemma 2 below provides the theoretical guarantee for the above algorithm when $k\precsim n^{1/2}$.

Consider the submatrix model (2), soft-thresholded spectral Algorithm 2 with thresholded level $\sigma t$, and assume $\min\{k_{m},k_{n}\}\precsim\max\{m^{1/2},n^{1/2}\}$. There exist a universal $C>0$ such that when

the spectral method succeeds in the sense that $\hat{R}_{m}=R_{m},\hat{C}_{n}=C_{n}$ with probability at least $1-m^{-c}-n^{-c}-2\exp\left(-c(m+n)\right)$. Further if we choose $t=\Theta(\sigma\sqrt{\log\frac{m\vee n}{k_{m}k_{n}}})$ as the optimal thresholding level, we have de-noised spectral algorithm works when

Combining the hidden clique hypothesis $\sf HC_{l}$ together with Lemma 2, we have the following theorem holds under the sparse regime when $k\precsim n^{1/2}$.

Consider the submatrix model (2) and assume $\max\{k_{m},k_{n}\}\precsim\min\{m^{1/2},n^{1/2}\}$. There exists a critical rate for the signal to noise ratio $\sf SNR_{c}$ between

such that for $\lambda/\sigma\succsim\sqrt{\log\frac{m\vee n}{k_{m}k_{n}}}$, the linear time Algorithm 2 will succeed in submatrix localization, i.e., $\hat{R}_{m}=R_{m},\hat{C}_{n}=C_{n}$, with high probability. For $\lambda/\sigma\precsim 1$, there is no polynomial time algorithm that will work under the hidden clique hypothesis $\sf HC_{l}$.

The upper bound achieved by the de-noised spectral Algorithm 2 is optimal in the two boundary cases: $k=1$ and $k\asymp n^{1/2}$. When $k=1$, both the information theoretic and computational boundary meet at $\sqrt{\log n}$. When $k\asymp n^{1/2}$, the computational lower bound and upper bound match in Theorem 4, thus suggesting the near optimality of Algorithm 2 within the polynomial time algorithm class. The potential logarithmic gap is due to the crudeness of the hidden clique hypothesis. Precisely, for $k=2$, hidden clique is not only hard for $G(n,p)$ with $p=1/2$, but also hard for $G(n,p)$ with $p=1/\log n$. Similarly for $k=n^{\alpha},\alpha<1/2$, hidden clique is not only hard for $G(n,p)$ with $p=1/2$, but also for some $0<p<1/2$.

5 Extension to Growing Number of Submatrices

The computational boundaries established in the previous sections for a single submatrix can be extended to non-overlapping multiple submatrices model (3). The non-overlapping assumption corresponds to that for any $1\leq s\neq t\leq r$, $R_{s}\cap R_{t}=\emptyset$ and $C_{s}\cap C_{t}=\emptyset$. The Algorithm 3 below is an extension of the spectral projection Algorithm 1 to address the multiple submatrices localization problem.

We emphasize that the following Proposition 3 holds even when the number of submatrices $r$ grows with $m,n$.

Consider the non-overlapping multiple submatrices model (3) and Algorithm 3. Assume

for all $1\leq s\leq r$ and $\min\{k_{m},k_{n}\}\succsim\max\{m^{1/2},n^{1/2}\}$. There exist a universal $C>0$ such that when

the spectral method succeeds in the sense that $\hat{R}_{m}^{(s)}=R_{m}^{(s)},\hat{C}_{n}^{(s)}=C_{n}^{(s)},1\leq s\leq r$ with probability at least $1-m^{-c}-n^{-c}-2\exp\left(-c(m+n)\right)$.

Under the non-overlapping assumption, $rk_{m}\precsim m,~{}rk_{n}\precsim n$ hold in most cases. Thus the first term in Equation (12) is dominated by the latter two terms. Thus a growing number $r$ does not affect the bound in Equation (12) as long as the non-overlapping assumption holds.

Statistical Boundary

In this section we study the statistical boundary. As mentioned in the introduction, in the Gaussian noise setting, the statistical boundary for a single submatrix localization has been established in Butucea et al., (2013). In this section, we generalize to localization of a growing number of submatrices, as well as sub-Gaussian noise, at the expense of having non-exact constants for the threshold.

We begin with the information theoretic lower bound for the localization accuracy.

Consider the submatrix model (2) with Gaussian noise $Z_{ij}\sim\mathcal{N}(0,\sigma^{2})$. For any fixed $0<\alpha<1$, there exist a universal constant $C_{\alpha}$ such that if

any algorithm $\mathcal{A}$ will fail to localize the submatrix with probability at least $1-\alpha-\frac{\log 2}{k_{m}\log(m/k_{m})+k_{n}\log(n/k_{n})}$ in the following minimax sense:

2 Combinatorial Search for Growing Number of Submatrices

Combinatorial search over all submatrices of size $k_{m}\times k_{n}$ finds the location with the strongest aggregate signal and is statistically optimal (Butucea et al.,, 2013; Butucea and Ingster,, 2013). Unfortunately, it requires computational complexity $\Theta\left(\binom{m}{k_{m}}+\binom{n}{k_{n}}\right)$, which is exponential in $k_{m},k_{n}$. The search Algorithm 4 was introduced and analyzed under the Gaussian setting for a single submatrix in Butucea and Ingster, (2013), which can be used iteratively to solve multiple submatrices localization.

For the case of multiple submatrices, the submatrices can be extracted with the largest sum in a greedy fashion.

Lemma 5 below provides a theoretical guarantee for Algorithm 4 to achieve the information theoretic lower bound.

Consider the non-overlapping multiple submatrices model (3) and iterative application of Algorithm 4 in a greedy fashion for $r$ times. Assume

for all $1\leq s\leq r$ and $\max\{k_{m},k_{n}\}\precsim\min\{m,n\}$. There exists a universal constant $C>0$ such that if

then Algorithm 4 will succeed in returning the correct location of the submatrix with probability at least $1-\frac{2k_{m}k_{n}}{mn}$.

To complete Theorem 1, we include the following Theorem 5 capturing the statistical boundary. It is proved by exhibiting the information-theoretic lower bound Lemma 4 and analyzing Algorithm 4.

Consider the submatrix model (2). There exists a critical rate

for the signal to noise ratio, such that for any problem with $\lambda/\sigma\succsim{\sf SNR_{s}}$, the statistical search Algorithm 4 will succeed in submatrix localization, i.e., $\hat{R}_{m}=R_{m},\hat{C}_{n}=C_{n}$, with high probability. On the other hand, if $\lambda/\sigma\precsim{\sf SNR_{s}}$, no algorithm will work (in the minimax sense) with probability tending to $1$.

Discussion

In this paper we established the computational and statistical boundaries for submatrix localization in the setting of a growing number of submatrices with subgaussian noise. The primary goals are to demonstrate the intrinsic gap between what is statistical possible and what is computationally feasible and to contrast the interplay between computational efficiency and statistical accuracy for localization with that for detection.

As pointed out in Section 1.4, for any $k=n^{\alpha},0<\alpha<1$, there is an intrinsic SNR gap between computational and statistical boundaries for submatrix localization. Unlike the submatrix detection problem where for the regime $2/3<\alpha<1$, there is no gap between what is computationally possible and what is statistical possible. The inevitable gap in submatrix localization is due to the combinatorial structure of the problem. This phenomenon is also seen in some network related problems, for instance, stochastic block models with a growing number of communities. Compared to the submatrix detection problem, the algorithm to solve the localization problem is more complicated and the techniques required for the analysis are much more involved.

Detection for Growing Number of Submatrices

The current paper solves localization of a growing number of submatrices. In comparison, for detection, the only known results are for the case of a single submatrix as considered in Butucea and Ingster, (2013) for the statistical boundary and in Ma and Wu, (2013) for the computational boundary. The detection problem in the setting of a growing number of submatrices is of significant interest. In particular, it is interesting to understand the computational and statistical trade-offs in such a setting. This will need further investigation.

Estimation of the Noise Level σ𝜎\sigma

Although Algorithms 1 and 3 do not require the noise level $\sigma$ as an input, Algorithm 2 does require the knowledge of $\sigma$. The noise level $\sigma$ can be estimated robustly. In the Gaussian case, a simple robust estimator of $\sigma$ is the following median absolute deviation (MAD) estimator due to the fact that $M$ is sparse:

Proofs

We prove in this section the main results given in the paper. We first collect and prove a few important technical lemmas that will be used in the proofs of the main results.

We start with two Lemmas 6 and 7 that due to perturbation theory.

Then suppose there is a number $\delta>0$ such that

Further, suppose there are numbers $\alpha,\delta$ such that

then for $2$-norm, or any unitarily invariant norm, we have

Let us use the above version of the perturbation bound to derive a lemma that is particularly useful in our case. Simple algebra tells us that

and similarly we have (since the operator norm of a whole matrix is larger than that of the submatrix)

Thus the following version of the Wedin’s Theorem holds.

and also the following holds for $2$-norm (or any unitary invariant norm)

We will then introduce some concentration inequalities. Lemmas 8 and 9 are concentration of measure results from random matrix theory.

where $C,c>0$ are some universal constants.

The proof of this lemma is a simple application of Theorem 2.1 inHsu et al., (2012) for the case that $\mathcal{P}$ is a rank $r$ positive semidefinite projection matrix.

The following two are standard Chernoff-type bounds for bounded random variables.

Let $X_{i},1\leq i\leq n$ be independent random variables. Assume $a_{i}\leq X_{i}\leq b_{i},1\leq i\leq n$. Then for $S_{n}=\sum_{i=1}^{n}X_{i}$

Let $X_{i},1\leq i\leq n$ be independent zero-mean random variables. Suppose $|X_{i}|\leq M,1\leq i\leq n$. Then

We will end this section stating the Fano’s information inequality, which plays a key role in many information theoretic lower bounds.

Let $\mathcal{P}_{0},\mathcal{P}_{1},\ldots,\mathcal{P}_{M}$ be probability measures on the same probability space $(\Theta,\mathcal{F})$, $M\geq 2$. If for some $0<\alpha<1$

where $p_{e,M}$ is the minimax error for the multiple testing problem.

2 Main Proofs

Recall the matrix form of the submatrix model, with the SVD decomposition of the mean signal matrix $M$

Using Weyl’s interlacing inequality, we have

And according to the definition of the canonical angles, we have

Define the projection operator to be $\mathcal{P}$, we start the analysis by decomposing

We invoke the union bound for all $1\leq j\leq n$ to obtain

In the second approach, the second term of (31) can be handled through perturbation Sin Theta Theorem 7:

This second approach will be used in the multiple submatrices analysis.

Combining all the above, we have with probability at least $1-n^{-c}-m^{-c}$, for all $1\leq j\leq n$

Similarly we have for all $1\leq i\leq m$,

Clearly we know that for $i\in R_{m}$ and $i^{\prime}\in[m]\backslash R_{m}$

and for $j\in C_{n}$ and $j^{\prime}\in[n]\backslash C_{n}$

hold, then we have learned a metric $d$ (a one dimensional line) such that on this line, data forms clusters in the sense that

In this case, a simple cut-off clustering recovers the nodes exactly.

the spectral algorithm succeeds with probability at least

The proof of the validity of thresholded spectral algorithm at level $\sigma t$ is easy based on the proof of Lemma 1. Firstly we have the following decomposition

Let us prove this fact. Clearly if $(i,j)\notin R_{m}\times C_{n}$, $B_{ij}=0$. If $(i,j)\in R_{m}\times C_{n}$, we have

where the last step uses $|\eta_{t}(y)-y|\leq t$, for any $y$. Let us bound the variance of each thresholded entry $\eta_{\sigma t}(Z_{ij})$,

for some universal constant $C$. Clearly after thresholding, $\eta_{\sigma t}(Z)$ still have i.i.d entries, but the variance has been significantly reduced as $t\rightarrow\infty$.

Via the perturbation analysis established in Proof of Lemma 1

as $B$ only have $k_{m}k_{n}$ non zero entries. Thus applying Lemma 7, we have

As usual, we continue the analysis by decomposing (following the steps as in Lemma 7, but with an additional bias term $B$)

for $1\leq j\leq n$. We know for $j\in C_{n}$ and $j^{\prime}\in[n]\backslash C_{n}$

hold, then we have learned a metric $d$ (a one dimensional line) such that on this line, data forms clusters. In this case, a simple cut-off clustering recovers the nodes exactly.

the thresholded spectral algorithm succeeds with probability at least

Computational lower bound for localization (support recovery) is of different nature than the computational lower bound for detection (two point testing). The idea is to design a randomized polynomial time algorithmic reduction to relate a an instance of hidden clique problem to our submatrix localization problem. The proof proceeds in the following way: we will construct a randomized polynomial time transformation $\mathcal{T}$ to map a random instance of $\mathcal{G}(N,\kappa)$ to a random instance of our submatrix $\mathcal{M}(m=n,k_{m}\asymp k_{n}\asymp k,\lambda/\sigma)$ (abbreviated as $\mathcal{M}(n,k,\lambda/\sigma)$). Then we will provide a quantitative computational lower bound by showing that if there is a polynomial time algorithm that pushes below the hypothesized computational boundary for localization in the submatrix model, there will be a polynomial time algorithm that solves hidden clique localization with high probability (a contradiction to $\sf HC_{l}$).

Denote the randomized polynomial time transformation as

There are several stages for the construction of the algorithmic reduction. First we define a graph $\mathcal{G}^{e}(N,\kappa(N))$ that is stochastically equivalent to the hidden clique graph $\mathcal{G}(N,\kappa(N))$, but is easier for theoretical analysis. $\mathcal{G}^{e}$ has the property: each node independently has the probability $\kappa(N)/N$ to be a clique node, and with the remaining probability a non-clique node. Using Bernstein’s inequality and the inequality (46) proved below. with probability at least $1-2N^{-1}$ the number of clique nodes $\kappa^{e}$ in $\mathcal{G}^{e}$

Consider a hidden clique graph $\mathcal{G}^{e}(2N,2\kappa(N))$ with $N=n$ and $\kappa(N)=\kappa$. Denote the set of clique nodes for $\mathcal{G}^{e}(2N,2\kappa(N))$ to be $C_{N,\kappa}$. Represent the hidden clique graph using the symmetric adjacency matrix $G\in\{-1,1\}^{2N\times 2N}$, where $G_{ij}=1$ if $i,j\in C_{N,\kappa}$, otherwise with equal probability to be either $-1$ or $1$. As remarked before, with probability at least $1-2N^{-1}$, we have planted $2\kappa(1\pm o(1))$ clique nodes in graph $\mathcal{G}^{e}$ with $2N$ nodes. Take out the upper-right submatrix of $G$, denote as $G_{UR}$ where $U$ is the index set $1\leq i\leq N$ and $R$ is the index set $N+1\leq j\leq 2N$. Now $G_{UR}$ has independent entries.

Due to the bootstrapping property, the matrices $\left[(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j)}\right]_{1\leq i,j\leq n}$, indexed by $0\leq s,t<l$ are independent of each other. Recall that $C_{N,\kappa}$ stands for the clique set of the hidden clique graph. We define the row candidate set $R_{l}:=\{i\in[n]:\exists~{}0\leq s<l,\psi^{(s)}(i)\in C_{N,\kappa}\}$ and column candidate set $C_{l}:=\{j\in[n]:\exists~{}0\leq t<l,\phi^{(t)}(j)\in C_{N,\kappa}\}$. Observe that $R_{l}\times C_{l}$ are the indices where the matrix $M$ contains signal.

Now let us discuss the independence issue in $M$ through our Bootstrapping construction. Clearly due to sampling with replacement and bootstrapping, condition on $\mathcal{G}^{e}$, we have independence among samples for the same location $(i,j)$

For the independence among entries in one Bootstrapped matrix, clearly

The only case where there might be a slight dependence is between $(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j)}$ and $(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j^{\prime})}$. The way to eliminate the slight dependence is through Vu, (2008)’s result on universality of random discrete graphs. Vu, (2008) showed random regular graph $\mathcal{G}(n,n/2)$ shares many similarities as Erdős-Rényi random graph $\mathcal{G}(n,1/2)$, for instance, top and second eigenvalues ($n/2$ and $\sqrt{n}$ respectively), limiting spectral distribution, sandwich conjecture, determinant, etc. Let us consider the case where the upper-right of the adjacency matrix $G$ consists of random bi-regular graph (see Deshpande and Montanari, (2013) for difficulty of clique problem under random regular graph) with degree $n/2$ instead of the Erdős-Rényi graph. The only thing we need to change is assuming hidden clique hypothesis is still valid for the following random graph: for a $n\times n$ adjacency matrix $G$, first find a clique/principal submatrix of size $k$ uniformly randomly and connect density, for the remaining part of the matrix, sample a random regular graph of $G(n-k,\frac{n-k}{2})$ and a random bi-regular graph of size $k\times(n-k)$ with left regular degree $n/2-k$ and right regular degree $k/2$ (here degree test will not work in this graph and spectral barrier still suggests $k\precsim\sqrt{n}$ is hard due to universality result of random discrete graphs). In the bootstrapping step, condition on the same row $\psi^{(s)}(i)$ being not a clique, $(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j)}\perp(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j^{\prime})}|\psi^{(s)}(i)$, and each one is a Rademacher random variable (regardless of the choice of $\psi^{(s)}(i)$), which implies $(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j)}\perp(G_{UR})_{\psi^{(s)}(i)\phi^{(t)}(j^{\prime})}$ holds unconditionally. Thus in the bootstrapping procedure, we have independence among entries within the matrix.

Let us move to verify the sub-Gaussianity of $M$ matrix. Note that for the index $i,j$ that is not a clique for any of the matrices, $M_{ij}$ is sub-Gaussian, due to Hoeffding’s inequality

For the index $i,j$ being a clique in at least one of the matrices, we claim the number of matrices has $(i,j)$ being clique is $O^{*}(1)$. Due to Bernstein’s inequality, we have $\max_{i}|\{0\leq s<l:\psi^{(s)}(i)\in C_{N,\kappa}\}|\leq\frac{\kappa l}{n}+\frac{8}{3}\log n$ with probability at least $1-n^{-1}$. This further implies there are at least $l^{2}-(\frac{\kappa l}{n}+\frac{8}{3}\log n)^{2}$ many independent Rademacher random variables in each $i,j$ position, thus

Let us estimate the corresponding $k$ in the submatrix model. We need to bound the order of the cardinality of $R_{l}$, denoted as $|R_{l}|$. The total number of positions with signal (at least one clique node inside) is

which is of the order $k:=\kappa l$. Let us provide a high probability bound on $|R_{l}|$. By Bernstein’s inequality

Thus if we take $u=\sqrt{4\kappa l\log n}$, as long as $\log n=o(\kappa l)$,

So with probability at least $1-2n^{-1}$, the number of positions that contain signal nodes is bounded as

Equation (47) implies that with high probability

The above means, in the submatrix parametrization, $k_{m}\asymp k_{n}\asymp\kappa l\asymp n^{\alpha}$, $\lambda/\sigma\asymp l^{-1}\asymp n^{-\beta}$, which implies $\kappa\asymp n^{\alpha-\beta}$.

Suppose there exists a polynomial time algorithm $\mathcal{A}_{M}$ that pushes below the computational boundary. In other words,

with the last inequality having a slack $\epsilon>0$. More precisely, $\mathcal{A}_{M}$ returns two estimated index sets $\hat{R}_{n}$ and $\hat{C}_{n}$ corresponding to the location of the submatrix (and correct with probability going to $1$) under the regime $\beta=\alpha-1/2+\epsilon$. Suppose under some conditions, this algorithm $\mathcal{A}_{M}$ can be modified to a randomized polynomial time algorithm $\mathcal{A}_{\mathcal{G}}$ that correctly identifies the hidden clique nodes with high probability. It means in the corresponding hidden clique graph $\mathcal{G}(2N,2\kappa)$, $\mathcal{A}_{\mathcal{G}}$ also pushes below the computational boundary of hidden clique by the amount $\epsilon$:

In summary, the quantitative computational lower bound implies that if the computational boundary for submatrix localization is pushed below by an amount $\epsilon$ in the power, the hidden clique boundary is correspondingly improved by $\epsilon$.

Now let us show that any algorithm $\mathcal{A}_{M}$ that localizes the submatrix introduces a randomized algorithm that finds the hidden clique nodes with probability tending to 1. The algorithm relies on the following simple lemma.

For the hidden clique model $\mathcal{G}(N,\kappa)$, suppose an algorithm provides a candidate set $S$ of size $k$ that contains the true clique subset exactly. If

then by looking at the adjacency matrix restricted to $S$ we can recover the clique subset exactly with high probability.

The proof of Lemma 13 is immediate. If $i$ is a clique node, then $\min_{i}\sum_{j\in C}G_{ij}\geq\kappa-C/2\cdot\sqrt{k\log N}$. If $i$ is not a clique node, then $\max_{i}\sum_{j\in C}G_{ij}\leq C/2\cdot\sqrt{k\log N}$. The proof is completed.

Algorithm $\mathcal{A}_{M}$ provides candidate sets $R_{l},C_{l}$ of size $k$, inside which $\kappa$ are correct clique nodes, and thus exact recovery can be completed through Lemma 13 since $\kappa\succsim(k\log N)^{1/2}$ (since $\kappa\asymp n^{1/2-\epsilon}\succsim k^{1/2}\asymp n^{\alpha/2}$ when $\epsilon$ is small). The algorithm $\mathcal{A}_{M}$ induces another randomized polynomial time algorithm $\mathcal{A}_{\mathcal{G}}$ that solves the hidden clique problem $\mathcal{G}(2N,2\kappa)$ with $\kappa\precsim N^{1/2}$. The algorithm $\mathcal{A}_{\mathcal{G}}$ returns the support $\hat{C}_{N,\kappa}$ that coincides with the true support $C_{N,\kappa}$ with probability going to $1$ (a contradiction to the hidden clique hypothesis $\sf HC_{l}$). We conclude that, under the hypothesis, there is no polynomial time algorithm $\mathcal{A}_{M}$ that can push below the computational boundary $\lambda\precsim\sqrt{\frac{m+n}{k_{m}k_{n}}}$.

where $M=\lambda\cdot UV^{T}\in\Theta$. The parameter space $\Theta$ is composed of all $M=\lambda\cdot UV^{T}$ where $U$ are sampled uniformly on the collection of vectors with $k_{m}$ ones and other coordinates being zero, and similarly $V$ are sampled uniformly with $k_{n}$ ones and the rest zero. The cardinality of the parameter space is

corresponding to that many probability measures on the same probability space. Put a uniform prior on this parameter space and invoke Fano’s lemma 12. To obtain the lower bound, we need to upper bound the Kullback-Leibler divergence $d_{\sf KL}(\mathcal{P}_{M}||\bar{\mathcal{P}})$ for any $M\in\Theta$, where

Invoke the simple bound on binomial coefficients $\left(\frac{n}{k}\right)^{k}\leq\binom{n}{k}\leq\left(\frac{ne}{k}\right)^{k}$. If we choose

then the condition (28) holds. Any submatrix localization algorithm translates into a multiple testing procedure that picks a parameter $M^{\prime}\in\Theta$. By Fano’s information inequality 12, the minimax error, which is also the localization error, is at least $1-\alpha-\frac{\log 2}{k_{m}\log(m/k_{m})+k_{n}\log(n/k_{n})}$. ∎

Recall the definition 4 of a sub-Gaussian random variable. Taking $Z_{ij},i\in I,j\in J$, we have the following concentration from the Chernoff’s bound for $\sum_{i\in I,j\in J}Z_{ij}$

such submatrices, so by a union bound, we have

If we take $t=2\sqrt{\frac{k_{m}\log(em/k_{m})+k_{n}\log(en/k_{n})}{c}}$, then with probability at least

then the maximum submatrix is unique and is the true one. Recollecting terms, we reach

To make the proof fully rigorous, we need the following monotonicity trick. Consider the submatrix of size $k_{m}k_{n}$ with $a$ rows to be in the correct set $R_{m}$ and $b$ columns to be in the correct set $C_{n}$, where $a<k_{m}$ and $b<k_{n}$. The cardinality of the set of such matrices is

Using the same calculation as before we want

Hence, if equation (58) is satisfied, (59) is satisfied up to a universal constant for all $a<k_{m}$ and $b<k_{n}$. Thus we have proved that if

with a suitable constant $C$, the statistical search algorithm picks out the correct submatrix. The sum of the probabilities of the bad events is bounded by

For the multiple non-overlapping submatrices case, as long as

then sequential application of Algorithm 4 will find the $r$-submatrices.

We are going to provide theoretical justification to the extension of the submatrix localization algorithm to multiple non-overlapping submatrices case as in Algorithm 3. Write out the matrix form of the submatrix model, with the SVD version of the signal matrix $M$

Due to the non-overlapping property, we have $1\leq s\neq t\leq r$, $1_{R_{s}}^{T}1_{R_{t}}=0$, so as to $C_{n}$. The singular values of $U\Lambda V^{T}$ are $\lambda_{s}\sqrt{k_{s}^{(m)}k_{s}^{(n)}},1\leq s\leq r$, and all the other singular values are .

According to the definition of the canonical angles, we have

Thus invoke the union bound for all $1\leq j\leq n$

Combining all the above, we have with probability at least $1-n^{-c}-m^{-c}$, for all $1\leq j\leq n$

Similarly we have for all $1\leq i\leq m$,

Clearly we know for any $1\leq s\leq r$ and $i\in R_{s}$ and $i^{\prime}\in[m]\backslash R_{s}$

and for any $1\leq s\leq r$ and $j\in C_{s}$ and $j^{\prime}\in[n]\backslash C_{s}$

Thus if $k_{s}^{(m)}\asymp k_{m},k_{s}^{(n)}\asymp k_{n}$ and $\lambda_{s}\asymp\lambda$ for all $1\leq s\leq r$

We have learned a metric $d$ (of intrinsic dimension $r$) such that under this metric, data forms into clusters in the sense that

Thus it satisfies the geometric separation property.

the spectral algorithm succeeds with probability at least

because $rk_{m}\precsim m,~{}rk_{n}\precsim n$ in most cases, the first term does not have an effect in. most cases. ∎

Proof of Theorem 3 is a direct result of Lemma 1 and Theorem 2. Proof of Theorem 4 is obvious based on Lemma 2 and the hidden clique hypothesis $\sf HC_{l}$. Proof of Theorem 5 combines the result of Lemma 5 and Lemma 4.

References

Appendix A Convex Relaxation Algorithm

In this section we will investigate a convex relaxation approach to the problem. The same algorithm has also been investigated in a parallel work of Chen and Xu, (2014). Our analysis is slightly different, with the explicit construction of the dual certificate using the idea in Gross, (2011). For the purposes of comparing to the spectral approach, we include the convex relaxation analysis in this section. Let us write the optimization problem

This problem is non-convex: the feasibility set is non-convex, and so is the optimization function (although it is bi-convex). However, we can relax the problem and transform it into a convex optimization problem. Of course, we need to ensure that the solution to the relaxed problem is the exact solution (with high probability) under appropriate conditions.

The matrix version of the submatrix problem suggests that the signal matrix is of the structure “low rank and sparsity on the singular vectors.” We recall from the low rank matrix recovery literature, see e.g. Candes and Plan, (2011) and Cai et al., (2014), that we can utilize the low rank structure and solve relaxed versions as follows.

Consider the constraint minimization relaxation,

Unfortunately, Relaxation 1 is only good in terms of estimation of the whole matrix. The stronger objective of localization requires simultaneous exploitation of sparsity and low rank-ness, as in Relaxation 2.

Relaxation 2

Let us expand the objective of the original non-convex optimization problem, drop the quadratic term to make the procedure adaptive in terms of $\lambda$, and convexify the feasibility set at the same time.

The time complexity to solve this convex optimization problem is at least $\Theta((m+n)^{3})$ implemented with alternating direction methods of multipliers (ADMM). The disadvantage is that the theoretical guarantee only holds for the exact solution $\hat{X}_{0}$; however, in reality we can only approximately find $\hat{X}_{0}$ through ADMM or some other optimization methods. We also remark that this algorithm requires the prior knowledge of the submatrix size $k_{m},k_{n}$, which means it is not fully adaptive.

Consider the submatrix model (2) and the Algorithm 5. There exists a universal $C>0$ such that when

the convex relaxation succeeds (in the sense that $\hat{R}_{m}=R_{m},\hat{C}_{n}=C_{n}$) with probability at least $1-2m^{-c}-2n^{-c}-2(mn)^{-c}-2\exp\left(-c(m+n)\right)$.

Let us construct the dual certificate to secure that the true solution is the unique solution. If we can construct a pair of a primal certificate $M^{*}=\sqrt{k_{m}k_{n}}UV^{T}$ and dual certificate $\Delta^{*},\Theta^{*},\mu^{*},\nu^{*}$ that satisfy

Here $W\in\mathcal{P}_{U^{\perp}\cap V^{\perp}},\|W\|_{2}\leq 1$, and $UV^{T}+W$ denotes the sub differential of $\|\cdot\|_{*}$ evaluated at $M^{*}$

We claim that any solution $\hat{M}$ to Relaxation 2 must satisfy $\hat{M}=M^{*}$. If not, we write $\hat{M}=M^{*}+H$ and find that

All of the above equations are due to the primal feasibility, and the second inequality also uses the convexity of $\|\cdot\|_{*}$. Note that (79) can be written in a more explicit form

Due to the optimality of $\hat{M}$ in terms of the objective function

We can see that if $H_{ij}>0$, then we must have $M^{*}_{ij}=0$, which through complimentary slackness implies $\Theta^{*}_{ij}=0$. This in turn means $-\Delta_{ij}<0$. When $H_{ij}<0$, we must have $M^{*}_{ij}=1$, which again means $\Delta_{ij}=0$, $\Theta^{*}_{ij}>0$, a contradiction. Thus $H_{ij}=0$ for all $i,j$.

Now let us see how to construct the dual certificate $-\Delta^{*}+\Theta^{*}$, $\mu^{*},\nu^{*}$ to satisfy the conditions (74) - (76). Expand (87) as

Thus if we choose $\mu^{*}\geq C\cdot\sigma(\sqrt{m}+\sqrt{n})$, it holds that

with probability at least $1-2\exp\left(-c(m+n)\right)$. Thus with the choice of $W$, Equation (93) becomes

Now let us write out the explicit form of the projection $\mathcal{P}_{U\cup V}(Z)$

Let us see the concentration property of $[\mathcal{P}_{U\cup V}(Z)]_{ij}$:

with probability at least $1-2n^{-c}$. For all the $1\leq i\leq m$

with probability at least $1-2m^{-c}$. For all $i,j$,

with probability at least $1-2(mn)^{-c}$. Now, pick the dual certificate variables in the following way

where the last two equation follows from (97) and (98). We conclude that the relaxation algorithm succeeds with probability at least

We have achieved the the same boundary as the spectral method upper bound. ∎

Appendix B Algorithmic Reduction for Detection

Consider the submatrix model (2) with parameter tuple $(m=n,k_{m}\asymp k_{n}\asymp n^{\alpha},\lambda/\sigma=n^{-\beta})$, where $\frac{1}{2}<\alpha<1,~{}\beta>0$. Under the hardness assumption $\sf HC_{d}$, if

it is not possible to detect the true support of the submatrix with probability going to $1$ for any polynomial algorithm.

We would like to build a randomized polynomial mapping from the hidden clique graph $\mathcal{G}(N,\kappa(N))$ to a matrix $M(m=n,k_{m}\asymp k_{n}\asymp k,\lambda/\sigma)$ for the submatrix model. Denote this transformation as

There are several stages of the construction. First, we define a graph that is stochastically equivalent to the hidden clique graph $\mathcal{G}$, but is easier for the analysis. Let us call it $\mathcal{G}^{e}$. $\mathcal{G}^{e}$ has the property: each node independently has the probability $\kappa(N)/N$ to be a clique node. By Bernstein’s inequality, with probability at least $1-2N^{-1}$, the number of cliques $\kappa^{e}$ in $\mathcal{G}^{e}$

Consider a double sized hidden clique graph $\mathcal{G}^{e}(2N,2\kappa(N))$ with $N=n^{1+\beta}$, and $\kappa(N)=k=n^{\alpha}$, $\frac{1}{2}<\alpha<1$. Denote the clique nodes set as $C_{N,\kappa}$. Connect the hidden clique graph to form a symmetric matrix $G\in\{-1,1\}^{2N\times 2N}$, where $G_{ij}=1$ if $i,j\in C_{N,\kappa}$, otherwise with equal probability to be either $-1$ or $1$. Take out the upper-right submatrix of $G$, $G_{UR}$ where $U$ is the index set $1\leq i\leq N$ and $R$ is the index set $N+1\leq j\leq 2N$.

For the $s,t$ with clique nodes inside, we know the maximum number of clique nodes is $\log n$ (due to Bernstein’s inequality that $\max_{1\leq s\leq n}|I_{s}\cap C_{N,\kappa}|\leq\frac{k}{n}+\frac{8}{3}\log n$ with probability at least $1-n^{-1}$), which means there are at least $n^{\beta}-\frac{k}{n}-\frac{8}{3}\log n$ many independent Rademacher random variables in each $s,t$ block, thus

Thus we can take sub-Gaussian parameter to be any $\sigma<1$ because $kn^{-1-\beta},n^{-\beta}\log n$ are both $o(1)$. Now this constructed $M(n,k)$ matrix satisfies the submatrix model with $\lambda=n^{-\beta}$ and sub-Gaussian parameter $\sigma=1-o(1)$.

Let us see how many elements in $1\leq s\leq n$ are such that $I_{s}\cap C_{N,\kappa}\neq\emptyset$. Namely, we want to estimate how many clique nodes there exist in the transformed submatrix model. We have the two sided bound

which is of the order $k$. Using Bernstein’s bound, we have with high probability

Thus the submatrix model $M(m=n,k_{m}\asymp k_{n},\lambda/\sigma)$ satisfies $k_{m}\asymp k_{n}\asymp k$.

Suppose there exists an polynomial time algorithm $\mathcal{A}_{M}$ that pushes below the computational boundary quantitatively by a small $\epsilon>0$ amount

where $\beta=2\alpha-1+\epsilon$. Namely, $\mathcal{A}_{M}$ detects below the boundary $\sigma\frac{m+n}{k_{m}k_{n}}$, then it naturally introduced a polynomial time detection algorithm for hidden clique problem $\mathcal{G}(N=n^{1+\beta},\kappa=n^{\alpha})$, which violates the $\sf HC_{d}$ because