Detecting positive correlations in a multivariate sample

Introduction

We are interested in testing whether the population covariance matrix is the identity matrix, or not, so the null hypothesis is

This testing problem is well studied in the classical regime where the dimension $n$ is fixed and the sample size $m$ increases to infinity. See, for example, Muirhead , Section 8.4, where the generalized likelihood ratio test (GLRT) – against the alternative hypothesis $H_{1}\colon\ \sigma_{i,j}\neq 0$ for some $i\neq j$ – is studied in detail, as well as some unbiased variant. When the dimension is large (i.e., $n\to\infty$), the GLRT may be degenerate. This is discussed in detail in Ledoit and Wolf , where other tests – including a new one – are examined for consistency in this high-dimensional regime. Their ideas are further explored in Srivastava , Fisher and Chen, Zhang and Zhong . All these tests are based on symmetric polynomials of the sample correlation coefficients. We discuss this in Section 3.1. These tests are shown to be consistent when $n/m\to c\in(0,\infty)$ under additional mild conditions.

While these papers focus on general consistency, our focus is on alternatives where the covariance matrix is sparse, meaning that even under the alternative hypothesis, only a few variables are substantially correlated. This sparse setting has been investigated in the last few years, with recent work on the estimation of sparse covariance matrices, see El Karoui , Bickel and Levina and Cai, Zhang and Zhou . To our knowledge, testing for sparse correlation structures in a multivariate sample has not been considered in detail before. This is what we study here.

We introduce sparse models of correlation matrices to test against. Though many more models are possible, we choose a few emblematic examples that are of interest in a much wider sense within the literature on sparse covariance estimation. In all cases, the null hypothesis is that the observed vector has identity covariance matrix. For the alternative hypothesis, we consider the following prototypical examples:

Block model. The covariance under the alternative hypothesis is the identity matrix except for a $k\times k$ block on the diagonal. Formally, given $\rho>0$, we assume here that there is a subset of indices of the form $S=\{i,\dots,i+k-1\}$ modulo $n$ – for aesthetic reasons – such that $\sigma_{i,j}\geq\rho$ if $i,j\in S,i\neq j$. The set $S$ is called the anomalous set.

Clique model. This model is defined as the block model with the possible anomalous set $S$ ranging over all the subsets of indices of size $k$.

Perfect matching model. Suppose $n$ is a perfect square with $n=k^{2}$. Here the components of the observed vector $X$ correspond to edges of the complete bipartite graph on $2k$ vertices. The alternative hypothesis is that the bipartite graph has a perfect matching such that $\sigma_{i,j}\geq\rho$ for all $i,j\in S,i\neq j$ where $S$ is the anomalous set of indices corresponding to the edges of the perfect matching.

The block model is closely related to the models used in Cai, Zhang and Zhou to obtain bounds on the minimax risk of estimating sparse matrices. Roughly speaking, Cai, Zhang and Zhou use the block model with $S=\{1,\dots,k\}$ and place nonzero entries in a (carefully designed) fashion within that block. The fraction of nonzero entries within the block is about one-half. We could also assume that only a fraction of the entries in the block are nonzero and it would only change constants later on. More importantly, to make the detection problem interesting, we need to consider all possible blocks. Note that the block model is parametric. The clique model is a natural generalization of the block model leading to a nonparametric model. The perfect matching model gives an example of a class of sets with a more intricate combinatorial structure which our approach is able to deal with.

2 Tests and their risks

and indeed all lower bounds derived in this paper start with this inequality. We will derive upper and lower bounds for the minimax risk,

The lower bounds will be obtained by putting a prior on model $\mathcal{C}$ and obtaining a lower bound on the corresponding Bayesian risk which never exceeds the worst-case risk. In all cases, we draw the set $S$ uniformly at random within the class $\mathcal{C}$. The upper bounds are obtained by studying the performance of specific tests.

We focus on the case where the dimension $n$ and the sample size $m$ are both large. Of course, such asymptotic statements only make sense if we define sequences of integers $m=m_{n}$, $k=k_{n}$, positive reals $\rho=\rho_{n}$, and classes $\mathcal{C}=\mathcal{C}_{n}$. This dependency in $n$ will be left implicit. In this asymptotic setting, we say that reliable detection is possible (resp. impossible) if $R_{*}^{\max}\to 0$ (resp. $\to 1$) as $n\to\infty$. Also we say that a sequence of tests $(f_{n})$ is asymptotically powerful (resp. powerless) if $R^{\max}(f_{n})\to 0$ (resp. $\to 1$).

3 A preview of results for the clique model

Among the models we consider, the clique model is perhaps the most compelling because of its relevance in applications and its complexity. Also, for a given value of $k$, the clique model is the richest possible and therefore for any given $\rho,n,m$, $R_{*}^{\max}$ is larger than for any other model. This makes the clique model an important benchmark.

Here we summarize our main findings for this special class. We discover various types of behavior in distinct ranges of the parameters $n,m,k,\rho$. Roughly speaking, and ignoring logarithmic factors, we arrive at the following conclusions. Two tests are competing for near-optimality. The first one is a “global” test akin to the classical test Muirhead , Section 8.4, and the refinements in Chen, Zhang and Zhong and Ledoit and Wolf . The second is a “local” test reminiscent to the generalized likelihood ratio test. The latter dominates the former when

