Detection and Recovery of Hidden Submatrices

We start by presenting our upper bounds. To that end, we propose three algorithms and analyze their performance. Define the statistics,

The statistics in (7) amounts to adding up all the elements of $\mathsf{X}$, while (8) and (9) enumerate all $k\times k$ submatrices of $\mathsf{X}$ in ${\cal K}_{k,1,n}$ and ${\cal K}_{k,1,n}^{\mathsf{con}}$, and take the submatrix with the maximal sum of entries, respectively. Fix $\delta>0$. Then, our tests are defined as,

where the thresholds are given by $\tau_{\mathsf{sum}}\triangleq\frac{mk^{2}\lambda}{2}$, $\tau^{\mathsf{SD}}_{\mathsf{scan}}\triangleq\sqrt{(4+\delta)k^{2}\log\binom{n}{k}}$, and $\tau^{\mathsf{CSD}}_{\mathsf{scan}}\triangleq\sqrt{(4+\delta)k^{2}\log{n}}$, and correspond roughly to the average between the expected values of each of the statistics in (7)–(9) under the null and alternative hypotheses. It should be emphasized that the tests in (10)–(11) were proposed in, e.g., [37, 11, 42], for the single planted submatrix detection problem.

A few important remarks are in order. First, note that in the scan test, we search for a single planted matrix rather than $m$ such matrices. Second, the sum test exhibits polynomial computational complexity, of $O(n^{2})$ operations, and hence efficient. The scan test in (11), however, exhibits an exponential computational complexity, and thus is inefficient. Indeed, the search space in (11) is of cardinality $|{{\cal K}}_{k,1,n}|=\binom{n}{k}^{2}$. On the other hand, the scan test ${\cal A}_{\mathsf{scan}}^{\mathsf{CSD}}$ for the consecutive setting is efficient because $|{\cal K}_{k,1,n}^{\mathsf{con}}|\leq n^{2}$. The following result provides sufficient conditions under which the risk of each of the above tests is asymptotically small.

As can be seen from Theorem 1, only the sum test performance barrier exhibits dependency on $m$. The scan test is, for both $\mathsf{SD}$ and $\mathsf{CSD}$, inherently independent of $m$. This makes sense because when summing all the elements of $\mathsf{X}$, as $m$ gets larger the mean (the “signal”) under the alternative hypothesis gets larger as well. On the other hand, since the scan test searches for a single planted submatrix, the number of planted submatrices does not play a role. One could argue that it might be beneficial to search for the $m$ planted submatrices in the scan test, however, as we show below, this is not needed, and the bounds above are asymptotically tight.

To present our lower bounds, we first recall that the optimal testing error probability is determined by the total variation distance between the distributions under the null and the alternative hypotheses as follows (see, e.g., [47, Lemma 2.1]),

The following result shows that under certain conditions the total variation between the null and alternative distributions is asymptotically small, and thus, there exists no test which can solve the above detection problems reliably.

Theorem 2 above shows that our upper bounds in Theorem 1 are tight up to poly-log factors. Indeed, item 1 in Theorem 2 complements Items 1-2 in Theorem 1, for the $\mathsf{SD}$ problem, while item 2 in Theorem 2 complements Item 3 in Theorem 1, for the $\mathsf{CSD}$ problem. In the sequel, we illustrate our results using phase diagrams that show the tradeoff between $k$ and $\lambda$ as a function of $n$. One evident and important observation here is that the statistical limit for the $\mathsf{CSD}$ problem is attained using an efficient test. Thus, there is no statistical computational gap in the detection problem in Definition 2, and accordingly, it is either statistically impossible to solve the detection problem or it can be solved in polynomial time. This is not the case for the $\mathsf{SD}$ problem. Note that both the efficient sum and the exhaustive scan tests are needed to attain the information-theoretic lower bound (up to poly-log factors). As discussed above, however, here the scan test is not efficient. We next give evidence that, based on the low-degree polynomial conjecture, efficient algorithms that run in polynomial-time do not exist in the regime where the scan test succeeds while the sum test fails.

