Distribution-free Detection of a Submatrix

arXiv:1604.07449v1 [] 25 Apr 2016

Ery Arias-Castro∗

Yuchao Liu†

Abstract. We consider the problem of detecting the presence of a submatrix with largerthan-usual values in a large data matrix. This problem was considered in (Butucea
and Ingster, 2013) under a one-parameter exponential family, and one of the test they
analyzed is the scan test. Taking a nonparametric stance, we show that a calibration
by permutation leads to the same (first-order) asymptotic performance. This is true for
the two types of permutations we consider. We also study the corresponding rank-based
variants and precisely quantify the loss in asymptotic power.

1

Introduction

Biclustering has emerged as an important set of tools in bioinformatics, in particular, in the analysis
of gene expression data (Cheng and Church, 2000). It comes in different forms, and in fact,
the various methods proposed under that umbrella may target different goals. See (Madeira and
Oliveira, 2004) for a survey. Here we follow (Shabalin et al., 2009), where the problem is posed as
that of discovering a submatrix of unusually large values in a (large) data matrix. For example,
in the context of a microarray dataset, the data matrix is organized by genes (rows) and samples
(columns). We let X = (Xij ) denote the matrix, M denote the number of rows and N denote the
number of columns, so the data matrix X is M -by-N .

1.1

Submatrix detection

In its simplest form, there is only one submatrix to be discovered. In that context, the detection
problem is that of merely detecting of the presence of an anomalous (or unusual) submatrix, which
leads to a hypothesis testing problem. This was considered in (Butucea and Ingster, 2013) from a
minimax perspective. Their work relies on parametric assumptions. For example, in the normal
model, they assume that the Xij ’s are independent and normal, with mean θij and unit variance.
Under the null hypothesis, θij = 0 for all i ∈ [M ] ∶= {1, . . . , M } and all j ∈ [N ]. Under the alternative,
there is a m-by-n submatrix indexed by Itrue ⊂ [M ] and Jtrue ⊂ [N ] such that
θij ≥ θ‡ ,

∀(i, j) ∈ Itrue × Jtrue ,

(1)

while θij = 0 otherwise. Here θ‡ > 0 controls the signal-to-noise ratio. In that paper, Butucea and
Ingster precisely establish how large θ‡ needs to be as a function of (M, N, m, n) in order for there
to exist a procedure that has (worst-case) risk tending to zero in the large-sample limit (i.e., as
the size of the matrix grows). They consider two tests which together are shown to be minimax
optimal. One is the sum test based on
sum(X) = ∑

∑ Xij .

i∈[M ] j∈[N ]
∗
†

University of California, San Diego — http://www.math.ucsd.edu/~eariasca/
University of California, San Diego — http://www.math.ucsd.edu/~yul085/

1

(2)

2
It is most useful when the submatrix is large. The other one is the scan test which, when the
submatrix size is known (meaning m and n are known) is based on
scan(X) =

max

I⊂[M ],∣I∣=m

max

J ⊂[N ],∣J ∣=n

∑ ∑ Xij .

(3)

i∈I j∈J

When m and n are unknown, one can perform a scan test for each (m, n) in some range of interest
and control for multiple testing using the Bonferroni method. From (Butucea and Ingster, 2013),
and also from our own prior work, we know that the resulting procedure achieves the same firstorder asymptotic performance.
To avoid making parametric assumptions, some works such as (Barry et al., 2005; Hastie et al.,
2000) have suggested a calibration by permutation. We consider two somewhat stylized permutation
approaches:
• Unidimensional permutation. The entries are permuted within their row. (One could permute
within columns, which is the same after transposition.)
• Bidimensional permutation. The matrix is vectorized, the entries are permuted uniformly at
random as one would in a vector, and the matrix is reformed.
The first method is most relevant when one is not willing to assume that the entries in different
rows are comparable. It is appealing in the context of microarray data and was suggested, for
example, in (Hastie et al., 2000). The second method is most relevant in a setting where all the
variables are comparable. In the parlance of hypothesis testing, the first method derives from a
model where the entries within each row are exchangeable under the null, while the second method
arises when assuming that all the entries are exchangeable under the null.
Contribution 1 (Calibration by permutation). We analyze the performance of the scan test when
calibrated using one of these two permutation approaches. We show that, regardless of the variant,
the resulting test is (first order) asymptotically as powerful as a calibration by Monte Carlo with
full knowledge of the parametric model. We prove this under some standard parametric models.
Remark 1. We focus on the scan statistic (3) and abandon the sum statistic (2) for at least two
reasons: 1) the sum statistic cannot be calibrated without knowledge of the null distribution; 2)
the sum statistic is able to surpass the scan statistic when it is impossible to locate the submatrix
with any reasonable accuracy, which is somewhat less interesting to the practitioner.
A calibration by permutation is computationally intensive in that it requires the repeated computation of the test statistic on permuted data. In practice, several hundred permutations are used,
which can cause the method to be rather time-consuming. A possible way to avoid this is to use
ranks, which was traditionally important before the availability of computers with enough computational power. (Hettmansperger, 1984) is a classical reference. In line with the two permutation
methods described above, we consider the corresponding methods for ranking the entries:
• Unidimensional ranks. The entries are ranked relative to the other entries in their row.
• Bidimensional ranks. The entries are ranked relative to the all other entries.
The use of ranks has the benefit of only requiring calibration (typically done on a computer nowadays) once for each matrix size M × N . It has the added benefit of yielding a method that is much
more robust to outliers.