is small, corresponding to smaller values of $k$.

The “local” test that achieves near-optimal behavior in a large range of the parameters is a scan statistic that requires the computation of a maximum over all ${n\choose k}$ subsets of components of size $k$. In its naive implementation, this test is computationally intractable, unless $k$ is very small. We also believe that computing this test is a fundamentally hard computational problem. We do not have a rigorous argument to prove such a hardness result but it is worth pointing out that the problem is quite similar, in spirit, to the notoriously difficult hidden clique problem, see Alon, Krivelevich and Sudakov .

What performance can we achieve with limited computational power? Such questions of trade-off between statistical performance and computational complexity are at the heart of high-dimensional statistics and machine learning. We probe this question and describe a family of tests that balances detection performance and computational complexity.

3.2 An application in the study of random geometric graphs

In Section 7, we apply the lower bound for the optimal risk in the clique model in a perhaps unexpected context and derive a new lower bound for the clique number of a high-dimensional random geometric graph. The setup is as follows.

4 More related work

As mentioned before, the literature on sparse covariance estimation has become quite extensive. In spite of this surge of interest in sparse high-dimensional models, not much has been done in terms of detection of correlations. We note the work of Verzelen and Villers , who consider the task of testing a given dependency structure. Our objective here is admittedly more modest and a more closely related is our own paper , which focuses entirely on the case where the sample size is equal to one (i.e., $m=1$). Our results here are seen to extend those in the one-sample case, with the regimes now partitioned according to the sample size.

Note that our work is different from Butucea and Ingster where the task is the detection of a submatrix with higher per-coordinate mean in a large matrix with i.i.d. Gaussian entries, which is more closely related to the literature on the detection of sparse nonzero entries in the mean of a random vector. Our work has parallels with that literature which, for the clique model, focuses on the “detection-of-means” problem (see Jin , Ingster , Baraud , Donoho and Jin , Hall and Jin , Arias-Castro, Candès, Helgason and Zeitouni , Addario-Berry, Broutin, Devroye and Lugosi ) defined as follows: Under the null hypothesis, the vectors $X_{t}$ are i.i.d. standard normal, while under the alternative hypothesis, there is a subset $S\subset\{1,\ldots,n\}$ in some class $\mathcal{C}$ of interest such that the $X_{t}$ are i.i.d. normal with mean $(\mu_{1},\dots,\mu_{n})^{T}$ and identity covariance, where $\mu_{i}\geq\mu$ for $i\in S$ and $\mu_{i}=0$ for $i\notin S$. Thus, $\mu>0$ is the minimum (per-coordinate) signal amplitude. Of course, one immediately reduces by sufficiency to the case $m=1$ by averaging over the sample. This explains why the literature focuses on the case $m=1$. The connection between the detection-of-means problem with the correlation detection problem studied here was detailed (for $m=1$) in our previous paper , where $\rho$ was found to correspond to $\mu^{2}$. The connection is based on the following simple representation of equi-correlated normal random variables.

Let $X_{1},\dots,X_{k}$ be standard normal random variables with $\operatorname{Cov}(X_{i},\allowbreak X_{j})=\rho>0$ for $i\neq j$. Then there are independent standard normal random variables $V,Y_{1},\dots,Y_{k}$ such that $X_{i}=\sqrt{\rho}V+\sqrt{1-\rho}Y_{i}$ for all $i$.

Thus, given $V$, the problem becomes that of detecting a subset of variables – here implicitly assumed to be indexed by $S=\{1,\ldots,k\}$ – with nonzero mean (equal to $\sqrt{\rho}V$) and with a variance equal to $1-\rho$ (instead of $1$). This representation was used in to obtain a general lower bound that seemed otherwise out of reach of more standard methods based on the second moment of the likelihood ratio.

This connection with the detection-of-means problem also applies in the case where $m>1$, but with a twist. Indeed, when detecting correlations one does not average the vectors $X_{t}$ but their covariances. So a simple reduction to the case $m=1$ does not apply. However, one may still apply the representation result Lemma 1 to each observation vector $X_{t}$, yielding $V_{t}$’s and $Y_{t,i}$’s that are independent standard normal random variables. By conditioning on $V_{1},\dots,V_{m}$, the problem becomes equivalent to detecting a subset of variables with means $\sqrt{\rho}V_{t},t=1,\dots,m$. What makes the situation more complex is that the signs of the $V_{t}$’s are random. Our approach to finding a general lower bound is based on this representation without which more standard methods seem to fail. The general lower bound, which is the key technical result of this paper, is given in Theorem 2 below.

5 Contribution and content of the paper

We obtain a general lower bound in Section 2 akin to, but not a straightforward extension of, the lower bound we obtained in . We then study a number of tests that are near optimal in the sense that they come close to achieving the detection lower bound for various models. This is done in Section 3. We then specialize these general results in Sections 4, 5 and 6, to the three models described in Section 1.1. We also discuss computational issues, particularly in the clique model. In Section 7, we apply our general lower bound to the problem of studying the size of the clique number of a random geometric graph on a high-dimensional sphere. We close the paper with a discussion in Section 8 of possible extensions and challenges.

Lower bounds