Note that the problem in Definition 1 exhibits a gap in terms of what can be achieved by the proposed polynomial-time algorithm and the computationally expensive scan test algorithm. In particular, it can be seen that in the regime where $\frac{1}{\sqrt{k}}\ll\lambda\ll\frac{n}{mk^{2}}$, while the problem can be solved by an exhaustive search using the scan test, we do not have a polynomial-time algorithm. Next, we give evidence that, in fact, an efficient algorithm does not exist in this region. To that end, we start with a brief introduction to the method of low-degree polynomials.

The premise of this method is to take low-degree multivariate polynomials in the entries of the observations as a proxy for efficiently-computable functions. The ideas below were first developed in a sequence of works in the sum-of-squares optimization literature [10, 28, 31, 30].

In the following, we follow the notations and definitions of [28, 18]. Any distribution $\mathbb{P}_{{\cal H}_{0}}$ on $\Omega_{n}=\mathbb{R}^{n\times n}$ induces an inner product of measurable functions $f,g:\Omega_{n}\to\mathbb{R}$ given by $\left\langle f,g\right\rangle_{{\cal H}_{0}}=\mathbb{E}_{{\cal H}_{0}}[f(\mathsf{X})g(\mathsf{X})]$, and norm $\left\|f\right\|_{{\cal H}_{0}}=\left\langle f,f\right\rangle_{{\cal H}_{0}}^{1/2}$. We Let $L^{2}(\mathbb{P}_{{\cal H}_{0}})$ denote the Hilbert space consisting of functions $f$ for which $\left\|f\right\|_{{\cal H}_{0}}<\infty$, endowed with the above inner product and norm. In the computationally-unbounded case, the Neyman-Pearson lemma shows that the likelihood ratio test achieves the optimal tradeoff between $\mathsf{Type}$-$\mathsf{I}$ and $\mathsf{Type}$-$\mathsf{II}$ error probabilities. Furthermore, it is well-known that the same test optimally distinguishes $\mathbb{P}_{{\cal H}_{0}}$ from $\mathbb{P}_{{\cal H}_{1}}$ in the $L^{2}$ sense. Specifically, denoting by $\mathsf{L}_{n}\triangleq\mathbb{P}_{{\cal H}_{1}}/\mathbb{P}_{{\cal H}_{0}}$ the likelihood ratio, the second-moment method for contiguity (see, e.g., [18]) shows that if $\left\|\mathsf{L}_{n}\right\|_{{\cal H}_{0}}^{2}$ remains bounded as $n\to\infty$, then $\mathbb{P}_{{\cal H}_{1}}$ is contiguous to $\mathbb{P}_{{\cal H}_{0}}$. This implies that $\mathbb{P}_{{\cal H}_{1}}$ and $\mathbb{P}_{{\cal H}_{0}}$ are statistically indistinguishable, i.e., no test can have both $\mathsf{Type}$-$\mathsf{I}$ and $\mathsf{Type}$-$\mathsf{II}$ error probabilities tending to zero.

We now describe the low-degree method. The idea is to find the low-degree polynomial that best distinguishes $\mathbb{P}_{{\cal H}_{0}}$ from $\mathbb{P}_{{\cal H}_{1}}$ in the $L^{2}$ sense. To that end, we let ${\cal V}_{n,\leq\mathsf{D}}\subset L^{2}(\mathbb{P}_{{\cal H}_{0}})$ denote the linear subspace of polynomials $\Omega_{n}\to\mathbb{R}$ of degree at most $\mathsf{D}\in\mathbb{N}$. We further define ${\cal P}_{\leq\mathsf{D}}:L^{2}(\mathbb{P}_{{\cal H}_{0}})\to{\cal V}_{n,\leq\mathsf{D}}$ the orthogonal projection operator. Then, the $\mathsf{D}$-low-degree likelihood ratio $\mathsf{L}_{n}^{\leq\mathsf{D}}$ is the projection of a function $\mathsf{L}_{n}$ to the span of coordinate-degree-$\mathsf{D}$ functions, where the projection is orthogonal with respect to the inner product $\left\langle\cdot,\cdot\right\rangle_{{\cal H}_{0}}$. As discussed above, the likelihood ratio optimally distinguishes $\mathbb{P}_{{\cal H}_{0}}$ from $\mathbb{P}_{{\cal H}_{1}}$ in the $L^{2}$ sense. The next lemma shows that over the set of low-degree polynomials, the $\mathsf{D}$-low-degree likelihood ratio have exhibit the same property.