3
Contribution 2 (Rank-based method). We analyze the performance of the scan test when the
entries are replaced by their ranks following one of the two methods just described. We show that,
regardless of the variant, there is a mild loss of asymptotic power, which we precisely quantify. We
do this under some standard parametric models.

1.2

More related work

The scan statistic (3) is computationally intractable and there has been efforts to offer alternative
approaches. We already mentioned (Shabalin et al., 2009), which proposes an alternate optimization
strategy: given a set of rows, optimize over the set of columns, and vice versa, alternating in this
fashion until convergence to a local maximum. This is the algorithm we use in our simulations.
It does not come with theoretical guarantees (other than converging to a local maximum) but
performs well numerically. A spectral method is proposed in (Cai et al., 2015) and a semidefinite
relaxation is proposed in (Chen and Xu, 2014). These methods can be run in time polynomial in
the problem size (meaning in M and N ). (Ma and Wu, 2015) establishes a lower bound based
on the Planted Clique Problem that strictly separates the performance of methods that run in
polynomial time from the performance of the scan statistic.
Our work here is not on the computational complexity of the problem. Rather we assume that
we can compute the scan statistic and proceed to study it. In effect, we contribute here to a long
line of work that studies permutation and rank-based methods for nonparametric inference. Most
notably, we continue our recent work (Arias-Castro et al., 2015) where we study the detection
problem under a similar premise but under much more stringent structural assumptions. The
setting there would correspond to an instance where the submatrix is in fact a block, meaning,
that Itrue and Jtrue are of the form Itrue = {i + 1, . . . , i + k} and Jtrue = {j + 1, . . . , j + l}. The
present setting assumes much less structure. The related applications are very different in the end.
Nevertheless, the technical arguments developed there apply here with only minor adaptation. The
main differences are that we consider two types of permutation and ranking protocols.

1.3

Content

The rest of the paper is organized as follows. In Section 2 we describe a parametric setting where
likelihood methods have been shown to perform well. This parametric setting will serve as benchmark for the nonparametric methods that ensue. In Section 3 we consider the detection problem
and study the scan statistic with each of the two types of calibration by permutation. In Section 4
we consider the same problem and study the rank-based scan statistic using each of the two types
of rankings. In Section 5 we present some numerical experiments on simulated data. All the proofs
are in Section 6.

2

The parametric scan

Following the classical line in the literature on nonparametric tests, we will evaluate the nonparametric methods introduced later on a family on parametric models. As in (Butucea and Ingster,
2013), and in our preceding work (Arias-Castro et al., 2015), we consider a one-parameter exponential family in natural form.
To define such a family, fix a probability distribution ν on the real line with zero mean and unit
variance, and with a sub-exponential right tail, specifically meaning that ϕ(θ) ∶= ∫R eθx ν(dx) < ∞
for some θ > 0. Let θ⋆ denote the supremum of all such θ > 0. (Note that θ⋆ may be infinite.) The

4
family is then parameterized by θ ∈ [0, θ⋆ ) and has density with respect to ν defined as
fθ (x) = exp{θx − log ϕ(θ)}.

(4)

By varying ν, we obtain the normal (location) family, the Poisson family (translated to have zero
mean), and the Rademacher family.
Such a parametric model is attractive as a benchmark because it includes these popular models
and also because likelihood methods are known to be asymptotically optimal under such a model.
Butucea and Ingster (2013) showed this to be the case for the problem of detection, where the
generalized likelihood ratio test is based on the scan statistic (3).
Under such a parametric model, the detection problem is formalized as a hypothesis testing
problem where ν plays the role of null distribution. In detail, suppose that the submatrix is known
to be m × n. The search space is therefore
Sm,n ∶= {S = I × J ∶ I ⊂ [M ], ∣I∣ = m and J ⊂ [N ], ∣J ∣ = n}.