In this section, we derive a general lower bound for the minimax risk $R_{*}^{\max}$. As mentioned in Section 1.2, the first step is to restrict the supremum in the definition of $R^{\max}(f)$ to covariance matrices in which all the nonzero entries are equal to $\rho>0$ and then lower bound the maximum by an average. In particular, we have $R_{*}^{\max}\geq R^{*}$ where $R^{*}=\inf_{f}R(f)$ and

Note that $R^{*}$ is just the Bayes risk for the uniform prior on the models $S\in\mathcal{C}$. It is well known that the test $f^{*}$ that achieves the infimum (i.e., $R(f^{*})=R^{*}$) is the likelihood ratio test with critical value 1, and there is a whole machinery that can be used to bound that risk from below.

The following lower bound has a similar flavor as the main result in our previous work . In particular, we make appear some moment of $Z$, a random variable that represents the size of the overlap of two index sets taken at random from the class $\mathcal{C}$. A straightforward adaptation of the arguments we used in leads to a lower bound in terms of the moment generating function of $Z$. Here, unfortunately, this quantity is too large to obtain sharp results in most regimes. The key contribution of the following result is to replace the exponential function by the hyperbolic cosine. This allows us to derive much sharper results, essentially because around $0$ one has $\exp(x)-1\sim x$ while $\cosh(x)-1\sim\frac{x^{2}}{2}$.

For any class $\mathcal{C}$, any $\rho\in(0,1)$, and any $a\geq\sqrt{3}$,

and where $\chi_{m}^{2}$ has chi-squared distribution with $m$ degrees of freedom, and $Z=|S\cap S^{\prime}|$ with $S,S^{\prime}$ i.i.d. uniform from $\mathcal{C}$.

where, for any positive integer $m$, $[m]:=\{1,\dots,m\}$, and $(Y_{t,i})_{i\in[n],t\in[m]},(V_{t})_{t\in[m]}$ are i.i.d. standard normal random variables.

Therefore, using the Cauchy–Schwarz inequality,

where $\varepsilon,\varepsilon^{\prime}$ are i.i.d. Rademacher vectors and $S,S^{\prime}$ are i.i.d. uniform in the class $\mathcal{C}$. We have

Let $Z=|S\cap S^{\prime}|$. We see that $H_{1}(X),H_{2}(X),H_{3}(X)$ are independent of each other under the null hypothesis with

For the latter, we used the fact that $\varepsilon_{t}^{2}={\varepsilon_{t}^{\prime}}^{2}=1$, to get

where the last line comes from a simple change of variables. Hence,

Let $\xi=\xi_{1}=\rho/(1-\rho^{2})$. Since $(\varepsilon_{t}\varepsilon^{\prime}_{t}\colon\ t=1,\dots,m)$ are i.i.d. Rademacher, we have

Holding $Z\geq 1$ fixed, we maximize this over $\|u\|^{2}=\sum_{t}u_{t}^{2}\leq a^{2}m$ using Lagrangian multipliers and checking the Karush–Kuhn–Tucker conditions, finding that at a local maximum all $u_{t}^{2}$ must be equal. Hence,

where the last equality comes from the fact that the function $h_{\rho}(c):=\cosh(c)\exp(-\rho c)$ is decreasing on $(0,\rho)$ and increasing on $(\rho,\infty)$, so that its maximum over $[0,\xi_{a}Z]$ is either at $c=0$ or $c=\xi_{a}Z$. Straightforward calculations lead to

Since $g(c)>1-\frac{1}{c}\log 2$, the maximum of $h_{\rho}(c)$ over $c\in[0,\xi_{a}Z]$ is at $c=\xi_{a}Z$ when

Since we consider $Z\geq 1$, the last inequality is true if $\rho\geq 1/2$ and $a^{2}\geq 3\log(2)$. This inequality is far off when $\rho$ is small, so we need to derive another bound. Noting that $g$ is seen to be strictly increasing on $(0,\infty)$ with range $(0,1)$, $w(\rho):=g^{-1}(\rho)$ is well defined and, as a function of $\rho$, is infinitely differentiable and strictly increasing. Elementary calculations show that

When $Z\geq 1$, the latter is true when $\rho\leq 1/2$ and $a^{2}\geq 3/2$. Hence, given that $a^{2}\geq 3$ by assumption, the maximum in (3) at $c=\xi_{a}Z$.

where in the first line we used $a^{2}\geq 1$ and in the second line the fact that $s+\frac{1-s}{2}\log(1-s)\geq 0$ for all $s\in(0,1)$. With this, we conclude. ∎

In Sections 4, 5, and 6, we specialize Theorem 2 to the different models we described in Section 1.1.

Tests

In this section, we introduce and briefly discuss two natural tests that will be seen to perform near optimally in various regimes of the parameters. This optimality property will be established in Sections 4, 5 and 6, by comparing simple performance bounds with the implications of Theorem 2.

The first test, that we call “squared-sum test”, is based on a global test statistic that does not take the class $\mathcal{C}$ into account at all.

The second test, a “localized” squared-sum test, is based on a simple scan statistic. It may also be interpreted as a simplified version of the generalized likelihood ratio test.

As we will see, one of the two tests above always has a near-optimal performance in all three specific classes we discuss. Thus, the story is essentially complete for the point of view of detection performance. Unfortunately, when the class $\mathcal{C}$ is large – as in the clique model –, the localized squared-sum test is computationally unfeasible, at least in its naive implementation. We discuss two possible substitutes. The first one is a simple “maximum correlation test” that turns out to be nearly optimal for very small values of $k$. In Section 4.5, we discuss another test in the context of the clique model that is both near-optimal and computationally feasible when the sample size is at most logarithmic in the dimension $n$. In Section 4.6, a conceptually different computationally efficient alternative is discussed.