As was mentioned above, in the computationally-unbounded regime, an important property of the likelihood ratio is that if $\left\|\mathsf{L}_{n}\right\|_{{\cal H}_{0}}$ is bounded, then $\mathbb{P}_{{\cal H}_{0}}$ and $\mathbb{P}_{{\cal H}_{1}}$ are statistically indistinguishable. The following conjecture states that a computational analog of this property holds, with $\mathsf{L}_{n}^{\leq\mathsf{D}}$ playing the role of the likelihood ratio. In fact, it also postulates that polynomials of degree $\approx\log n$ are a proxy for polynomial-time algorithms. The conjecture below is based on [28, 31, 30], and [28, Conj. 2.2.4]. We give an informal statement of this conjecture, which appears in [18, Conj. 1.16]. For a precise statement, we refer the reader to [28, Conj. 2.2.4] and [18, Sec. 4].

In the sequel, we will rely on Conjecture 1 to give evidence for the statistical-computational gap observed for the problem in Definition 1 in the regime where $\frac{1}{\sqrt{k}}\ll\lambda\ll\frac{n}{mk^{2}}$. At this point we would like to mention [28, Hypothesis 2.1.5], which states a more general form of Conjecture 1 in the sense that it postulates that degree-$\mathsf{D}$ polynomials are a proxy for $n^{O(D)}$-time algorithms. Note that if $\left\|\mathsf{L}_{n}^{\leq\mathsf{D}}\right\|_{{\cal H}_{0}}=O(1)$, then we expect detection in time $\mathsf{T}(n)=e^{\mathsf{D}(n)}$ to be impossible.

Together with Conjecture 1, Theorem 3 implies that if we take degree-$\log n$ polynomials as a proxy for all efficient algorithms, our calculations predict that an $n^{O(\log n)}$ algorithm does not exist when $\frac{1}{\sqrt{k}}\ll\lambda\ll\frac{n}{mk^{2}}$. This is summarized in the following corollary.

These predictions agree precisely with the previously established statistical-computational tradeoffs in the previous subsections. A more explicit formula for the computational barrier which exhibits dependency on $\mathsf{D}$ and $\lambda$ can be deduced from the proof of Theorem 3; to keep the exposition simple we opted to present the refined result above.

We note that numerical and theoretical evidence for the existence of computational-statistical gaps were observed in other statistical models that are also inspired by cryo-EM, including heterogeneous multi-reference alignment [7, 49] and sparse multi-reference alignment [16].

Using Theorems 1–3 we are now in a position to draw the obtained phase diagrams for our detection problems. Specifically, treating $k$ and $\lambda$ as polynomials in $n$, i.e., $k=\Theta(n^{\beta})$ and $\lambda=\Theta(n^{-\alpha})$, for some $\alpha\in(0,1)$ and $\beta\in(0,1)$, we obtain the phase diagrams in Figure 2(a), for a fixed number of submatrices $m=O(1)$. Specifically,

1. 1.

   Computationally easy regime (blue region): there is a polynomial-time algorithm for the detection task when $\alpha<2\beta-1$.

2. 2.

   Computationally hard regime (red region): there is an inefficient algorithm for detection when $\alpha<\beta/2$ and $\alpha>2\beta-1$, but the problem is computationally hard (no polynomial-time algorithm exists) in the sense that the class of low-degree polynomials fails in this region.

3. 3.

   Statistically impossible regime: detection is statistically impossible when $\alpha>\frac{\beta}{2}\vee(2\beta-1)$.

When the number of submatrices grows with $n=\omega(1)$, we get different phase diagrams depending on its value. For example, if $m=\Theta(n^{1/4})$, we get Figure 2(b). Specifically,

1. 1.

   Computationally easy regime (blue region): there is a polynomial-time algorithm for the detection task when $\alpha<2\beta-\frac{3}{4}$.

2. 2.

   Computationally hard regime (red region): there is an inefficient algorithm for detection when $\alpha<\beta/2$ and $\alpha>2\beta-\frac{3}{4}$, but the problem is computationally hard (no polynomial-time algorithm exists) in the sense that the class of low-degree polynomials fails in this region.