(5)

We assume that the Xij ’s are independent with Xij ∼ fθij , and the testing problem is
H0 ∶ θij = 0,

∀(i, j) ∈ [M ] × [N ],

(6)

versus

⎧
⎪
⎪θij ≥ θ‡ , ∀(i, j) ∈ Strue ,
H1 ∶ ∃Strue ∈ Sm,n such that ⎨
⎪
⎪
⎩θij = 0, otherwise.
Here θ‡ controls the signal-to-noise ratio is assumed to be known in this formulation.
In this context, we have the following.

(7)

Theorem 1 (Butucea and Ingster (2013)). Consider an exponential model as described above, with
ν having finite fourth moment. Assume that
m n
log(M ∨ N )
,
→ 0,
→ 0.
M N
m∧n
Then the sum test based on (2), at any fixed level α > 0, has limiting power 1 when
mn
θ‡ √
→ ∞.
MN
M, N, m, n → ∞,

Then the scan test based on (3), at any fixed level α > 0, has limiting power 1 when
√
θ‡ mn
lim inf √
> 1.
N
+
n
log
)
2(m log M
m
n

(8)

(9)

(10)

Conversely, the following matching lower bound holds. Assume in addition that log M ≍ log N and
m ≍ n. Then any test at any fixed level α > 0 has limiting power at most α when
√
θ‡ mn
mn
θ‡ √
→ 0 and lim inf √
< 1.
(11)
N
MN
2(m log M
+
n
log
)
m
n
We note that Butucea and Ingster (2013) derived their lower bound under slightly weaker
assumptions on M, N, m, n.
Remark 2. Proper calibration in this context is based on knowledge of the null distribution ν. In
more detail, consider a test that rejects for large values of a statistic T (X). Assuming a desired
level of α > 0 and that ν is either diffuse or discrete (for simplicity), the critical value for T is set
at tα , where tα = inf{t ∶ ν(T (X) ≥ t) ≤ α}. The test is then I{T (X) ≥ tα }. In practice, tα may be
approximated by Monte Carlo sampling.

5

3

Permutation scan tests

In the previous section we described the work of Butucea and Ingster (2013), who in certain
parametric models show that the sum test (2) and scan test (3) are jointly optimal for the problem
of detecting a submatrix. This is so if they are both calibrated with full knowledge of the null
distribution (denoted ν earlier).
What if the null distribution is unknown? A proven approach is via permutation. This is
shown to be optimal in some classical settings (Lehmann and Romano, 2005) and was recently
shown to also be optimal in more structured detection settings (Arias-Castro et al., 2015). We
prove that this is also the case in the present setting of detecting a submatrix. We consider the two
types of permutation, unidimensional and bidimensional, described in Section 1.1. More elaborate
permutation schemes have been suggested, e.g., in (Barry et al., 2005), but these are not considered
here, in part to keep the exposition simple. Indeed, we simply aim at showing that a calibration
by permutation performs very well in the present context.
Let Π be a subgroup of permutations of [M ] × [N ], identified with [M N ]. Then a calibration
by permutation of the scan statistic (or any other statistic) yields the P-value
P(X) =

#{π ∈ Π ∶ scan(Xπ ) ≥ scan(X)}
,
∣Π∣

(12)

where Xπ = (Xπ(i,j) ) is the matrix permuted by π. The permutation scan test at level α is the test
I{P(X) ≤ α}. It is well-known that this this a valid P-value, in the sense that, under the null, it
dominates the uniform distribution on [0, 1] (Lehmann and Romano, 2005). (This remains true of
a Monte Carlo approximation.)
The set of unidimensional permutations, denoted Π1 , is that of all permutations that permute
within each row, while the set of bidimensional permutations, denoted Π2 , is simply the set of all
permutations. Obviously, Π1 ⊂ Π2 with ∣Π1 ∣ = (N !)M and ∣Π2 ∣ = (M N )!, and they are both groups.1
Theorem 2. Consider an exponential model as described in Section 2. In addition to (8), assume
log3 (M ∨ N )/(m ∧ n) → 0,

(13)

and that either (i) ν has support bounded from above, or (ii) maxi,j θij ≤ θ̄ for some θ̄ < θ⋆ fixed.
Let the group of permutations Π be either Π1 or Π2 ; if Π = Π1 , we require that ϕ(θ) < ∞ for some
θ < 0. Then the permutation scan test based on (12), at any fixed level α > 0, has limiting power 1
when (10) holds.
The additional condition (on ν or the nonzero θij ’s) seems artificial, but just as in (Arias-Castro
et al., 2015), we are not able to eliminate it. Other than that, in view of Theorem 1 we see that
the permutation scan test — just like the parametric scan test — is optimal to first-order under a
general one-parameter exponential model.

4

Rank-based scan tests

Rank tests are classical special cases of permutation tests (Hettmansperger, 1984). Traditionally,
when computers were not as readily available and not as powerful, permutation tests were not
practical, but rank tests could still be, as long as calibration had been done once for the same (or
1

The group structure is important. See the detailed discussion in (Hemerik and Goeman, 2014).

6
a comparable) problem size. Another well-known advantage of rank tests is their robustness to
outliers.
We consider the two ranking protocols described in Section 1.1. After the observations are
ranked, the distribution under the null is the permutation distribution, either uni- or bi-dimensional
depending on the ranking protocol. This is strictly true under an appropriate exchangeability
condition, which holds in the null model we consider here where all observations are IID. In fact, the
unidimensional rank scan test is a form of unidimensional permutation test, and the bidimensional
rank scan is a form of bidimensional permutation test, each time, the statistic being the rank scan
scan(R) =

max

I⊂[M ],∣I∣=m

max

J ⊂[N ],∣J ∣=n

∑ ∑ Rij ,

(14)

i∈I j∈J

where R = (Rij ) is the matrix of ranks.
Rank tests have been studied in minute detail in the classical setting (Hájek and Sidak, 1967;
Hettmansperger, 1984). Typically, this is done, again, by comparing their performance with the
likelihood ratio test in the context of some parametric model. Typically, there is some loss in
efficiency, unless one tailors the procedure to a particular parametric family.2 Such a performance
analysis was recently carried out for the rank scan in more structured settings (Arias-Castro et al.,
2015). We again extend this work here and obtain the following.
Define
1
Υ = E(Z 1(Z>Y ) ) + E(Z 1(Z=Y ) )
2
where Y, Z are IID with distribution ν. (This is the same constant introduced by Arias-Castro
et al. (2015).)
Theorem 3. Consider an exponential model as described in Section 2. Assume that (8) holds. Let
the group of permutations Π be either Π1 or Π2 . The rank scan test at any fixed level α > 0 has
limiting power 1 when
√
θ‡ mn
1
lim inf √
(15)
> √ .
N
2 3Υ
2(m log M
+
n
log
)
m
n
The proof is omitted as it is entirely based on an adaptation of that of Theorem 2 and arguments
given in (Arias-Castro et al., 2015) to handle the rank moments.
Compared to the (optimal) performance of the parametric and permutation scan tests in the
same setting (Theorem 1 and Theorem 2), we see that there is a loss in power. However, the loss
can √
be quite√small. For example, as argued in (Arias-Castro et al., 2015), in the normal model
1/(2 3Υ) ≤ π/3 ≈ 1.023.

5

Numerical experiments

We performed some numerical experiments3 to assess the accuracy of our asymptotic theory. To do
so, we had to deal with two major issues in terms of computational complexity. The first issue is
the computation of the scan statistic defined in (3). There are no known computationally tractable
method for doing so. As Butucea and Ingster (2013) did, we opted instead for an approximation in
the form of the alternate optimization (or hill-climbing) algorithm of Shabalin et al. (2009). Since
in principle this algorithm only converges to a local maximum, we run the algorithm on several
2
3

Actually, Hajek (1962) proposes a more complex method that avoids the need for knowing the null distribution.
In the spirit of reproducible research, our code is publicly available at https://github.com/nozoeli/NPDetect

7
random initializations and take the largest output. The second issue is that of computing the
permutation P-value defined in (12). (This is true for the permutation test and also for the special
case of the rank test.) Indeed, examining all possible permutations in Π (either Π1 or Π2 ) is only
feasible for very small matrices. As usual, we opted for Monte Carlo sampling. Specifically, we
picked π1 , . . . πB IID uniform from Π with B = 500 in our setup. We then estimate the permutation
P-value by
#{b ∈ [B] ∶ scan(Xπb ) ≥ scan(X)}
P̂(X) =
.
(16)
B+1
We mention that when rank methods are applied, the ties in the data are broken randomly.
Simulation setup Our simulation strategy is as follows. A data matrix X of size M × N is
generated with the anomaly as [m] × [n]. All the entries of X are independent with distribution
f0 (same as ν) except for the anomalous ones which have distribution fθ‡ for some θ‡ > 0. We
compare the permutation tests and rank tests (unidimensional and bidimensional) with the scan
test calibrated by Monte Carlo (using 500 samples), which serves as an oracle benchmark as it has
full knowledge of the null distribution f0 . By construction, all tests have the prescribed level. As
we increase θ‡ , the P-values of the different tests are recorded. Each setting is repeated 200 times.
As one of the main purposes of our simulations is to confirm our theory, we zoom in on the
region near the critical value
¿
Á 2(m log M + n log N )
Á
À
m
n
θcrit =
,
(17)
mn
which comes from (10). Specifically, we increase θ‡ from 0.5 × θcrit to 1.5 × θcrit with step size
0.125 × θcrit to explore the behavior of P-values around the critical value.
The Normal Case Here we generate data from normal family, where fθ is N (θ, 1). We used
two setups, (M, N, m, n) = (200, 100, 10, 15) and (M, N, m, n) = (200, 100, 30, 10), to assess the
performance of the tests under different anomaly sizes. The resulting boxplots of the averaged
P-values are shown in Figure 1.
From the plots we see that the P-values are generally very close to 0 when θ‡ exceeds θcrit .
When (m, n) = (10, 15) the convergence towards 0 is slower, which may be due to the small size
of the anomalous submatrix. As expected, the (oracle) Monte Carlo test is best, followed by the
bidimensional permutation test, followed by the unidimensional permutation test. That said, the
differences appear to be minor, which confirms our theoretical findings.
For the rank tests, we observe a similar behavior of the P-values, with the bidimensional showing
superiority over the unidimensional rank test, but the loss of power with respect
√ to the oracle test
is a bit more substantial, as predicted by the theory. As shown before, 1/(2 3Υ) ≈ 1.03 for the
standard normal, so that we should place the critical threshold approximately at 1.03 × θcrit . This
appears to be confirmed in the setting where (m, n) = (30, 10). While the P-values for the rank
tests converge relatively slowly when (m, n) = (10, 15) (for unidimensional rank test the P-value is
close to 0 at θ = 1.5 × θcrit ), this may be due to the relatively small size of the anomaly.
The Poisson Case As another example, we consider the Poisson family, where fθ corresponds
to Poisson(eθ ) − 1. The data matrix and anomaly sizes are the same as they are in the normal
case. The resulting boxplots of the P-values are shown in Figure 2. Overall, we observe a similar
behavior of the P-values.

8

Figure 1: P-values of various forms of scan tests in the normal model

9

Figure 2: P-values of various forms of scan tests in the Poisson model

10

6
6.1

Proofs
Preliminaries

We start with some preliminary results that already appear, one or another, in our previous work
(Arias-Castro et al., 2015). First, for any one-parameter exponential family (fθ ∶ θ ∈ Θ) with a
standardized base distribution ν, as we consider to be here,
Eθ (X) ≥ θ,

∀θ ∈ Θ.

(18)

Next, in the same context, if sup Θ > 0 (which we assume throughout), then fθ has a sub-exponential
right tail, which is uniform in θ ∈ θ̄ if θ̄ ∈ Θ. In particular, there is γ̄ that depends on θ̄ > 0 such
that, if X1 , . . . , Xk are independent, with Xj ∼ fθj with θj ≤ θ̄, then
max Xj ≤ γ̄ log k,

with probability tending 1 as k → ∞.

(19)

j∈[k]

By symmetry, if inf Θ < 0 (which we assume in the case of unidimensional permutations), the same
is true on the left. In particular, ν itself (corresponding to θ = 0) has a sub-exponential left tail in
this case, and this is all that will be used below. In particular, there is a constant γ0 > 0 such that,
if X1 , . . . , Xk are IID ν, then
min Xj ≥ −γ0 log k,

with probability tending 1 as k → ∞.

(20)

j∈[k]

6.2

Proof of Theorem 2

The proof arguments are parallel to those of Arias-Castro et al. (2015), derived in the context of
more structured settings. For that reason, we only detail the proof of Theorem 2 in the case of
unidimensional permutations, which, compared to bidimensional permutations, is a bit more different from the setting considered in (Arias-Castro et al., 2015) and requires additional arguments.
Therefore, in what follows, we take Π = Π1 . Recall that in this case we assume in addition that
ϕ(θ) < ∞ for some θ < 0. This implies that ν as sub-exponential tails
Case (i) We first focus on the condition where ν has support bounded from above and let b0
denote such an upper bound. (Necessarily, b0 > 0.) Thus, regardless of the θij ’s,
P(maxi,j Xij ≤ b0 ) = 1.

(21)

The permutation scan test has limiting power 1 if and only if P(P(X) ≤ α) → 1 under the alternative. We show that by proving the stronger claim that P(X) → 0 in probability under the
alternative.
We first work conditional on X = x, where x = (xij ) denotes a realization of X = (Xij ). We
may equivalently center the rows of X before scanning, and the resulting test remains unchanged.
Therefore, we may assume that all the rows of x sum to 0. Let ζ = scan(x) for short. We have
P(x) = P(scan(xπ ) ≥ ζ),

(22)

where the randomness comes solely from π, uniformly drawn from Π. Using the union bound, we
get
P(x) ≤ ∣Sm,n ∣ max P ( ∑(i,j)∈S xπ(i,j) ≥ ζ).
(23)
S∈Sm,n

11
For each i ∈ [N ], let (Aij ∶ j ∈ [n]) be a sample from (xij ∶ j ∈ [N ]) without replacement and let
Ai = ∑j∈[n] Aij . Note that A1 , . . . , AM are independent and, for S = I × J , we have
∑ xπ(i,j) ∼ ∑ Ai .
(i,j)∈S

(24)

i∈I

Fix I ⊂ [M ] of size m. Using Markov’s inequality and the independence of the Ai ’s, we get
P ( ∑i∈I Ai ≥ ζ) ≤ e−cζ ∏i∈I φi (c),

(25)

where φi is the moment generating function of Ai . The key is (Hoeffding, 1963, Th 4), which
implies that φi ≤ ψi , where ψi is the moment generating function of Bi , where Bi = ∑j∈[n] Bij and
(Bij ∶ j ∈ [n]) is a sample from (xij ∶ j ∈ [N ]) with replacement, meaning that these are IID random
variables uniformly distributed in (xij ∶ j ∈ [N ]). By (21), we have Bij ≤ b0 with probability one,
and the usual arguments leading to the (one-sided) Bernstein’s inequality yield the usual bound
ψi (c) ≤ exp [

nc2 σi2 ecb0 − 1 − cb0
],
2
c2 b20 /2

where σi2 is the variance of Bi1 , meaning, σi2 =
mean. Letting σ 2 = maxi∈[M ] σi2 , we derive
e−cζ ∏ φi (c) ≤ e−cζ ∏ exp [
i∈I

i∈I

1
N

(26)

∑j∈[N ] (xij − x̄i )2 , with x̄i =

1
N

∑j∈[N ] xij being the

nc2 σi2 ecb0 − 1 − cb0
mnc2 σ 2 ecb0 − 1 − cb0
−cζ
]
≤
e
exp
[
],
2
2
c2 b20 /2
c2 b20 /2

(27)

the latter being the usual bound that leads to Bernstein’s inequality. The same optimization over
c yields
P (∑i∈I Ai ≥ ζ) ≤ exp [−

ζ2
].
2mnσ 2 + 32 b0 ζ

(28)

We now emphasize the dependency of ζ and σ 2 on x by adding x as a subscript. Noting that this
bound is independent of I (of size m), we get
P(x) ≤ ∣Sm,n ∣ exp [−

ζx2
].
2mnσx2 + 23 b0 ζx

(29)

2
We now free X and bound ζX from below, and σX
from above. When doing so, we need to take
into account that we assumed the rows summed to 0. When this is no longer the case, ζX denotes
the scan of X after centering all the rows. Let X̄i denote the mean of row i. By definition of the
scan in (3),

ζX ≥ ζtrue ∶= ∑

∑ (Xij − X̄i ) = (1 −

i∈Itrue j∈Jtrue

n
N)

∑

∑ Xij −

i∈Itrue j∈Jtrue

n
N

∑

∑ Xij .

(30)

i∈Itrue j∉Jtrue

For the expectation, by (8) and (18), we have
E(ζtrue ) ≥ (1 −

n
N)

∑

∑ θij ≥ (1 − o(1))mnθ‡ .

(31)

i∈Itrue j∈Jtrue

For the variance, we have Var(Xij ) = 1 when (i, j) ∉ Strue (since ν has variance 1) and Var(Xij ) ≤
E(Xij2 ) ≤ b20 always. Using this, we derive
n 2
Var(ζtrue ) ≤ mnb20 + ( N
) mN = mn(b20 +

n
N)

= O(mn).

(32)

12
√
Because of (8) and (10), E(ζtrue ) ≫

Var(ζtrue ), and thus by Chebyshev’s inequality,

ζtrue = (1 + oP (1)) E(ζtrue ) ≥ (1 + oP (1))mnθ‡ .

(33)

We now bound σx2 . For i ∈ Itrue , we have
σi2 (X) ≤

1
1
2
∑ X =
N j∈[N ] ij N

2
∑ Xij +
j∈Jtrue

For i ∉ Itrue ,
σi2 (X) ≤
Therefore

1
N

2
∑ Xij ≤
j∉Jtrue

nb20 1
+
N
N

2
∑ Xij .

1
2
∑ X .
N j∈[N ] ij
1
∑ Tij ,
i∈[M ] N j∈[N ]

sto

(34)

j∉Jtrue

2
σX
≤ 1 + o(1) + max

(35)

(36)

where (Tij ∶ (i, j) ∈ [M ] × [N ]) are IID with distribution that of X 2 − 1 when X ∼ ν. Note that
E(Tij ) = 0 since ν has variance 1 and
max Tij ≤ t̄ ∶= b20 ∨ (γ0 log(M N ))2 ,
i,j

(37)

by (21) and when the following event holds
A ∶= { min Xij ≥ −γ0 log(M N )},
i,j

(38)

which by (20) happens with probability tending to 1. Let PA be the probability conditional on A
and EA the corresponding expectation. Let µA = EA (Tij ) and τA2 = VarA (Tij ) < ∞, because ν has
finite fourth moment. By Bernstein’s inequality, for any c > µA ,
N (c − µA )2
1
].
∑ Tij > c) ≤ exp [ −
N j∈[N ]
2τA2 + 23 t̄c

(39)

1
N (c − µA )2
].
∑ Tij > c) ≤ M exp [ −
i∈[M ] N j∈[N ]
2τA2 + 23 t̄c

(40)

PA (
Then using a union bound

PA ( max

Taking logs, noting that µA → 0 and τA2 → τ 2 ∶= Var(Tij ), as well as t̄ = O(log(M N )), and using (8)
and (13), we see that the RHS tends to 0 for any c > 0 fixed. Therefore maxi∈[M ] N1 ∑j∈[N ] Tij = oP (1)
conditional on A, and since P(A) → 1, also unconditionally. Coming back to (36), we conclude that
2
σX
= 1 + oP (1).

(41)

2
The upper bound on ζX and the lower bound on σX
, combined, imply by monotonicity that
2
mnθ‡2
ζX
≥ (1 + oP (1))
.
2 + 2b ζ
2mnσX
2 + 23 b0 θ‡
3 0 X

(42)

)(N
), so that
We also have ∣Sm,n ∣ = (M
m
n
)
(N )
log ∣Sm,n ∣ = log (M
m + log n ≤ (1 + o(1))Λ,

(43)

13
with
N
],
Λ ∶= [m log M
m + n log n

(44)

)
where in the last inequality we used (8) and the fact that log (K
k ≤ k log(K/k) + k for all integers
1 ≤ k ≤ K.
Coming back to (29) and collecting all the bounds in between, we find that
log P(X) ≤ (1 + o(1))Λ − (1 + oP (1))

mnθ‡2
2 + 23 b0 θ‡

.

(45)

Under (10), there is ε > 0 such that, eventually,
√
θ‡ ≥ (1 + ε) 2Λ/(mn).

(46)

When that’s the case, we get
log P(X) ≤ (1 + o(1))Λ − (1 + oP (1))

1+

(1 + 2ε)Λ
√
.
2Λ/(mn)

(47)

1
3 b0 (1 + ε)

Noting that Λ/(mn) = o(1) and Λ → ∞ under (8), we get
log P(X) ≤ −(1 + oP (1))2εΛ → −∞,

(48)

which is what we needed to prove.
Case (ii) We now consider the case where θij ≤ θ̄ for all (i, j) ∈ [M ] × [N ] for some θ̄ < θ∗ .
Although (21) may not hold for any b0 , we redefine b0 = γ̄ log(M N ), where γ̄ depends on θ̄, and
condition on the event
B ∶= { max Xij ≤ b0 },
(49)
i,j

which holds with probability tending to 1 by (19). The bound (29) holds unchanged (assuming
2
are handled, now that we conditioned on
that maxi,j xij ≤ b0 ). What is different is how ζX and σX
B. Let PB and EB denote the probability and expectation conditional on B.
We have
EB (ζtrue ) ≥ (1 −
≥ (1 −

n
N)
n
N)

∑

∑ EB (Xij ) −

i∈Itrue j∈Jtrue

∑

n
N

∑ E(Xij ∣Xij ≤

i∈Itrue j∈Jtrue

≥ (1 + o(1))mnθ‡ .

∑ ∑ EB (Xij )
i∈Itrue j∉Jtrue
n
b0 ) − N
∑ ∑ E(Xij ∣Xij
i∈Itrue j∉Jtrue

(50)
≤ b0 )

(51)
(52)

In the last inequality, for j ∉ Jtrue we used the fact that E(Xij ) = 0, which implies that EB (Xij ) ≤ 0
in that case. And for j ∈ Jtrue we used the fact that E(Xij ∣Xij ≤ b0 ) → θij ≥ θ‡ combined with a
Cèsaro-type argument. On the other hand, in a similar way, we also have
VarB (ζtrue ) = O(mnb20 ) = O(mn log2 (M N )).
(53)
√
So we still have EB (ζtrue ) ≫ VarB (ζtrue ), by (8) and (10), and in addition (13). In particular,
(33) holds under B. In very much the same way, one can verify that the same is true of (41).
From there we get to (47) in exactly the same way, conditional on
√ B, and then unconditionally
since P(B) → 1. Then, to conclude, we only need to check that b0 Λ/(mn) = o(1), which is the
case by (13).

14

Acknowledgements
This work was partially supported by a grant from the US Office of Naval Research (N00014-13-10257) and a grant from the US National Science Foundation (DMS 1223137).

References
Arias-Castro, E., R. M. Castro, E. Tánczos, and M. Wang (2015). Distribution-free detection of
structured anomalies: Permutation and rank-based scans. arXiv preprint arXiv:1508.03002 .
Barry, W. T., A. B. Nobel, and F. A. Wright (2005). Significance analysis of functional categories in
gene expression studies: a structured permutation approach. Bioinformatics 21 (9), 1943–1949.
Butucea, C. and Y. I. Ingster (2013). Detection of a sparse submatrix of a high-dimensional noisy
matrix. Bernoulli 19 (5B), 2652–2688.
Cai, T. T., T. Liang, and A. Rakhlin (2015). Computational and statistical boundaries for submatrix localization in a large noisy matrix. arXiv preprint arXiv:1502.01988 .
Chen, Y. and J. Xu (2014). Statistical-computational tradeoffs in planted problems and submatrix
localization with a growing number of clusters and submatrices. arXiv preprint arXiv:1402.1267 .
Cheng, Y. and G. M. Church (2000). Biclustering of expression data. In Proceedings of the Eighth
International Conference on Intelligent Systems for Molecular Biology, pp. 93–103. AAAI Press.
Hajek, J. (1962). Asymptotically most powerful rank-order tests. The Annals of Mathematical
Statistics 33 (3), 1124–1147.
Hájek, J. and Z. Sidak (1967). Theory of rank tests. Academic Press, Academia Publishing House
of the Czechoslovak Acad.
Hastie, T., R. Tibshirani, M. B. Eisen, A. Alizadeh, R. Levy, L. Staudt, W. C. Chan, D. Botstein,
and P. Brown (2000). ‘gene shaving’ as a method for identifying distinct sets of genes with similar
expression patterns. Genome Biology 1 (2), 1–21.
Hemerik, J. and J. Goeman (2014). Exact testing with random permutations. arXiv preprint
arXiv:1411.7565 .
Hettmansperger, T. P. (1984). Statistical inference based on ranks. Wiley Series in Probability
and Mathematical Statistics: Probability and Mathematical Statistics. New York: John Wiley
& Sons, Inc.
Hoeffding, W. (1963). Probability inequalities for sums of bounded random variables. J. Amer.
Statist. Assoc. 58, 13–30.
Lehmann, E. L. and J. P. Romano (2005). Testing statistical hypotheses (Third ed.). Springer
Texts in Statistics. New York: Springer.
Ma, Z. and Y. Wu (2015). Computational barriers in minimax submatrix detection. The Annals
of Statistics 43 (3), 1089–1116.
Madeira, S. C. and A. L. Oliveira (2004). Biclustering algorithms for biological data analysis: a
survey. IEEE/ACM Transactions on Computational Biology and Bioinformatics (TCBB) 1 (1),
24–45.
Shabalin, A. A., V. J. Weigman, C. M. Perou, and A. B. Nobel (2009). Finding large average
submatrices in high dimensional data. The Annals of Applied Statistics 3 (3), 985–1012.