All performance bounds derived below are in terms of the average correlation

Let $\Sigma$ denote the covariance matrix of the distribution of $X_{1}$. When the alternative hypothesis is simply $H_{1}\colon\ \Sigma\neq I$ (where $I$ is the identity matrix), without any sign restriction on the entries of $\Sigma$, one of the simplest tests is that of Nagao , which is based on the Frobenius norm of the difference between the sample covariance matrix $\widehat{\Sigma}$ and the identity matrix $I$. Nagao’s test is based on the test statistic

We also refer to Schott , who (like us) assumes that the variables have unit variance under the alternative hypothesis. Ledoit and Wolf show that this test is not always consistent against fixed alternatives. They, and others including Srivastava , Fisher and Chen, Zhang and Zhong , suggest variants based on consistent estimates for the Frobenius norm $\operatorname{tr}[(\Sigma-I)^{2}]$.

Given that we know that the variances are equal to 1 and the correlations are non-negative under the alternative hypothesis, it is more natural to consider the test that rejects for large

values of $\sum_{i<j}\widehat{\sigma}_{i,j}$, where $\widehat{\Sigma}=(\widehat{\sigma}_{i,j})$. For simplicity, we consider instead the squared-sum test that rejects for large values of the test statistic

The two tests are thus closely related. In fact, one may easily check that they have similar asymptotic power properties. Our preference for the second test is only for convenience.

The following result gives a simple characterization of the performance of the squared-sum test. Since the test does not use information about the class $\mathcal{C}$, its minimax risk does not depend on the model either.

Using the assumptions on $a$ and $b$, we have

and therefore the test with critical value $n(m+a\sqrt{m})$ is asymptotically powerful.

Suppose that $b\sqrt{m}\to 0$. We still have that $Y/n\sim\chi^{2}_{m}$ under $H_{0}$ while $Y/n\sim(1+b)\chi^{2}_{m}$ under $H_{1}$. If $m$ is fixed, $Y/n$ is asymptotically $\chi^{2}_{m}$ under the alternative hypothesis since $b\to 0$ in this case. If $m\to\infty$, $(Y/n-m)/\sqrt{2m}$ is asymptotically standard normal under both the null and the alternative hypotheses, since under $H_{1}$,

2 A localized squared-sum test

When $k$ is smaller, global tests such as the squared-sum test are not very powerful. The generalized likelihood ratio test “scans” over all subsets $S$ in the class $\mathcal{C}$. Instead of studying the generalized likelihood ratio test, we consider a localized version of the squared-sum test that has similar power and is a little easier to analyze. The localized squared-sum test rejects the null hypothesis for large values of the test statistic

Other consequences of this proposition will be discussed in the next sections. {pf*}Proof of Proposition 2 Observe that under the null hypothesis $Y_{S}\sim k\chi^{2}_{m}$ for all $S\in\mathcal{C}$. By a simple Chernoff bound for the chi-square distribution, for all $b>1$,

where $H(b):=b-1-\log b$ for $b>1$. Hence, by the union bound,

When $(\log|\mathcal{C}|)/m\to 0$, using the fact that $H(b)\sim(b-1)^{2}/2$ when $b\to 1$, we see that the right-hand side in (8) tends to zero when $b\geq 1+\sqrt{5\log|\mathcal{C}|/m}$. When $(\log|\mathcal{C}|)/m\to\infty$, using the fact that $H(b)\sim b$ when $b\to\infty$, so the right-hand side in (8) tends to zero when $b\geq 3\log|\mathcal{C}|/m$.

The case when $\log|\mathcal{C}|/m\asymp 1$ can be dealt with in the same way, yielding that, with a proper choice of threshold, the localized squared-sum test is asymptotically powerful when

3 Maximum correlation test

Finally, we mention the possibly simplest test that one would think of when confronted with testing $H_{0}$ in the sparse regime. This is the test that rejects for large values of the maximum pairwise empirical correlation

In fact, this test does have some power in the sparse regime, and is actually near-optimal when $k$ is fixed as the following result shows. However, one cannot expect a good performance of this test for large values of $k$. An advantage of this test is that it may be computed efficiently in a straightforward manner.

The maximum correlation test that rejects $H_{0}$ when $Y_{\max}>\sqrt{5m\log n}$ is asymptotically powerful when

From this, the result follows immediately.

Clique model

In this section, we discuss the implications of the general results of the previous sections for the clique model. We derive a lower bound based on Theorem 2 in various ranges of the parameters and compare it with the performance bounds for the squared-sum test and the scan statistics-based test considered in Section 3. We also propose a goodness-of-fit test for the case where $\rho\to 1$, and consider two alternative tests that take computational considerations into account. A digest is provided at the end of the section.

In order to apply Theorem 2, note that in the clique model, $Z$ has hypergeometric distribution with parameters $(n,k,k)$, which is stochastically bounded by the binomial distribution with parameters $(k,p)$, where $p:=k/(n-k)$. In particular, for all $\xi\geq 0$,