3. 3.

   Statistically impossible regime: detection is statistically impossible when $\alpha>\frac{\beta}{2}\vee(2\beta-3/4)$.

Finally, for the consecutive problem, we get the phase diagram in Figure 3, independently of the value of $m$. Here, there are only two regions where the problem is either statistically impossible or easy to solve.

We start by presenting our upper bounds for both exact and correlated types of recovery for the consecutive problem in Definition 3. To that end, we propose the following recovery algorithm. It can be shown that the maximum-likelihood (ML) estimator, minimizing the error probability, is given by (see Subsection 4.4 for a complete derivation),

The computational complexity of the exhaustive search in (19) is of order $n^{2m}$. Thus, for $m=O(1)$, the ML estimator runs in polynomial time, and thus, is efficient. However, if $m=\omega(1)$ then the exhaustive search is not efficient anymore. Nonetheless, the following straightforward modification of (19) provably achieves the same asymptotic performance of the ML estimator above, and at the same time computationally efficient.

Before we present this algorithm, we make a simplifying technical assumption on the possible set of planted submatrices, and then explain how this assumption can be removed. We assume that each pair of submatrices in the underlying planted submatrices $\mathsf{K}^{\star}$ are at least $k$ columns and rows far way. In other words, there are at least $k$ columns and $k$ rows separating any pair of submatrices in $\mathsf{K}^{\star}$. Similar assumptions are frequently taken when analyzing statistical models inspired by cryo-EM, see, for example [6]. We will refer to the above as the separation assumption.

Our recovery algorithm works as follows: in the $\ell\in[m]$ step, we find the ML estimate of a single submatrix using,

where $\mathsf{X}^{(\ell)}$ is defined recursively as follows: $\mathsf{X}^{(1)}\triangleq\mathsf{X}$, and for $\ell\geq 2$,

where $\mathsf{E}(\hat{\mathsf{K}}_{\ell-1})$ is an $n\times n$ matrix such that $[\mathsf{E}(\hat{\mathsf{K}}_{\ell-1})]_{ij}=-\infty$, for $(i,j)\in\hat{\mathsf{K}}_{\ell-1}$, and $[\mathsf{E}(\hat{\mathsf{K}}_{\ell-1})]_{ij}=1$, otherwise. To wit, in each step of the algorithm we “peel” the set of estimated indices (or, estimated submatrices) in previous steps from the search space. This is done by setting the corresponding entries of $\mathsf{X}$ to $-\infty$ so that the sum in (20) will not be maximized by previously chosen sets of indices. We denote by $\hat{\mathsf{K}}_{\mathsf{peel}}(\mathsf{X})=\{\hat{\mathsf{K}}_{\ell}\}_{\ell=1}^{m}$ the output of the above algorithm.

We have the following result.

The following result shows that under certain conditions, exact and correlated recoveries are impossible.

Thus, similarly to the detection problem, the consecutive recovery problem is either statistically impossible or easy to solve. The corresponding phase diagram for exact and correlated types of recoveries is given in Figure 4. Roughly speaking, exact recovery is possible if $\lambda=\omega(k^{-1/2})$ and impossible if $\lambda=o(k^{-1/2})$. Correlated recovery is possible if $\lambda=\omega(k^{-1})$ and impossible if $\lambda=o(k^{-1})$.

A few remarks are in order. First, note that there is a gap between detection and exact recovery; the barrier for $\lambda$ for the former is at $k^{-1}$, while for the latter at $k^{-1/2}$. In the context of cryo-EM, this indicates a gap between the ability to detect the existence of particle images in the data set, and the ability to perform successful particle picking (exact recovery). Recently, new computational methods were devised to elucidate molecular structures without particle picking, thus bypassing the limit of exact recovery, allowing constructing structures in very low SNR environments, e.g.,  [6, 36, 38]. This in turn opens the door to recovering small molecular structures that induce low SNR [29]. Second, there is no gap between detection and correlated recovery, and these different tasks are asymptotically statistically the same. The same gap exists between correlated and exact recoveries, implying that exact recovery is strictly harder than correlated recovery.