where $H(t):=t(\log t-1)+1$. Note that $H(t)\sim t^{2}/2$ when $t\to 0$ and $H(t)\sim t\log t$ when $t\to\infty$.

Case 1: large $k$. Suppose that $k$ is so large and $\rho$ is so small that

We first note that these conditions imply that $\rho\to 0$. Let $\zeta=\rho\sqrt{m}k^{2}/n$ and choose $a,b\to\infty$ such that $a^{2}b\zeta\to 0$. When $Z\leq bk^{2}/n$, we use the fact that $\rho a^{2}bZ\to 0$ and $\cosh(x)\leq 1+x^{2}$ when $x\in(0,1)$ to get that for all sufficiently large $n$,

We now show that, if, in addition to (14), we have either $\rho m\to 0$ or $\rho^{2}mk\to 0$, then

This implies that reliable detection is impossible in this range of the parameters.

We choose $a$ such that $a^{2}\rho m\to 0$. We use the bound $\cosh(x)\leq\exp(x)$ and (13), to get

uniformly over $z>bk^{2}/n$. Hence, eventually,

We may assume that $k\leq m$ for otherwise $\rho m\to 0$, which we already covered. We choose $a$ such that $a^{2}\rho\sqrt{mk}\to 0$ – which implies in particular that $a^{2}\rho k\to 0$. We use the bounds $\cosh(x)\leq 1+x^{2}$ for $x\in$, the fact that $Z\leq k$ – since $Z=|S\cap S^{\prime}|$ with $|S|=|S^{\prime}|=k$ – and the fact that $a^{2}\rho k\to 0$, to get

Let $\zeta=\rho\sqrt{mk/\log(n/k^{2})}$ and choose $a\to\infty$ such that $a\zeta\to 0$ and $a^{2}\rho k\to 0$. The latter is possible because (17) implies that $\rho k\to 0$. Then, as in Case 1(b),

and again, reliable detection is impossible by Theorem 2.

Hence, there is some $\varepsilon>0$ fixed such that $\rho<1-\varepsilon$. Let $\zeta=\rho m/\log(n/k^{2})$ and choose $a\to\infty$ such that $a^{2}\zeta\to 0$. We use the fact that $\cosh(x)\leq\exp(x)$, and use the same bound on the moment generating function of $Z$, to get (eventually)

implying that reliable detection is impossible.

This is the only situation where we bound

The discussion of these various regimes leads to the following.

In the clique model, under either (14) with (15) or (16), (17), (18), or (19), $R^{*}\to 1$.

2 Localized squared-sum test

Next, we take a closer look at the performance of the localized squared-sum test for the clique model. In this case, we have $|\mathcal{C}|={n\choose k}$ so $\log|\mathcal{C}|\sim k\log(n/k)$. Plugging this into (10), we see that the local squared-sum test is asymptotically powerful when

and the constant $A$ is large enough. Based on this and Corollary 3, we conclude that the test is near-optimal in regimes (17) and (18), though only up to a logarithmic factor if $k^{2}/n\to 0$ slower than any power of $n$. It is also near-optimal up to a logarithmic factor in regime (14) when neither (15) nor (16) is satisfied.

However, we do not have such a guarantee in the regime (14) (with either (15) or (16)). In this range of parameters, it is the squared-sum test that yields an optimal performance up to a logarithmic factor. Also, comparing Proposition 1 and Proposition 2, we see that the local test dominates when $\max(1,(k/m)^{1/2})k^{3/2}/n$ tends to zero faster than $1/\log(n/k)$.

We make some progress in this direction in two ways. In Section 4.5, we suggest a test that has good performance and that is efficiently computable if $m$ is only logarithmic in $n$. In Section 4.6, inspired by recent work of Berthet and Rigollet , we consider a convex relaxation of the problem following d’Aspremont, El Ghaoui, Jordan and Lanckriet .

3 The case of ρ\rho constant

The situation changes dramatically when the sample size $m$ becomes at least logarithmic in the dimension $n$. Indeed, even for $k=2$, both the localized squared-sum test and the maximum correlation test have a vanishing risk for any constant value of $\rho$ when $\log(n)/m\to 0$. This reveals an interesting “phase transition” occurring when the sample size is about logarithmic in the dimension.

4 The case of ρ\rho tending to 1

The regime in (19) does not have a match in either the squared-sum test or the localized squared-sum test. It is instead met by a goodness-of-fit test which is a variant of test proposed and analyzed in our previous work in the same regime with $m=1$, although the construction here is slightly different.

Here we assume that $\rho\to 1$ so fast that

Let $\zeta$ denote the last term tending to zero in (4.4), and choose $a\to\infty$ such that $a\max(\zeta,\allowbreak 1/\log n)\to 0$.

The test we propose is based on the idea that the variables that are positively correlated are closer together than the other variables that are independent of each other. Take $\eta\to 0$ such that

Under the null hypothesis, $B_{s}\sim\operatorname{Bin}(n,p_{s})$ where

A simple combination of Bernstein’s inequality and the union bound gives

Under the alternative hypothesis where $S$ is anomalous, we use the representation (2), to get

by Markov’s inequality. Hence, by another application of Markov’s inequality,

Finally, when $\eta=2a(1-\rho)^{1/2}\sqrt{m}\geq 2(\log(n)/n)^{1/m}$, we have

5 Balancing detection ability and running time

Given the often enormous size of data sets that statisticians need to handle as an every-day practice, it is of great interest to design computationally efficient, yet near-optimal tests. In the case of the clique model, this is a highly non-trivial task, because the class $\mathcal{C}$ has size exponential in $k$ and computing the localized squared-sum test (or other versions of the generalized likelihood ratio test and scan statistics) involves a non-trivial optimization problem over all ${n\choose k}$ elements of $\mathcal{C}$. In fact, often it seems that small testing risk and computational efficiency are contradicting terms. In this section, we show that in at least one non-trivial instance, it is possible to design a computationally efficient (i.e., computable in time quadratic in $n$) test that has near optimal risk.

This is the case when the sample size $m$ is (at most) logarithmic in $n$ and $k\sim n^{a}$ for some $a\in(0,1)$. (Recall from Section 4.3 that this is a quite interesting range of parameters.)

is asymptotically powerful in the clique model when

Under the alternative hypothesis where $S$ is anomalous, we have

where $Z_{(1)}\leq\cdots\leq Z_{(m)}$ are the ordered values of

In the regime of (18) with $m\sim C\log n$, we see that the test is optimal up to a constant factor in $\rho$ when $k\sim n^{a}$ for some $a<1/2$. In this range of parameters, it seems hopeless to compute (or even approximate) the local squared-sum test.

However, when $m$ is much larger than logarithmic in $n$, this test also requires super-polynomial computational time and therefore it is not useful in practice. In such cases, one may have to resort to sub-optimal tests such as the maximum correlation test described in Section 3.3. It is an important and difficult challenge to find out the possibilities and limitations of powerful detection taking computational constraints into account.

6 A convex relaxation

In parallel to our work, Berthet and Rigollet study a related problem of detecting a sparse principal component. The setting there is the same, except for the alternative hypothesis, where the covariance matrix is of the form $\Sigma=I+\theta vv^{\top}$, with $\theta>0$ and $v$ a unit vector with at most $k$ non-zero components. They study a test based on $\lambda_{k}^{\max}(\widehat{\Sigma})$, the largest $k$-sparse eigenvalue of $\widehat{\Sigma}$, defined as

where $A_{S}$ denotes the principal submatrix of $A$ indexed by $S$ and $\lambda^{\max}(A)$ the largest eigenvalue of $A$. This test is, in fact, intimately related to our localized squared-sum test, as we shall see in the analysis below. Berthet and Rigollet prove that this test – with a proper choice of critical value – is near-optimal. However, just like our localized squared-sum test, the test of Berthet and Rigollet is also computationally unfeasible due to the maximization over ${n\choose k}$ sets. For computational reasons, turn to the convex relaxation of d’Aspremont, El Ghaoui, Jordan and Lanckriet , for which they also establish a performance bound.

Berthet and Rigollet show that, when $n,m\to\infty$, with high probability under the null hypothesis,

for a universal constant $C$. Under the alternative hypothesis where $S$ is anomalous, we have

which matches the performance of the localized squared-sum test (20) up to a multiplicative constant.

The semidefinite relaxation of d’Aspremont, El Ghaoui, Jordan and Lanckriet for $\lambda_{k}^{\max}$ is

for a universal constant $C$, while under the alternative hypothesis,

Hence, the test based on the statistic $\operatorname{SDP}_{k}(\widehat{\Sigma})$ is asymptotically powerful when

This rate matches (23) when $(k/m)\log(n/k)\to\infty$, and is otherwise comparable to what the maximum correlation test achieves (11). Thus, the relaxed test of Berthet and Rigollet is computationally efficient and near-optimal when the sample size is of smaller order that $k\log(n/k)$. Note that this allows one to handle larger values of $m$ than for the test introduced in Section 4.5 where $m$ had to be at most a constant multiple of $\log n$.

Interestingly, Berthet and Rigollet also show that their analysis of the relaxed test is optimal in the following sense. If one can improve the rate given in (24) for the relaxed test, then one obtains an algorithm that improves by an order of magnitude upon known results for the hidden clique problem, see Alon, Krivelevich and Sudakov . We refer to for a more precise statement.

7 Digest

In this section, we briefly summarize our findings for the case of the clique model in simplified regimes. We consider combinations of the following settings:

The sparse regime corresponds to $k\ll n^{1/2-\epsilon}$ for some $\epsilon>0$. The non-sparse regime corresponds to $k\gg\sqrt{n}$.

The ultra-high dimensional setting corresponds to $m\ll k\log n$ , while the (potentially) high dimensional setting corresponds to $m\gg k\log n$.

Our results can be summarized as follows:

In the sparse high dimensional setting, detection is impossible when $\rho\ll\sqrt{(\frac{\log n}{mk})}$

(see (17)). On the other hand when $\rho\gg\sqrt{(\frac{\log n}{mk})}$, the localized squared-sum test is asymptotically powerful. Furthermore in the extremely sparse case where $k$ is a constant, the same performance is achieved by the maximum correlation test.

In the sparse ultra-high dimensional setting, detection is impossible when $\rho\ll\min(\frac{\log n}{m},1)$ (see (18)). This rate is matched by the localized squared-sum test, and by the computationally efficient test of Berthet and Rigollet based on $\operatorname{SDP}_{k}$ (see Section 4.6). Furthermore, when $m\sim C\log n$ the test of Section 4.5 can also be computed in polynomial time (in $n$) and is asymptotically powerful for $\rho\gg\frac{\log n}{m}+\frac{1}{k}$. When $m\ll\log n$, detection is impossible when $(1-\rho)^{1/2}(\frac{n}{k^{2}})^{1/m}\to\infty$ (see (19)) and the goodness-of-fit test of Section 4.4 matches that bound up to a sub-logarithmic factor.

In the non-sparse regime, the squared sum test is asymptotically powerful for $\rho\gg\frac{n}{k^{2}\sqrt{m}}$. This rate is optimal (see (14)) if either $m$ or $k$ is not too large (that is either (15) or (16) is satisfied). In case neither (15) nor (16) is satisfied, the localized squared-sum test is asymptotically powerful. Note that in the non-sparse case we can differentiate two regimes for the value of $k$. If $n^{1/2}\ll k\ll n^{2/3}$, then condition (16) is more demanding than the last part of condition (14), while for $n^{2/3}\ll k\ll n$ it is the other way around. As a referee kindly pointed out, in the former case one may further tighten the lower bound and recover missing logarithmic factors by a more careful bounding of the moment generating function of $Z^{2}$.

Block model

Next, we discuss the consequences of our main results for the block model which serves as a prototypical example of a “small” or “parametric” class. We focus on the case where $\rho$ is bounded away from 1. Specifically, we assume that $\rho\leq\rho_{0}<1$, and define $C_{0}=(1-\rho_{0}^{2})^{-1}$.

We distinguish between two main regimes and we show that $R^{*}\to 1$ in both cases.

Let $\zeta=\rho k\sqrt{m/\log(n/k)}$ and choose $a\to\infty$ such that $a^{2}\zeta\to 0$ and $a^{2}\rho k\to 0$. The latter is possible because (25) implies that $\rho k\to 0$. We use the bound $\cosh(x)\leq 1+x^{2}$ for $x\in(0,1)$, to get

by our assumptions. Theorem 1 now implies that reliable detection is impossible in this range of the parameters.

Let $\zeta=\rho km/\log(n/k)$ and choose $a\to\infty$ such that $a^{2}\zeta\to 0$. We use the bound $\cosh(x)\leq\exp(x)$ to get

In the block model, under either (25) or (26), $R^{*}\to 1$.

In view of Corollary 4, the squared-sum test is near-optimal for the block model only when $k\asymp n$. However, the localized squared-sum test has a much better performance. We have $|\mathcal{C}|=n$, and plugging this into (10), we see that the localized squared-sum test is asymptotically powerful when

for a large enough constant $A>0$. With Corollary 3, we conclude that the test is near-optimal except in the case where $k/n\to 0$ slower than any negative power of $n$, where the test is optimal up to a logarithmic factor.

Perfect matching model

Here we work out the corollaries of our main results for the perfect matching model. This model illustrates how one may proceed when the model in question has a non-trivial combinatorial structure. In order to use Theorem 2, one needs to use the specific properties of the class. We focus on the case where $\rho$ is bounded away from 1. Specifically, we assume that $\rho\leq\rho_{0}<1$, and define $C_{0}=(1-\rho_{0}^{2})^{-1}$.

In the perfect matching model, $Z$ is distributed as the number of fixed points in a random permutation over $\{1,\ldots,k\}$. It is well known that

We prove that $R^{*}\to 1$ in two main regimes with the help of Theorem 2. To simplify notation, we assume that $k$ is even and recall that $n=k^{2}$ in this model.

We choose $a\to\infty$ such that $a^{2}\rho\sqrt{k\max(k,m)}\to 0$. We use the bounds $\cosh(x)\leq 1+x^{2}$ for $x\in(0,1)$ and $Z\leq k$, and the fact that $a^{2}\rho k\to 0$, to get, for $n$ sufficiently large,

Now let $c=C_{0}^{2}a^{4}\rho^{2}mk$. Using (27), one obtains

because $c\to 0$ and $\log[(k/2+1)!]\sim(k/2)\log k$ as $k\to\infty$.

We choose $a\to\infty$ such that $a^{2}\rho m/\log(\min(k,m))\to 0$. Using (27), one obtains

Now we take care separately of these last two terms. First, note that

For the other term, the situation is slightly more subtle. Let $Y$ be a sum of $m$ independent Rademacher random variables. Using the binomial identity, it is easy to prove that

Now thanks to Hoeffding’s inequality, we obtain for any $t>0$,

Consider the class of perfect matchings on the complete bipartite graph. Under either of (28), or (29), $R^{*}\to 1$.

It is easy to derive upper bounds for the performance of the localized squared-sum test in this model. All we need to observe is that $|\mathcal{C}|=k!$ and therefore $\log|\mathcal{C}|\sim k\log k$ when $k\to\infty$. Plugging this into (10), we see that the local squared-sum test is asymptotically powerful when

The clique number of random geometric graphs

In this section we describe a, perhaps unexpected, application of Theorem 2. We use this theorem to derive a lower bound for the clique number of random geometric graphs on high-dimensional spheres.

(i.e., the probability that an edge is present equals $p$). The clique number $\omega(n,d,p)$ is the size of the largest clique of $G(n,d,p)$ (i.e., the largest fully connected subset of vertices). In Devroye, György, Lugosi and Udina the behavior of the random variable $\omega(n,d,p)$ is studied for fixed values of $p$ when $n$ is large and $d=d_{n}$ grows as a function of $n$. The rate of growth of $\omega(n,d,p)$ is shown to depend in a crucial way of how fast $d_{n}$ increases with $n$. Specifically, the following results are established (and hold with probability converging to $1$ as $n\to\infty$):

The above-mentioned results leave open the question of where exactly the “phase transition” occurs, and whether the upper bound in the regime $d_{n}\sim\log^{2}n$ is sharp. In this section we are able to answer both of these questions. Below we establish a general lower bound for the clique number which implies that, perhaps surprisingly, the phase transition occurs around $\log^{2}n$ and that the upper bounds above cannot be improved in an essential way. We show that the median of the clique number $\omega(n,d,p)$ is bounded from below by $\exp(\kappa\log^{2}n/d)$ where $\kappa$ is a positive constant that depends on $p$ only. This implies, for example, that if $d\sim c\log n$ for some $c>0$, then $\omega(n,d,p)$ grows as a positive power of $n$. On the other hand, even when $d\sim\log^{2-\epsilon}n$ for any fixed $\epsilon>0$, then $\omega(n,d,p)$ is much larger than any power of $\log n$. For the sake of simplicity, we only state the result for the case of $p=1/2$. The argument is identical for other values of $p$.

There exist universal constants $c_{1},c_{2},c_{3},c_{4}>0$ such that for all $n,d$ such that $d\geq c_{1}\log(c_{2}n)$, the median of the clique number $\omega(n,d,1/2)$ satisfies

One may take $c_{1}=7/16$, $c_{2}=16\log 2$, $c_{3}=1/16$, and $c_{4}=49/5120$. In particular,

The basic idea of the proof is to define a test that works well whenever the median clique number is small. But then the lower bound of Theorem 2 implies that the clique number cannot be small.

Under the null hypothesis (when $\rho=0$), the $Z_{i}$’s are i.i.d. uniform on the sphere $S_{d-1}$ implying that $G\sim G(n,d,1/2)$ and, consequently, $\omega\sim\omega(n,d,1/2)$. Devroye, György, Lugosi and Udina show that, under the alternative hypothesis, with probability at least $7/8$, the graph contains a clique of size $k$ whenever

where we used Markov’s inequality in the last line.

Combining the bounds on the probabilities of type I and type II errors, we conclude that $R^{*}\leq 1/4$. Put it another way,

We conclude that, for any $\rho\in(0,1)$,

Discussion

We close this paper by discussing some open problems and directions of further research.

Sharper bounds. The cornerstone of our analysis is the lower bound stated in Theorem 2. It is powerful enough that we can deduce useful bounds in many different models, which are seen to be optimal up to constant or logarithmic factors. While a considerable effort has been devoted in the related detection-of-means problem for finding the right constants, one wonders if it is possible to obtain results that fine here, at least in some regimes. One possible avenue is via the truncated second moment approach, which underlies the lower bounds in Ingster , Donoho and Jin , Hall and Jin , Butucea and Ingster . The computations are rather daunting in the setup of this paper and we decided not to take this route. Note that the second moment approach (without truncation) has limited applicability, though it is a little more useful here than it is in the case where $m=1$.

Comparison with the detection-of-means setting. Our results reveal that the dependence on the sample size is different here. In the detection-of-means setting, one reduces by sufficiency to the case where $m=1$ by simply averaging the multiple observations and working with $\bar{X}_{1},\dots,\bar{X}_{n}$, where $\bar{X}_{i}=\frac{1}{m}\sum_{t}X_{t,i}$. Therefore, if initially $X_{t,i}\sim\mathcal{N}(\mu,1)$ when $i$ is anomalous, we now have $\bar{X}_{i}\sim\mathcal{N}(\mu,1/m)$. Therefore, we reduce the problem to where $m=1$ and $\mu$ is replaced by $\sqrt{m}\mu$. From this, we know that reliable detection is possible if either

where $C$ is a large enough constant. In our previous work, we argued that, at least when $\rho$ is bounded away from 1, the parameter $\rho$ in the correlation-detection problem played a similar role as $\mu^{2}$ in the detection-of-means setting. The case when the sample size $m\to\infty$ is, however, quite different both in the “dense” regime $\rho\sqrt{m}\gg n/k^{2}$, and in the “intermediary” regime $\rho\sqrt{m}\succ\sqrt{\log(n/k)/k}$. We also note that this regime does seem to have an equivalent in the detection-of-means setting.

General correlations. More generally, the problem of detecting correlations of arbitrary sign – not just positive correlations like we do here – remains open. Even though one can design natural tests akin to our squared-sum and local squared-sum tests for that situation, the challenge is in deriving tight lower bounds. We mention that our approach to obtaining a lower bound in Section 2 does not apply here, since the representation (2) is not valid when the correlations may be negative.

Acknowledgements

We would like to thank the anonymous referees for helpful comments and suggestions. We also thank Quentin Berthet and Philippe Rigollet for shedding some light on their results at the Nonparametric and High-dimensional Statistics conference, held in December of 2012, in Luminy, France. EAC was partially supported by NSF grant DMS-11-20888 and ONR grant N00014-09-1-0258. GL was supported by the Spanish Ministry of Science and Technology grant MTM2012-37195.

References