Optimal Bipartite Network Clustering

arXiv:1803.06031v1 [math.ST] 15 Mar 2018

Zhixin Zhou and Arash A. Amini
March 19, 2018

Abstract
We consider the problem of bipartite community detection in networks, or more generally the network biclustering problem. We present a fast two-stage procedure based on
spectral initialization followed by the application of a pseudo-likelihood classifier twice.
Under mild regularity conditions, we establish the weak consistency of the procedure (i.e.,
the convergence of the misclassification rate to zero) under a general bipartite stochastic
block model. We show that the procedure is optimal in the sense that it achieves the optimal convergence rate that is achievable by a biclustering oracle, adaptively over the whole
class, up to constants. The optimal rate we obtain sharpens some of the existing results and
generalizes others to a wide regime of average degree growth. As a special case, we recover
the known exact recovery threshold in the log n regime of sparsity. To obtain the general
consistency result, as part of the provable version of the algorithm, we introduce a subblock partitioning scheme that is also computationally attractive, allowing for distributed
implementation of the algorithm without sacrificing optimality. The provable version of the
algorithm is derived from a general blueprint for pseudo-likelihood biclustering algorithms
that employ simple EM type updates. We show the effectiveness of this general class by
numerical simulations.
Keywords: Bipartite networks; stochastic block model; community detection; biclustering; network analysis; pseudo-likelihood, spectral clustering.

1

Introduction

Network analysis has become an active area of research over the past few years, with applications and contributions from many disciplines including statistics, computer science, physics,
biology and social sciences. A fundamental problem in network analysis is detecting and
identifying communities, also known as clusters, to help better understand the underlying
structure of the network. The problem has seen rapid advances in recent years with numerous breakthroughs in modeling, theoretical understanding, and practical applications [FH16].
In particular, there has been much excitement and progress in understanding and analyzing
the stochastic block model (SBM) and its variants. We refer to [Abb17] for a recent survey
of the field. Much of this effort, especially on the theoretical side has been focused on the
univariate (or symmetric) case, while the bipartite counterpart, despite numerous practical
applications, has received comparatively less attention. Of course, there has been lots of activity in terms of modeling and algorithm development for bipartite clustering both in the
context of networks [ZRMZ07; LCJ14; WFL14; Roh15; RAL17] as well as other domains, such
as topic modeling and text mining [Dhi01; Dhi03], as well as in biological applications [CC00;
1

MTSCO10]. But much of this work either lacks theoretical investigations or has not considered
the issue of statistical optimality.
In this paper, we provide a unified treatment of the community detection, or clustering, in
the bipartite setting with a focus on deriving fundamental theoretical limits of the problem.
The main goal is to propose computationally feasible algorithms for bipartite network clustering
that exhibit provable statistical optimality. We will focus on the bipartite version of the SBM
which is a natural model for bipartite networks with clusters. SBM is a stochastic network
model where the probability of edge formation depends on the latent (unobserved) community
assignment of the nodes, often referred to as node labels. The goal of the community detection
problem is to recover these labels given an instance of the network. This is a non-trivial task
since, for example, maximum likelihood estimation involves a search over exponentially many
labels.
Community detection in bipartite SBM is closely related to the biclustering problem, for
which many algorithms have been developed over the years [Har72; CC00; TSS02; GLMZ16].
On the other hand, in recent years, many algorithms have been proposed for clustering in univariate SBMs, including global approaches such as spectral clustering [RCY11; Krz+13; LR13;
Fis+13; Vu14; Mas14; YP14; BLM15; GLM17; PZ17] and convex relaxations via semidefinite programs (SDPs) [AL14; HWX16a; Ban15; GV16; MS16; RTJM16; ABKK17; PW17],
as well as local methods such as belief propagation [DKMZ11], Bayesian MCMC [NS01] and
variational Bayes [CDP12; BCCZ13], greedy profile likelihood [BC09; ZLZ+12] and pseudolikelihood [ACBL+13] maximization, among others. A limitation of spectral clustering approaches is that they are often not optimal on their own, and the SDPs have the drawback of
not being able to fit the full generality of SBMs. Various algorithms can further improve the
clustering accuracy, and adapt to the generality of SBM. Profile likelihood maximization was
proposed and analyzed in [BC09], but the underlying optimization problem is computationally
infeasible and the approach only applicable to networks of limited size. Pseudo-likelihood ideas
were used in [ACBL+13] to derive EM type updates to maximize a surrogate to the likelihood
of the SBM based on a block compression which is computed using initial labels obtained by
spectral clustering.
The pseudo-likelihood approach belongs to the general class of algorithms based on “Good
Initialization followed by Fast Local Updates” (GI-FLU) which has been a staple of recent
developments in devising optimal clustering algorithms, as was pointed out by [GMZZ17]. The
GI-FLU strategies often use a spectral initialization and due to their often cheap local updates
are scalable to very large networks. In cases where they are accompanied by optimality guarantees they seem to occupy the ideal sweet spot in the computational versus statistical trade-off.
We build on these ideas and especially the approach in [ACBL+13] to extend the algorithms to
the bipartite settings. Moreover, we provide modifications to the general blueprint suggested
by [ACBL+13]—in addition to the natural modification required for the bipartite setup—which
allows us to demonstrate optimality of the procedure under the full generality of the bipartite
SBM.
In the univariate setting, there has been interesting recent advancements in understanding
optimal recovery in what we refer to as the semi-sparse regime, where the (expected) average
network degree dav is allowed to grow to infinity but rather slowly, as the number of nodes n
increases to infinity. In a series of papers [MNS15; ABH16; HWX16a; HWX16b] the thresholds
for optimal exact recovery, also known as strong consistency, were established in the context
of simple planted partition models, and SDPs where shown to achieve the optimal threshold.
2

In [AS15], the problem of strong consistency was considered for a general SBM and the optimal threshold for strong consistency was established. In subsequent work [ZZ+16; GMZZ17;
GMZZ16], the results were extended to include weak consistency, i.e., requiring the fraction of
misclassified nodes to go to zero, rather than drop to exactly zero (as in strong consistency),
and rates of optimal convergence where established (with some slack in the exponent). To
achieve the more relaxed consistency results, [GMZZ17] limited the model to what we refer to
as strongly assortative SBM (see [AL18] for a definition).
Our work from a theoretical perspective is mostly inspired by the insightful work of [AS15;
GMZZ17]. We extend these ideas by presenting results that are sharper and more general that
what has been obtained so far. In short, we only assume that the clusters are distinguishable (in
the sense of Chernoff divergence) and the network is not very dense. Relative to [GMZZ17], our
results are shaper, removing the o(1) term in the exponents of the rates, and hold for the full
generality of the SBM (i.e., no assortativity assumption needed). Compared to [AS15], our work
greatly relaxes the assumptions on degree growth. It is well-known that for strong consistency
one needs at least dav = Ω(log n), i.e., average degree should grow at least logarithmically,
and this is the regime considered in [AS15]. In our work, dav could grow to infinity arbitrarily
√
slowly or at the other extreme as fast as dav = O( n). Thus, our results establish optimal rates
√
of weak consistency below the sequare-root regime dav = O( n) and above the sparse regime
dav = O(1), and in particular, between the log n and O(1) regimes, where weak consistency is
possible but not the strong consistency. In addition, in contrast to [AS15] (and some results
in [GMZZ17]), we allow all the parameters of the model, including the number of communities,
the spread of mean parameters (ω) and the community balance parameter (β) to grow subject to
two compact conditions (namely (A3) and (A4)); these conditions encapsulate in a simple way
how much cumulative growth these auxiliary parameters can exhibit relative to the information
col ) for the optimal rate to still be achieved by the algorithm
growth of the model (I = Imin ∧ Imin
we present. We make more detailed comparisons with the work of [AS15; GMZZ17] in Section 3.
Contributions. Establishing these results require a fair amount of technical and algorithmic
novelty over the previous work. Here, we highlight some of these features and point out to the
relevant parts of the paper for details:
1. We introduce an efficient sub-block (or sub-graph) partitioning scheme for the provable version of the algorithm, Algorithm 3, which for example replaces the edge splitting idea
of [AS15] (cf. Section 3.1 and Remark 5 for the shortcomings of edge splitting in our setting).
The subblock partitioning is mainly introduced to generate enough independence for the technical arguments to go through. However, the idea turns out to be computationally appealing
as well. The computational bottleneck of GI-FLU approaches discussed earlier is often the
spectral initialization. The subsequent (often likelihood based) local updates are usually quite
fast, O(n), computations. Our block partitioning scheme allows one to break the costly initial
step into the application of spectral methods—as well as likelihood ratio classifiers—on smaller
subblocks, without losing optimality. If done in parallel, spectral clustering on subblocks will
be in fact computationally cheaper than performing a spectral decomposition of the entire
matrix. Although, from a theoretical perspective breaking into Q = 4 blocks is enough (in the
bipartite settings), our results allow for the number of blocks Q to even grow slowly to infinity.
Thus, the provable version of our algorithm has computational appeal, esp. in distributed settings and for very large networks where it is prohibitive to perform eigendecomposition of the
entire adjacency matrix. The algorithm is naturally parallelizable since it proceeds in stages
3

and in each stage the operations on the underlying subblocks can be performed in parallel; see
Section 5 for details.
2. Our algorithms being an extension of [ACBL+13], are modifications of a natural EM
algorithm on mixtures of Poisson vectors, hence very familiar from a statistical perspective. In
other words, they are not tailor-made to the community detection problem and are derived from
fairly general well-known principals. Although other (optimal) algorithms in the literature are
more or less preforming similar operations, the link to EM algorithms and mixture modeling
is quite clear in our work. We provide in Section 4.2 the general blueprint of the algorithms
based on the pseudo-likelihood idea and block compression (Algorithm 1). We then show how
a provable version can be constructed by combining with the sub-block partitioning ideas in
Section 5. It is worth noting that although we can provide no guarantees for the general
algorithms of Section 4.2, empirically they perform very well, as illustrated in Section 8.
3. In order to get the sharper rate, analyzing a single step of an EM type algorithm is not
enough, and thus we analyze the second step as well. We will show that the first step gets us
from a good (but crude) initial rate γ1 to the fast rate ≈ exp(−I/Q) which is in the vicinity
of the optimal rate, and then repeating the iteration once more, with the more accurate labels
obtained in the first step—hence more accurate parameter estimates—gets us to the optimal
rate ≈ exp(−I).
4. Among the technical contributions, are a uniform consistency result (Lemma 2) for the
likelihood ratio classifier (LRC) over a subset of the parameters close to the truth, sharp approximations for the Poisson-binomial distributions (Section 9.4), and extension (and elucidation)
of a novel technique of [AS15] in deriving error exponents to general exponential families (cf.
Section 9.3). The uniform consistency result for LRCs lets us tolerate some degree of dependence among the statistics from iteration to iteration, allowing the sub-block partioning idea to
go through. That is, we can run LRCs on the same blocks used in estimating their parameters;
see Sections 5 and 6.1 for details.
5. The bipartite clustering setup (as opposed to the symmetric case) allows us to introduce
an oracle version of the problem which helps in understanding the nature of the optimal rates
observed in community detection and their relation to classical hypothesis testing and mixuture
modeling. That is, we try to answer the curious question of why or how the Chernoff exponent
of a (simple binary) hypothesis testing problem seems to control the misclassification rate in
community detection and network clustering. The oracle also provides a lower bound on the
performance of any algorithm. See Section 1 and Proposition 1 for details.
The rest of the paper is organized as follows. We introduce the model and the biclustering
oracle in Section 2, and then present our main results in Section 3. The general algorithms based
on the pseduo-likelihood idea are presented in Section 4 and a provable version in Section 5. The
proofs of the consistency results will appear in Sections 6 and 7. In Section 8, we demonstrate
the numerical performance of the methods.

2

Network biclustering

We start by introducing the network biclustering problem based a stochastic block modeling,
and set up some notation that will be used throughout the paper. We then discuss how a

4

biclustering oracle with side information can optimally recover the labels. These ideas will be
the basis of the algorithms discussed in this paper.

2.1

Bipartite block model

We will be working with a bipartite network which can be represented by a biadjacency matrix
A ∈ {0, 1}N ×M , where for simplicity we assume that the nodes on the two sides are indexed
by the sets [N ] and [M ]. We assume that there are K and L communities for the two sides
respectively, and the membership of the nodes to these communities are given by two vectors
y = (yi ) ∈ [K]N and z = (zj ) ∈ [L]M . Thus, yi = k if node i on side 1 belongs to community
k ∈ [K]. We call yi and zj the labels of nodes i and j respectively. We often treat these labels
as binary vectors as well, using the identification [K] ' {0, 1}K via the one-hot encoding, that
is yi = k ⇐⇒ yik = 1, yik0 = 0, k 0 6= k.
Given the labels y and z, and a connectivity matrix P ∈ [0, 1]K×L (also known as the edge
probability matrix), the general bipartite stochastic block model (biSBM) assumes that: Aij
are Bernoulli variables, independent over (i, j) ∈ [N ] × [M ] with mean parameters,
E[Aij ] = yiT P zj = Pk` ,

if yi = k, zj = `.

(1)

We denote this model compactly as A ∼ SBM(y, z, P ). It is helpful to consider the Poisson
version of the model as well which is denoted as A ∼ pSBM(y, z, P ). This is the same model
as the Bernoulli SBM, with the exception that each entry Aij is drawn (independently) from a
Poisson variate with mean given in (1). These two models behave very closely when the entries
of P are small enough. Throughout, we treat y, z and P as unknown deterministic parameters.
The goal of network biclustering is to recover these three sets of parameters given an instance
of A.
In fact, as we will see, the parameters P themselves are not that important. What matters is
the set of (Poisson) mean parameters which are derived from P and the sizes of the communities.
In order to define these parameters, let n(z) = (n1 (z), . . P
. , nL (z)) ∈ NL , be the
PMnumber of nodes
1{z
=
`}
=
in each of the communities of side 2. That is, n` (z) = M
j
j=1 zj` . A similar
j=1
K
notation, namely n(y) ∈ N , denotes the community sizes of side 2. The row mean parameters
are defined as
Λ = (λk` ) = (Pk` n` (z)) = P diag(n(z)) ∈ RK×L

(2)

where diag(v) for a vector v = (vk ) is a diagonal matrix with diagonal entries vk . The column
mean parameters can be defined in a similar fashion, namely,

ΓT = nk (y)Pk` = diag(n(y))P ∈ RK×L .
(3)
Note the transpose in the above definition, i.e., Γ ∈ RL×K , and this convention allows us to
define information measures based on rows of matrices Λ and Γ in a similar fashion, as will be
discussed in Section 3.

2.2

Biclustering oracle with side information

The key idea behind the algorithms discussed in this paper, as well as our consistency arguments
is the following simple observation: Assume that we have prior knowledge of P and the column
5

labels z, but not the row labels y. For each row, we can sum the columns of A according
to
P
their column memberships, i.e., we can perform the (ideal) block compression b∗i` := j Aij zj` .
The vector b∗i∗ = (b∗i1 , . . . , b∗iL ) contains the same information for recovering the community
of i, as the original matrix A—i.e., it is a sufficient statistic. Assume that we are under the
pSBM model. Then, b∗i∗ has the distribution of a vector of independent Poisson variables. More
precisely,
b∗i∗ ∼ Qk :=

L
Y

Poi(λk` ),

if,

yi = k,

(4)

`=1

where λk` are the row mean parameters defined in (2). Note that the distributions Qk , k =
1, . . . , K are known under our simplifying assumptions. The problem of determining the row
labels thus reduces to deciding from which of these K known distributions it comes from.
Whether node i belongs to a particular community k can be decided optimally by performing
a likelihood ratio (LR) test of Qk against each of Qr , r 6= k.
The above LR test is the heart of the algorithms discussed in Sections 4 and 5. The
difficulty of the biclustering problem (relative to a simple mixture modeling) is that in practice,
we do not know in advance either y or Λ—hence neither the exact test statistics (b∗i∗ ) nor the
distributions {Qk } are known. We thus proceed by a natural iterative procedure: Based on the
initial estimates of y and z, we obtain estimates of (b∗i∗ ) and {Qk }, perform the approximate
LR test to obtain better estimates of z, and then repeat the procedure over the columns to
obtain better estimates of y. These new label estimates lead to better estimates of (b∗i∗ ) and
{Qk }, and we can repeat the process.
We refer to the algorithm that has access to the true column labels z and parameters Λ,
and performs the optimal LR tests, as the oracle classifier. Note that the performance of this
oracle gives a lower bound on the performance of any biclustering algorithm in our model. The
performance of the oracle in turn is controlled by the error exponent of the simple hypothesis
testing problems Qk versus Qr , r 6= k, as detailed in Proposition 1. This line of reasoning reveals
col —defined in (8) and (9)—that control the
the origin of the information quantities Ikr and I`r
optimal rate of the biclustering problem. Note that the bipartite setup has the advantage
of disentangling the row and column labels, so that a non-trivial oracle exists. It does not
make much sense to assume known column labels in the unipartite SBM, since by symmetry
we then know the row labels as well, hence nothing left to estimate. On the other hand, due
to the close relation between the bipartite and unipartite problems, the above argument also
sheds light on why the error exponent of a hypothesis test is the key factor controlling optimal
misclassification rates of community detection in unipartite SBM.

2.3

Notation on misclassification rates

Let Πn the set of permutations on [n]. The (average) misclassification rate between two sets
of (column) labels yb and y is given by
n
1X 
Mis(b
y , y) := min
1 σ(b
yi ) 6= yi .
σ ∈ Πn n
i=1

6

(5)

Letting σ ∗ be a minimizer in (5), the misclassification rate over cluster k is
X 
|i : σ ∗ (b
yi ) 6= k, yi = k|
1
1 σ ∗ (b
yi ) 6= yi =
,
Misk (b
y , y) :=
nk (y)
nk (y)

(6)

i:yi =k

using the cardinality notation to be discussed shortly. Note that (6) is not symmetric in its
arguments. We will also use the notation σ ∗ (b
y → y) to denote an optimal permutation in (5).
When Mis(b
y , y) is sufficiently small, this optimal permutation will be unique. It is also useful
to define the direct misclassification rate between yb and y, denoted as dMis(b
y , y), which is
∗
∗
obtained by setting the permutation in (5) to the identity. With σ = σ (b
y → y), we have
Mis(b
y , y) = dMis(σ ∗ (b
y ), y). We note that
X
Mis(b
y , y) =
πk (y) Misk (b
y , y) ≤ max Misk (b
y , y),
(7)
k∈K

k∈K

as well as maxk∈K Misk (b
y , y) ≤ Mis(b
y , y)/ mink0 πk0 (y). We can similarly define the misclassification rate of an estimate zb relative to z. Our goal is to derive efficient algorithms to obtain yb
and zb that have minimal misclassification rates asymptotically (as the number of nodes grow).

Other notation. We write w.h.p. as an abbreviation for “with high probability”, meaning
that the event holds with probability 1 − o(1). To avoid ambiguity, we assume all parameters,
including m, are functions of n. For example, f (n) = o(g(n)) denotes limn→∞ f (n)/g(n) = 0.
We write ZQ = Z/QZ to denote a cyclic group of order Q. Our convention regarding solutions
of optimization problems, whenever more than one exist is to choose one uniformly at random.
We use the shorthand notation |i : yi = k| := |{i : yi = k}| for cardinality of sets, where i ∈ [n]
n
is implicit, assuming the
P y is a vector of length n. For example, if yb, y ∈ [K] , we have the
identity |i : ybi 6= yi | = k∈[K] |i : yi = k, ybi 6= k|. It is worth nothing that we use community
and cluster interchangeably in this paper, although some authors prefer to use community for
the assorative clusters, and use “cluster” to refer to any general group of nodes. We will not
follow this convention and no assortativity will be implicitly assumed.

3

Main results

Let us start with some assumptions on the mean parameters. Recall the row and column mean
parameter matrices Λ and Γ defined in (2) and (3). Let Λmin and kΛk∞ be the minimum and
maximum value of the entries of Λ, respectively, and similarly for Γ. We assume
kΛk∞ kΓk∞
∨
≤ ω,
Λmin
Γmin

(A1)

for some ω > 0. That is, ω measures the deviation of the entries of the mean matrices from
uniform. We assume that the sizes of the clusters are bounded as
1
β
1
β
≤ πk (y) ≤
and
≤ π` (z) ≤
(A2)
βK
K
βL
L
for all k ∈ [K] and ` ∈ [L]. The following key quantity controls the misclassification rate:
Ikr := Ikr (Λ) := sup

L
X

s ∈ (0,1) `=1

s
(1 − s)λk` + sλr` − λ1−s
k` λr` ,

7

(8)

for k, r ∈ [K]. We can think of I(Λ) := (Ikr (Λ)) ∈ RK×K
, as an operator acting on pairs of
+
rows of a matrix Λ ∈ RK×L
,
say
λ
and
λ
,
producing
a
K
× K pairwise information matrix.
r∗
k∗
+
We often refer to the function of s being maximized in (8) as s 7→ Is , with some abuse of
notation assuming k and r are fixed, and we note that this function is strictly concave over R
whenever λk∗ 6= λr∗ , and we have I0 = I1 = 0.
Recalling the product Poisson distributions {Qk }, (8) is the Chernoff exponent in testing
the two hypothesis Qk and Qr [Che52]. The difference with the classical setting, in which the
Chernoff exponent appears, is that we work in the regime where we are effectively testing based
on a sample of size of 1 and instead of the sample size, we let Ikr → ∞. We define the column
information matrix similarly
col
I``
0 := I``0 (Γ) = sup

K
X

s∈(0,1) k=1

1−s s
(1 − s)Γ`k + sΓ`0 k − Γ`k
Γ`0 k ,

(9)

for all `, `0 ∈ [L]. Another set of key quantities in our analysis are:
εkr := max
`∈[L]

λ

k`

λr`

∨

λr` 
− 1,
λk`

εk := min εkr ,
r∈[K]

and

ε := min εk .

(10)

k∈[K]

The relation with hypothesis testing is formalized in the following proposition:
Proposition 1. Consider the likelihood ratio (LR) testing of the null hypothesis Qk against
2×K
Qr , based on a sample of size 1. Let Λ = [λk∗ ; λr∗ ] ∈ R+
. Assume that as Λmin → ∞, (a)
lim inf εkr > 0, and (b) ω = O(1). Then, there exist constants C and C 0 such that


≤ C exp −Ikr − 1 log Λmin ,
2
P(Type I error) + P(Type II error)
(11)

≥
exp −Ikr − L2 (log Λmin + C 0 ) .

See Corollary 6 and Appendix Section A.6 for the proof. Any hypothesis testing procedure
can be turned into a classifier, and a bound on the error of the hypothesis test (for a sample
of size 1) translates into a bound on the misclassification rate for the associated classifier.
This might not be immediately obvious, and we provide a formal statement in Lemma 3.
Proposition 1 thus provides a precise bound on the misclassification rate of the LR classifier
for deciding between Qk and Qr .
The significance of the Chernoff exponent of the hypothesis test in controlling the rates
is thus natural, given the full information about the {Qk } and the test statistics. What is
somewhat surprising is that almost the same bound holds when no such information is available
a priori. Our main result below is a formalization of this claim. In our assumptions, we include
a parameter Q ∈ N that controls the number of sub-blocks when partitioning, the details of
which are discussed in Section 5. Under the following two assumptions:
(Q2 log Q)β 2 ω 3 KL(K ∨ L) log(K ∨ L)(kΛk∞ ∨ kΓk∞ )2 = O(n ∧ m), and
2 3 2

3

(Q log Q) β ω (K ∨ L) (α ∨ α

−1

)(kΛk∞ ∨ kΓk∞ ) = o (Imin ∧

col 2
Imin
)



(A3)
,

where α := m/n, there is an algorithm that achieves almost the same rate as the oracle:

8

(A4)

Theorem 1 (Main result). Consider a bipartite SBM (Section 2.1) satisfying (A1)–(A4).
col → ∞ and Λ
Then, as Imin ∧ Imin
b output by Algorithm 3 in Section 5
min → ∞, the row labels y
satisfies for some ζ = o(1),

1


 X

1 
exp −Ikr −
− ζ log Λmin
Misk yb, y = O ω
1+
εkr
2

(12)

r6=k

for every k ∈ [K], with high probability. Similar bounds holds for zb w.r.t. z.

One can replace the big O with the small o in (12) to obtain an equivalent result (due to the
presence of ζ = o(1)). Let us discuss the assumptions of Theorem 1. The only real assumptions
are (A3) and (A4). The other two, namely (A1) and (A2) can be more or less thought of as
definitions of ω and β. For example, (A2) only imposes the mild constraint that no cluster is
empty. Similarly (A1) imposes the mild assumption that no entry of Λ or Γ is zero. The main
constraints on ω and β are encoded in (A3) and (A4) in tandem with other parameters of the
model.

Remark 1. In the first reading, one can take β, ω, Q = O(1), n  m and kΛk∞  kΓk∞ .
In this setting, (A3) is a very mild sparsity condition, implying that the degrees should not
√
grow faster than n. (A4) guarantees that the information quantities grow fast enough so
that the clusters are distinguishable. We only need (A4) for Algorithm 3 which uses a spectral
initialization. In Section 5.2.1, we present Theorem 2, for the likelihood-based portion of the
algorithm, assuming that a good initialization is provided. Theorem 2 only requires a weakened
version, (A40 ), of assumption (A4).
Depending on the behavior of εkr , the rate obtained in Theorem 1 can exhibit different
regimes which are summarized in Corollary 1 below. Consider the additional assumption:

1 
max ω 1 +
= O(1).
(A5)
εkr
k,r ∈ [K]
Corollary 1. Under the same assumptions as Theorem 1, w.h.p., for all k ∈ [K],
X


Misk yb, y = o
exp(−Ikr ) .

(13)

r6=k

If in addition we assume (A5), then for some ζ = o(1), w.h.p., for all k ∈ [K],
X

1



exp −Ikr −
Misk yb, y = O
− ζ log Λmin .
2

(14)

r6=k

Remark 2. Consider the oracle version of the biclustering problem where the connectivity
matrix P and the true column labels z are given. Then, the optimal row clustering reduces to
the likelihood ratio tests in Proposition 1. That is, given the row sums within blocks as sufficient
statistics, we compare the likelihoods at two different mean parameters. By Proposition 1, the
optimal misclassification rate for the oracle problem is
X




1
E Misk yb, y = O
exp −Ikr − log Λmin ,
(15)
2
r6=k

9

where the sum over r is due to the need to compare against all other clusters. The gap
between 1/2 and 1/2 − ζ is not avoidable when stating high probability results, due to the
Markov inequality; see Lemma 3 for the details. This error rate coincides with (14), which
merely loses a constant due to the unknown mean parameters and column labels. The rate is
sharp up to a factor of exp(− 12 (L − 1) log Λmin ) according to the lower bound in Proposition 1.
In order to understand the rates in Corollary 1 better, let us consider some examples which
also clarify our results relative to the previous literature.
Example 1. Consider a simple planted partition model where
n = m,

K = L,

Pkk =

a
,
n

Pk` =

b
, ∀k 6= `.
n

a βa
b βb
Then, λkk ∈ [ βK
, K ] and λk` ∈ [ βK
, K ] when k 6= `. Applying (8) with s = 1/2,
√
√
p
1X p
( a − b)2
2
Ikr ≥
( λk` − λr` ) ≥
.
2
βK
`

Assume that (A3) and (A4) hold, that is (using kΛk∞ ≤ βa/K)
β 4 ω 3 (K log K)a2 = O(n ∧ m)

√
√ 
and β 6 ω 2 K 4 a = o ( a − b)4 .

and further assume that βω 2 K 3 = o(b). Then w.h.p., we have

 (√a − √b)2 

.
Misk yb, y = o exp −
βK

(16)

For the details of (16), see Section 7.4. In particular, if
√
√
( a − b)2
lim inf
≥ 1,
(17)
n→∞
βK log n

we have Misk yb, y = o(1/n) w.h.p., that is, we have the exact recovery of the labels by
Algorithm 3. (Whenever misclassification rate drops below 1/n, it should be exactly zero.)
Note that this result holds without any assumption of assortativity, i.e., it holds whether a > b
or b > a.
Example 2. Suppose that P := Pe(log n)/n where Pe is a symmetric constant matrix, n = m,
K = L, and y = z. Then K, ω and εkr are constants. Then,
ek` log n,
λk` = λ

ek` := Pek` πk (y),
where λ

and Ikr = I˜kr log n

ek` and λ
er` as in (8). Assuming in addition that π(y) is constant,
where I˜kr is defined based on λ
ekr and I˜kr are constants. In this regime, our key assumptions (A3) and (A4) are satisfied.
both λ
By Corollary 1, w.h.p., we have






˜
Misk yb, y = o exp − min I˜kr log n = o n− minr6=k Ikr .
r6=k

(18)


As a consequence if mink6=r I˜kr ≥ 1, then Misk yb, y = o(1/n) w.h.p., that is we have exact
recovery by Algorithm 3.
10

3.1

Comparison with existing results

Let us now compare with [GMZZ17] and [AS15] whose results are closest to our work. Both
papers consider the symmetric (unipartite) SBM, but the results can be argued to hold in the
bipartite setting as well. The setup of Example 1 is more or less what is considered in [GMZZ17].
They have shown that there is an algorithm with misclassification error bounded by
 (1 − o(1))(√a − √b)2 
exp −
.
(19)
βK
We have sharpened this rate to (16) under assumption (A3) (i.e., assuming the average degree
√
grows slower than O( n)). Bound (19) implies that when
√
√
( a − b)2
lim inf
> 1,
n→∞
βK log n
one has exact recovery. Our bound on the other hand, imposes the relaxed condition (17).
We note that the results in [GMZZ17] are derived for a more general class of (assortative)
models than that of Example 1, namely, the class with connectivity matrix satisfying Pkk ≥ a/n
and Pk` ≤ b/n for k 6= `. The rate obtained in [GMZZ17] uniformly over this class is dominated
by that of the hardest within this class which is the model of Example 1. For other members
of this class, neither their rate (19) or the one we gave in (16) is optimal. The optimal rate in
those cases is given by the general form of Theorem 1 and is controlled by the general form of
Ikr in (8). In other words, Algorithm 3 that we present is rate adaptive over the class considered
in [GMZZ17], achieving the optimal rate simultaneously for each member of the class.
A key in our approach is to apply the likelihood-type algorithm twice, in contrast to the
single application in [GMZZ17]. After the second stage we obtain much better estimates of
the labels and parameters relative to the initial values, allowing us to establish the sharper
forms of the bounds. Another key is the result in Lemma 2(b) which provides a better error
rate than the classical Chernoff bound, using a very innovative technique introduced in [AS15].
Moreover, we keep track of the balance parameter β in (A2) throughout, allowing it to go to
infinity slowly. Last but not least, assortativity is a key assumption in [GMZZ17], while our
result does not rely on it. Besides consistency, our provable algorithm is more computationally
efficient in a practical sense. To obtain initial labels, we will apply spectral clustering on very
few subgraphs (8 to be exact). However, the provable version of the algorithm in [GMZZ17]
applies spectral clustering for each single node on the rest of the graph excluding that node. If
the cost of running the spectral clustering on a network of n nodes is Cn , then our approach costs
≈ 8Cn/8 while that of [GMZZ17] costs roughly nCn−1 . Our algorithm thus has a significant
advantage in computational complexity when n → ∞. To be fair, the algorithm introduced
in [GMZZ17] was for the symmetric SBM, which has the extra complication of dependency
in A due to symmetry. Our comparison here is mostly with Corollary 3.1 in [GMZZ17]. In
addition, [GMZZ17] have a result (their Theorem 5) for when ω grows arbitrarily fast which is
not covered by our result. See the discussion below for comments on the symmetric case and
dependence on ω.
The problem of exactly recovery for a general SBM has been considered in [AS15], again for
the case of a symmetric SBM, though the results are applicable to the bipartite setting (with
y = z) as well. The model and scaling considered in [AS15] is the same as that of Example 2,
and they show that exact recovery of all labels is possible if (and only if) mink,r:k6=r I˜kr ≥ 1
11

which is the same result we obtain in Example 2 for Algorithm 3. Thus, our result contains that
of [AS15] as a special case, namely in the log n-degree regime with other parameters (such as K
and the normalized connectivity matrix) kept constant. The results and algorithms of [AS15]
do not apply to the general model in our paper; consider the following two points:
1. Only the regime P ∼ log n/n, i.e., the degree grows as fast as log n, is investigated
in [AS15], while we allow the degree to grow in the range from “arbitrarily slowly” up to “as
√
fast as O( n)”.
2. One needs independent versions of the adjacency matrix in different stages of the algorithm. To achieve this goal, edge splitting was introduced in [AS15]. The idea is that one
can regard the two (or more) graphs obtained from edge splitting to be nearly independent.
To be specific, let P1 be the joint probability measure corresponding to a pair of graphs G1
and G2 generated independently with connectivity matrices qP and (1 − q)P . Let P2 be the
joint probability measure on G1 and G2 obtained by edge splitting from a single SBM with
connectivity matrix P , assigning every edge independently to either G1 or G2 with probabilities q and 1 − q. Then, P1 and P2 have the same marginal distributions. Having a vanishing
total variation between P1 and P2 is necessary for further analysis which, as was pointed out
by [AS15, pp. 46-47], is equivalent to showing that under P1 , G1 and G2 do no share any edge,
with high probability. Letting Pemin = mink` Pek` ,

2 (log n)2 n2
(1 − q)q Pemin
P1 (G1 and G2 do not share edges) ≤ 1 −
n2

2 (log n)2 = o(1), that is, the connecwhich is strictly bounded away from 1 unless (1 − q)q Pemin
tivity matrix of either G1 or G2 should vanish faster than 1/n. Our consistency result will not
hold in this regime. Thus, edge splitting cannot be used to derive the results in this paper, and
we introduce the block partitioning idea to supply us with the independent copies necessary
for analysis. Another technical issue about edge splitting is discussed in Remark 5.

3.2

Discussion

Our results do not directly apply to the symmetric case, due to the dependence between the
upper and lower triangular parts of the adjacency matrix A. However, a more sophisticated
two-stage block partitioning scheme can be used to derive similar bounds under mild extra
assumptions. One starts with an asymmetric partition into blocks of sizes {qn, (1 − q)n} ×
{qn, (1−q)n}, for q = 1/Q → 0 very slowly. In the first stage, one applies a similar procedure as
described in Algorithm 3 on the upper triangular portion of the large subblock (1−q)n×(1−q)n,
followed by the application of the LR classifier on the fat block qn × (1 − q)n to obtain very
accurate row labels of the small block qn × qn.. One then repeats the process using the “leaveone-out” of [GMZZ17], but applied to small blocks qn × qn rather than individual nodes. We
leave the details for a future work.
It was also shown by [GMZZ17, Theorem 5] that their equivalent of condition (A1) can be
removed by modifying the algorithm. In their setting, without assuming a  b, a misclassification rate of exp(−(1−ε)I) is achievable, where ε ∈ (0, 1) is a variable in the new version of their
algorithm. If those arguments can be extended to the general block model, it will be possible to
relax the requirements on ω in (A3) and (A4). When K, L = O(1), one can completely remove
sparsity condition (A3) using a much sharper Poisson-binomial approximation than what we
have used in this paper. Finally, we suspect that our result could be generalized beyond SBMs
12

to biclustering arrays where the row and column sums over clusters follow Poissonian central
limit theorems. We will explore these ideas in the future.

4

Pseudo-likelihood approach

In this section, after introducing the local and global mean parameters which will be used
throughout the paper, we present our general pseudo-likelihood approach to biclustering.

4.1

Local and global mean parameters

Let us define the following operator that takes an adjacency matrix A and row and column
labels ỹ and z̃, and outputs the corresponding (unbiased) estimate of its mean parameters:
n

[L (A, ỹ, z̃)]k` =

m

1 XX
Aij 1{ỹi = k, z̃j = `},
nk (ỹ)
i=1 j=1

k ∈ [K], ` ∈ [L].

(20)

Note that L (A, ỹ, z̃) is a K × L matrix with nonnegative entries. In general, we let
Λ̂ = (λ̂k` ) := L (A, ỹ, z̃),

(21)

Λ(ỹ, z̃) = (λk` (ỹ, z̃)) := L (E[A], ỹ, z̃),

(22)

for any row and column labels ỹ and z̃. Here Λ̂ is the estimate of the true row mean matrix.
Its expectation is E[Λ̂] = Λ(ỹ, z̃) due to the linearity of L . We call Λ(ỹ, z̃), the (global) row
mean parameters associated with labels ỹ and z̃. (We do not explicitly show the dependence
of Λ̂ on the labels, in contrast to the mean parameters.) We have the following key identity
Λ(ỹ, z̃) |ỹ=y, z̃=z = Λ

(23)

where Λ is the true (global) row mean parameter matrix defined in (2). In words, (23) states
that the global row mean parameters associated with the true labels y and z, are the true such
parameters. We will also use parameters such as Λ(y, z̃) which are obtained based on the true
row labels y and generic column labels z̃.
We also need local versions of all these definitions which are obtained based on submatrices
0
0
of A. More precisely, let A(q ,q) be a submatrix of A, and let y (q ) and z (q) be the corresponding
subvectors of z and y (i.e., corresponding to the same row and column index sets used to
extract the submatrix). Here q, q 0 range over [Q] = {1, . . . , Q} creating a partition of A into
Q2 subblocks. We call
(q 0 ,q)

0

Λ(q ,q) (ỹ, z̃) := (λk`

0

0

(ỹ, z̃)) := L (E[A(q ,q) ], ỹ (q ) , z̃ (q) ),
0

(24)
0

the local row mean parameters associated with submatrix A(q ,q) and sublabels y (q ) and z (q) .
0
0
The corresponding estimates are defined similarly (by replacing E[A(q ,q) ] with A(q ,q) ). We will
0
mostly work with submatrices obtained from a partition A(q ,q) , q 0 , q ∈ [Q] of A into (nearly)
equal-sized blocks—the details of which are described in Section 5. In such cases,
0

Λ(q ,q) (ỹ, z̃) ≈

1
Λ(ỹ, z̃),
Q
13

∀q 0 , q ∈ [Q]

assuming the each subblock in the partition has nearly similar cluster proportions: n(z (q) ) ≈
n(z). This is the case, for example, for a random partition as we show in Section 6.2. Of special
0
interest is when we replace both ỹ and z̃ with true labels y and z. In such cases, Λ(q ,q) does
not depend on q 0 . More precisely, we have for any q ∈ [Q],
(q 0 ,q)

λk`

(y, z) = Pk` n` (z (q) ),

∀q 0 ∈ [Q],

(25)

where n` (z (q) ) is the number of labels in class ` in z (q) , consistent with our notation for the
0
(q)
full label vectors. We often write Λ(q) = (λk` ) as a shorthand for Λ(q ,q) (y, z) which is justified
by the above discussion. These will be called the true local row mean parameters (associated
with column q subblock in the partition).

4.2

General pseudo-likelihood algorithm

Let us now describe our main algorithm based on the pseudo-likelihood (PL) idea, which is
a generalization of the approach in [ACBL+13] to the bipartite setup. The pseudo-likelihood
algorithm (PLA) is effectively an EM algorithm applied to the approximate mixture of Poissons
obtained from the block compression of the adjacency matrix A. It relies on some initial
estimates of the row and column labels to perform the first block compressions (for both rows
and columns). The initialization is often done by spectral clustering and will be discussed
in Section 5.2.2 once we introduce the provable version of the algorithm. Subsequent block
compressions are performed based on the label updates at previous steps of PLA.
Let us assume that we have obtained labels ỹ and z̃ as estimates of the true labels y and z.
We focus on the steps of PLA for recovering the row labels. Let us define an operator B(A; z̃)
that takes approximate columns labels and produces the corresponding column compression
of A:
m
X

n×L
B(A; z̃) := b(z̃) := bi` (z̃) ∈ Z+
, bi` (z̃) :=
Aij 1{z̃j = `}.
(26)
j=1

The distribution of bi` (z̃) is determined by the row class of i. It is not hard to see that
E[bi` (z̃)] = λk` (y, z̃) = λk` (y, z̃)|z̃=z ,

if yi = k,

(27)

where λk` (y, z̃) are the (global) row mean parameters defined in (22).
Now consider an operator L(b; ỹ) that given the column compression b and the initial
estimate of the row labels ỹ, produces estimates of the (row) mean parameters λk` (y, z̃):
n

L(b; ỹ) := Λ̂ := [λ̂k` ] ∈ RK×L
,
+

λ̂k` :=

1 X
bi` 1{ỹi = k}.
nk (ỹ)

(28)

i=1

Note that if ỹ = y, we have E[λ̂k` ] = λk` (y, z̃). The definition of the estimates in (28) are
consistent with those of (21) due to the following identity:
L(B(A; z̃); ỹ) = L (A, ỹ, z̃)
which holds for any row labels ỹ and column labels z̃. Let us write
n

π(ỹ) := (πk (ỹ)),

πk (ỹ) :=

1X
1{ỹi = k}
n
i=1

14

(29)

Algorithm 1 Pseudo-likelihood biclustering, meta algorithm
1: Initialize row and column labels ỹ and z̃.
2: while ỹ and z̃ have not converged do
3:
b ← B(A; z̃)
4:
while Λ̂ and π̂ not converged (optional) do
5:
Λ̂ ← L(b; ỹ)
6:
Option 1: π̃ ← 1, or option 2: π̃ ← π(ỹ)
7:
ỹ ← F(b, Λ̂, π̃)
8:
(Optional) Convert ỹ to hard labels.
9:
end while
10:
Repeat lines 3–7 with appropriate modifications to update z̃ and columns parameters
(by changing A to AT and swapping z̃ and ỹ.)
11:
(Optional) Convert ỹ and z̃ to hard labels if they are not.
12: end while
for the estimate of (row) class priors based on ỹ. We note that the operation B and L remain
valid even if ỹ and z̃ are soft labels with a minor modification. By a soft row label z̃j ∈ [0, 1]L
P
we mean a probability vector on [L]: z̃j` ≥ 0 and L
`=1 z̃j` = 1, which denotes a soft assignment
to each row cluster. To extend (26) to soft row labels, it is enough to replace 1{zj = `} with
zj` . Extending (28) to soft column labels ỹ is done similarly.
Now, given any block compression b = (bi` ) and any estimate Λ̂ of the (row) mean parameters and any estimate π̃ ∈ [0, 1]K of the (row) class prior, consider the operator
that
Q outputs
P
the (row) class posterior assuming that the rows of bi approximately follow k π̃k ` Poi(λ̂k` ):
n×K

F(b, Λ̂, π̃) := (π̂ik ) ∈ [0, 1]

,

QL

`=1 ϕ(bi` , λ̂k` )
QL
`=1 ϕ(bi` , λ̂k0 ` )
k0 =1 π̃k0

π̂ik := PK

π̃k

(30)

where ϕ(x, λ) = exp(x log λ − λ) is the Poisson likelihood (up to constants). In practice, we
only use π(ỹ) or a flat prior 1 as the estimated prior π̃ in this step; similarly, we only use a block
compression which is based on estimates of row labels, i.e., bi` = bi` (z̃) for some z̃ ∈ [n]L . Note
that F outputs soft-labels which can be considered our new estimates of y. We can convert
(π̂ik ) to hard labels if needed.
Algorithm 1 summarizes the general blueprint of PLA, which proceeds by iterating the
three operators (26), (28) and (30). Optional conversion from soft to hard labels is performed
by MAP assignment per row. With option 2 in step 6, the inner loop on lines 4–9 is the EM
algorithm for a mixture of Poisson vectors. We can also remove the inner loop and perform
iterations 5–8 only once. In total, Algorithm 1 has (at least) 6 possible versions, depending on
whether we include each of the steps 8 or 11 (for the soft to hard label conversion) and whether
to implement the inner loop till convergence or only for one step. We provide empirical results
for two of these versions in Section 8. In practice, we recommend to keep soft labels throughout,
and only run the inner loop for a few iterations (maybe even one if the computational cost is
of significance).
Remark 3 (PL naming). We have borrowed the name pseudo-likelihood (PL) from [ACBL+13]
based on which the algorithms in this paper are derived. In [ACBL+13], the setup is that of the
symmetric SBM, and in order to treat the full likelihood as the product of independent (over
15

Algorithm 2 Simplified pseudo-likelihood clustering
1:
2:
3:
4:
5:

Input: Initial column labels z̃, and Λ̃ that estimates Λ.
Output: Estimate of row labels yb.
b ← B(A; z̃)
yb ← F(b, Λ̃, 1)
Convert yb to hard labels, by computing MAP estimates.

nodes i = 1, . . . , n) of the mixture of Poisson vectors, one has to ignore the dependence among
the upper and lower triangular parts of the adjacency matrix, making the PL naming more
inline with the traditional use of the term. In our bipartite setup, there is no such dependence
to ignore, but we have kept the name PL for consistency with [ACBL+13] and ease of use.
We interpret the “pseudo” nature of the likelihood as the approximation used in the block
compression stage (with imperfect labels) and in replacing Poisson-binomial distribution with
the Poisson.

4.3

Likelihood ratio classifier

A basic simplified building block of the PLA is given in Algorithm 2. This operation—which
will play a key role in the development of the provable version of the algorithm in Section 5—can
be equivalently described as a likelihood ratio classifier (LRC). Let us write the joint Poisson
likelihood (up to a constant) as:
Φ(x, λ) =

L
Y

ϕ(x` , λ` ) =

L
Y

`=1

`=1

exp(x` log λ` − λ` ),

x ∈ RL , λ ∈ RL
+,

(31)

x ∈ RL , λ, λ0 ∈ RL
+.

(32)

and the corresponding likelihood ratio as:
L

Ψ(x; λ | λ0 ) = log

X
Φ(x, λ)
λ`
x` log 0 + λ0` − λ` ,
=
0
Φ(x, λ )
λ`
`=1

Recalling the column compression (26), the likelihood ratio classifier, based on initial row labels
z̃ and an estimate Λ̃ of the row mean parameter matrix, is
[LR(A, Λ̃, z̃)]i ∈ argmax log Φ(bi∗ (z̃), λ̃r∗ ),
r ∈ [K]

i ∈ [n].

(33)

which gives us a refined estimate of the row labels (i.e., y). It is not hard to see that the output
of Algorithm 2 is yb = LR(A, Λ̃, z̃)].

5

Provable version

When analyzing Algorithm 2, we need the initial labels to be independent of the adjacency
matrix. Hence, we cannot apply the initialization method (e.g., the spectral clustering) and
the likelihood ratio classifier (Algorithm 2) on the same adjacency matrix A, iteratively. In
this section, we introduce an algorithm, namely Algorithm 3, that partitions A into sub-blocks
and operates iteratively on collections of these blocks to maintain the desired independence.
For this version of the pseudo-likelihood algorithm, our main result, Theorem 1, holds.
16

Algorithm 3 Provable version
1:

Randomly partition the rows into 2 groups of equal size (n/2), so that
A = [Atop ; Abottom ]

2:

Randomly partition the rows and columns of Atop into 4 groups of equal size, so that we
have 16 sub-adjacency matrix with dimension (n/8) × (m/4), i.e.

Abottom


A(1,1)

A(2,1)

=
A(3,1)

A(4,1)

A(1,2) A(1,3) A(1,4)




A(2,2) A(2,3) A(2,4) 

.
(3,2)
(3,3)
(3,4)

A
A
A

A(4,2) A(4,3) A(4,4)

In each of the following steps, perform the stated operation for every q ∈ Z4 :


[ỹ (q−1) ; ỹ 0(q) ] ← rowSC [A(q−1,q) ; A(q,q) ] , ∀q.

3:

Obtain initial row labels:

4:

Obtain initial column labels:

5:

Get consistent (global) labels:

6:

Update (local) row mean parameters:

7:

Update row labels:

8:

Similarly update column labels z̃ as in steps 6 and 7.

9:

Update (local) row mean parameters:

10:

Obtain (global) row mean parameters:

13:
14:

ỹ ← Match(ỹ, ỹ 0 ) and z̃ ← Match(z̃, z̃ 0 ).
Λ̂(q+2) ← L (A(q,q+2) , ỹ (q) , z̃ (q+2) ), ∀q.

ỹ (q) ← LR(A(q,q+2) , Λ̂(q+2) , z̃ (q+2) ), ∀q.
Λ̂(q+3) ← L (A(q,q+3) , ỹ (q) , z̃ (q+3) ), ∀q.
P
Λ̂ ← q Λ̂(q) .
ybtop ← LR(Atop , Λ̂, z̃).

11:
12:

[z̃ (q) z̃ 0(q+1) ] ← colSC([A(q,q) A(q,q+1) ]), ∀q.

Swap Atop and Abottom , then repeat steps 2–9 to obtain ybbottom .

Apply step 1 to 10 on

AT

to obtain zb.

yb ← [b
ytop ; ybbottom ].

Let us assume that n and m are divisible by 2Q = 8. This assumption is not necessary but
helps simplify the notations. Let us write
yb = rowSC(A),

zb = colSC(A)

to denote labels obtained by applying the spectral clustering, respectively, on rows and columns
of the adjacency matrix A, the details of which are discussed in Section 5.2.2 below. We have
colSC(A) = rowSC(AT ). We also recall the LR classifier defined in (33). For matrices (or
vectors) A and B, we use [A; B] to denote column concatenation and [A B] to denote row
concatenation.
The general idea behind the partitioning scheme used in Algorithm 3, which is done by
sequential sampling without replacement, is to ensure that in each step where the LR classifier
17









A11
A21
A31
A41

A12
A22
A32
A42

A13
A23
A33
A43

A14
A24
A34
A44

(a) Gcol
1 (Step 3)

















A11
A21
A31
A41

A12
A22
A32
A42

A13
A23
A33
A43

A14
A24
A34
A44

(b) Grow
(Step 4)
1

















A11
A21
A31
A41

A12
A22
A32
A42

A13
A23
A33
A43

A14
A24
A34
A44









(c) G2 (Steps 6, 7)









A11
A21
A31
A41

A12
A22
A32
A42

A13
A23
A33
A43

A14
A24
A34
A44

(d) G3 (Step 9)









Figure 1: The four stages of partitioning in Algorithm 3. In each case, the collection of submatrices in
0
the partition is given a name which is used in the text. We have used to shorthand Aqq0 = A(q,q ) for simplicity. Block used in obtaining initial labels (a–b), in obtaining the first local parameter estimates (c),
and in the first application of LR classifier (d).

is applied, the initial labels used are independent of the sub-block of the adjacency matrix
under consideration. We do not require, however, that the initial labels be independent of the
estimates of the mean parameters Λ̂, since—as will be seen in Section 6.1—we have uniform
consistency of the LR classifier over all Λ̂ close to the truth. For example, in step 7, that is,
in the assignment ỹ (q) ← LR(A(q,q+2) , Λ̂(q+2) , z̃ (q+2) ), the claim is that z̃ (q+2) —at that stage
in the algorithm—is independent of A(q,q+2) but not necessarily of Λ̂(q+2) . This will become
clear in the following discussion where we keep track of the dependence of various estimates
through the algorithm. Note that in the description of Algorithm 3, we are using the computer
coding convention for in-place assignments, e.g., z̃ q gets updated in place and refers to different
objects at different points in the algorithm.
Figure 1 illustrates the partitions used in steps 2–9 of the algorithm. The collection of the
submatrices in the partition is given a name in each case. For example, Gcol
1 consists of the four
submatrices in Figure 1(a). Note that {Gcol
,
G
,
G
}
form
a
complete
partition
of the matrix
2
3
1
col
row
into disjoint blocks. Also, G1 and G1 involve the same elements of the matrix, i.e. they
cover the same portion of A. Thus, {Grow
1 , G2 , G3 } is also a complete cover of A with disjoint
row
blocks. Let us write G1 for the common portion of A covered by Gcol
1 and G1 .
row respectively, producing initial row and colSteps 3 and 4 operate on blocks in Gcol
1 and G1
umn labels. For example, in step 3, we apply row SC on each submatrix specified in Figure 1(a)
and obtain the label vectors (from the leftmost submatrix to the rightmost one):
[ỹ 0(1) ; ỹ (4) ], [ỹ (1) ; ỹ 0(2) ], [ỹ (2) ; ỹ 0(3) ], [ỹ (3) ; ỹ 0(4) ].

(34)

As a result of these steps, we obtain two sets of row labels ỹ = (ỹ (q) : q ∈ Z4 ) and ỹ 0 = (ỹ 0(q) :
q ∈ Z4 ), and similarly for the columns labels. Neither of ỹ or ỹ 0 is necessarily a consistent set
of labels for the whole matrix, since the cluster labels for individual pieces y (q) and ỹ 0(q) need
not match (e.g., cluster 1 in one piece could be labeled cluster 2 in another piece.). However,
if the sub-block labels (34) are sufficiently close to the truth, we can use the overlap among
them to find a global set of labels that are consistent with each block of ỹ and ỹ 0 . This is what
the Match operator in step 5 does, as will be detailed in Section 5.1. The resulting updated
global row and column labels only depend on G1 portion of A. Steps 6–13 go through the
following phases:
First local parameter estimates (step 6): Having obtained good initial (global) row
and column labels, in Step 5, we obtain estimates of the local mean parameters Λ̂(q+2) for
the submatrices in G2 as in Figure 1(c). Note for example, that Λ̂(q+2) computed in this step
depends on blocks A(q,q+2) and on G1 through z̃ (q+2) . Collectively, the estimates {Λ̂(q+2) : q ∈
18

Z4 } in Step 6 depend on G1 ∪ G2 portion of A.
First LR classifier (steps 7–8): Using the estimates of the (local) row mean parameters, in Step 6, we apply the LR classifier, ỹ (q) ← LR(A(q,q+2) , Λ̂(q+2) , z̃ (q+2) ) to each of the
submatrices in G2 (in Figure 1(c)). Here, Λ̂(q+2) depends on the same block A(q,q+2) on which
we apply LR classifier, but the dependence is not problematic due the uniform consistency of
LR classifier in parameters (Lemma 2). However, we note that z̃ (q+2) is a function of G1 blocks
of A, hence independent of A(q,q+2) which is key in our arguments. We will similarly apply
the LR classifier on the columns of G2 , and obtain z̃ (q) . By the end of step 8, the updated
labels ỹ and z̃ will depend on blocks in G1 ∪ G2 ; these labels will be much more accurate
(Mis ≈ exp(−I/Q)) than the initial labels obtained by spectral clustering.
Second parameter estimates (steps 9–10): Using the more accurate labels of step 8,
we obtain the local mean parameters Λ̂(q+3) in step 9 for the submatrices in G3 (Figure 1(d)).
This step is similar to step 6, but due to the much more accurate labels, the parameter estimates
are much more accurate
P mean parameter is the sum of local mean
P as well. Since the global
parameters, i.e. Λ = q∈[Q] Λ(q) , we use Λ̂ := q Λ̂(q) to estimate Λ in step 10. It is worth
recalling that the true local mean parameters, do not depend on the block row index; see (25).
Second LR classifier (step 11): Using the more accurate estimates of (global) row mean
parameters Λ̂ from step 10 and the more accurate labels z̃ in step 8, in step 11 we apply the
LR classifier ybtop ← LR(Atop , Λ̂, z̃) on Atop . We note that Atop in this step is independent of
z̃ (as well as Λ̂). This second LRC application is what brings us from very accurate labels
(Mis ≈ exp(−I/Q)) to almost optimal (Mis ≈ exp(−I)), as argued in Section 7.
Bottom half (steps 12–13): The same process is repeated in step 12, after swapping
the top and bottom halves of A, to get the bottom portion of the row labels. No matching is
required in step 13 when concatenating the top and bottom pieces to form a global set of row
labels yb. This is because the LR classifiers produce the same cluster labels; see Section 7.

5.1

Matching step

Let us describe the details of the matching step in Algorithm 3. Although, the idea is intuitively
clear, formally describing the procedure is fairly technical. In order to understand the idea,
consider the two-block labels ỹ (q−1,q) := [ỹ (q−1) ; ỹ 0(q) ], for q = 2, 3, that is,
ỹ (1,2) := [ỹ (1) ; ỹ 0(2) ],

ỹ (2,3) := [ỹ (2) ; ỹ 0(3) ].

We will detail how these two sets of labels can be fused together to generate a set of consistent
labels for the three-block true label vector y (1,2,3) := [y (1) ; y (2) ; y (3) ]. The two (overlapping)
two-blocks of the true label vector are also denoted as
y (1,2) := [y (1) ; y (2) ],

y (2,3) := [y (2) ; y (3) ].

More generally, we let y (q−1,q) = [y (q−1) ; y (q) ], similar to the notation for estimated blocks.
Recall our notation σ ∗ (· → ·) for (an) optimal permutation between two sets of labels (cf.
Section 2.3). Let us define

σq−1,q := σ ∗ ỹ (q−1,q) → y (q−1,q) , σq := σ ∗ (ỹ (q) → y (q) ), σq0 := σ ∗ (ỹ 0(q) → y (q) ).
(35)
19

ỹ (1)
ỹ

(1,2)

=
ỹ

0 (2)

ỹ (2)
ỹ (2,3) =
ỹ

0 (3)

σ1
σ12
σ20
σ2
σ23
σ30

ỹ

y (1)

0 (2)

σ20

= y (1,2)
y (2)

σ ∗ (ỹ 0(2) → ỹ (2) )

σ2−1

= y (2,3)
y (3)

y (2)

ỹ (2)

(a)

(b)

Figure 2: Pictorial depiction of the matching step. (a) Two-block and sub-block optimal permutations
to the truth. When ỹ (1,2) ≈ y (1,2) , we have σ1 = σ20 = σ1,2 and similarly ỹ (2,3) ≈ y (2,3) implies
σ2 = σ30 = σ2,3 . (b) Commutative diagram depicting how the missing permutation σ2−1 ◦ σ20 can be
obtained by matching observed labels ỹ 0(2) and ỹ (2) . See Section 5.1 for details.

Thus, for example we have
σ1,2 = σ ∗ (ỹ (1,2) → y (1,2) ),

σ2 = σ ∗ (ỹ (2) → y (2) ),

σ30 = σ ∗ (ỹ 0(3) → y (3) ),

and so on, as depicted in Figure 2(a). In other words, each of these permutations is the
optimal permutation from the corresponding block of the underlying estimated label to that of
the truth. Let us write ỹ (1,2) ≈ y (1,2) to mean that the two sets of labels are sufficiently close
(to be made precise later).
The first claim is that ỹ (1,2) ≈ y (1,2) implies that the underlying sub-blocks have the same
optimal permutation to the truth as the original two-block label, i.e.,
ỹ (1,2) ≈ y (1,2) =⇒ σ1 = σ20 = σ1,2
and similarly ỹ (2,3) ≈ y (2,3) =⇒ σ2 = σ30 = σ2,3 . The second claim is that each sub-block has
“almost” the same misclassification error as the bigger two-block. To see this, recall the direct
misclassification rate introduced in Section 2.3, i.e., misclassification rate without applying any
permutation (or equivalently with the identity permutation). We have


dMis σ2,3 (ỹ (2,3) ), y (2,3) = Mis ỹ (2,3) , y (2,3) ≤ ε.
(36)

where the inequality is by assumption (ε being the rate achieved by the spectral clustering
algorithm). A similar expression holds with (2, 3) replaced with (1, 2). Now (36) implies


dMis σ2 (ỹ (2) ), y (2) = dMis σ2,3 (ỹ (2) ), y (2) ≤ 2ε = ε0
(37)
where the equality uses σ2 = σ2,3 . To see the inequality, let n2,3 , n2 and n3 be the lengths of
y (2,3) , y (2) and y (3) . Then,



n3
n2
dMis σ2,3 (ỹ (2,3) ), y (2,3) =
dMis σ2,3 (ỹ (2) ), y (2) +
dMis σ2,3 (ỹ 0(3) ), y (3)
n2,3
n2,3

and the result follows since we have n2 = n3 = n2,3 /2 by construction. Note that dMis has the
property of being easily distributed over sub-blocks as opposed to Mis. Similarly to (37), we
20


obtain dMis σ30 (ỹ (3) ), y (3) ≤ ε0 considering the second component of ỹ (2,3) and y (2,3) . Applying

the same argument to indices (1, 2), we conclude similarly that dMis σ1 (ỹ (1) ), y (1) ≤ ε0 and

dMis σ20 (ỹ 0(2) ), y (2) ≤ ε0 .
Now consider the following three block vector undergoing transformation
 −1
 −1


 −1



σ2 ◦ σ1 (ỹ (1) )
σ2 ◦ σ20 (ỹ (1) )
σ1 (ỹ (1) )
σ2 ◦ σ1 (ỹ (1) )



 = 
 = 

ỹ (2)
ỹ (2)
 σ2 (ỹ (2) ) → σ2−1 ◦ σ2 (ỹ (2) ) → 
→
.
σ30 (ỹ (3) )

σ2−1 ◦ σ30 (ỹ (3) )

ỹ (3)

ỹ (3)

The leftmost vector has dMis of at most ε0 relative to y (1,2,3) by the previous arguments, and
since Mis ≤ dMis, we have the same bound on Mis rate for the leftmost vector. The first
transformation keeps the same Mis rate since we are applying a single permutation σ2−1 to all
elements. The second transformation is in fact an equality, using σ30 = σ2 established earlier.
The third transformation/equality follows similarly by σ1 = σ20 . Thus, if we can recover σ2−1 ◦σ20
from data, we can construct a consistent three-block label having Mis ≤ ε0 .
The third and final claim is that this is possible, and in fact we have
σ2−1 ◦ σ20 = σ ∗ (ỹ 0(2) → ỹ (2) )

(38)

that is, σ2−1 ◦ σ20 can be obtained (assuming ε0 is sufficiently small) by optimally matching ỹ 0(2)
to ỹ (2) , both of which we observe in practice. See the commutative diagram in Figure 2(b). In
order to make the above argument precise, we need to justify the first and third claims. We
will discuss the details in Section 6.4. The above matching process can be repeated over all the
two-blocks ỹ (q−1,q) to get a consistent set of global labels whose overall misclassification rate is
no more than twice that of the original two-blocks (cf. ε0 versus ε).

5.2
5.2.1

Results for Algorithm 3
General initialization

Before studying the spectral initialization, let us give a general bound on the misclassification
rate of Algorithm 3, assuming sufficiently good quality initial labels. In particular, assume that
the initial labels obtained in steps 3 and 4 of the algorithm are γ1 -good in the sense of (B3),
with γ1 satisfying
h 1  I
i 1
col
Imin
min
γ1 ≤
∧
∧ .
(39)
384β 2 ω 8LkΛk∞ 16KkΓk∞
4

Any other initialization algorithm besides spectral clustering can be used, as long as the above
guarantee on its output holds. We also need the following weaker version of (A4):
h I
ia 
col
min ∧ Imin
βω(kΛk∞ ∨ kΓk∞ ) = o
, for some a > 0.
(A40 )
Q log Q(K ∨ L)

col → ∞, Λ
Theorem 2. Assume that the model parameters satisfy Imin ∧ Imin
min → ∞, (A3)
0
and (A4 ), and the initial labels satisfy (39). Then, for some ζ = o(1), yb output by Algorithm 3
satisfies
 X
1



1 
Misk yb, y = O ω
1+
exp −Ikr −
− ζ log Λmin
(40)
εkr
2
r6=k

for every k ∈ [K] with probability 1 − o(1).

21

Algorithm 4 SC-RRE
Apply degree-reduction regularization A to obtain Are , as discussed in [ZA18].
(K∧L)
b1 Σ̂Z
bT , the (K ∧ L)-truncated SVD of Are .
2: Obtain Are
=Z
2
b1 Σ̂) where K is an isometry-invariant κ-approximate kmeans algorithm.
3: Output K (Z

1:

We refer to Section 3 for the definition of the parameters involved in the rate given in (40).
5.2.2

Spectral initialization

Theorem 2 requires initial labels that satisfy (39). A spectral algorithm, namely SC-RRE
given in Algorithm 4 can provide such initialization. (The acronym stands for reduced rank
efficient spectral clustering.) The algorithm is presented and analyzed in [ZA18]. One performs
a truncated SVD of rank r := K ∧ L on a regularized version of the adjacency matrix Are to
obtain Zb1 Σ̂Zb2T , where Σ̂ is the diagonal matrix retaining the top r largest singular values of Are ,
b1 ∈ Rn×r and Z
b2 ∈ Rn×r are the matrices of the corresponding singular vectors. One then
and Z
b1 Σ̂, rather than Zb1 which is the more common approach
runs a k-means type algorithm on Z
in spectral clustering. This allows one to state consistency results without a reference to the
gap in the spectrum of Are , while still retaining the attractive feature of the latter approach,
namely, the computational efficiency of running k-means on a matrix of reduced dimension.
We refer to [ZA18] for a detailed discussion. In particular, one has the following bound on the
misclassification rate of SC-RRE:
Theorem 3 ([ZA18]). Let α = m/n and Λ2∧ := mint6=s kΛs∗ − Λt∗ k2 . Consider the spectral
algorithm given in Algorithm 4, assume that for a sufficiently small C1 > 0,
β 2 KL(K ∧ L) α

kΛk∞
≤ C1 (1 + κ)−2 .
Λ2∧

(41)

Then the algorithm outputs estimated row labels ỹ satisfying
Mis(ỹ, y) ≤ C1−1 (1 + κ)2 βL(K ∧ L) α

 kΛk 
∞
2
Λ∧

.

Here, Λs∗ refers to the sth row the mean parameter matrix Λ (cf. Section 2.1) Combining
Theorems 2 and 3, one obtains Theorem 1. Some work is required to translate the bound of
the Theorem 3 to be applicable to sub-blocks. The details of are left to Section 7.2.

6

Preliminary analysis

We start by analyzing the properties of the operators introduced in Sections 4.1 and 4.2, for
some fixed (deterministic) initial labels ỹ and z̃. We assume that these labels satisfy:
Mis(ỹ, y) ≤

γ
,
βK

Mis(z̃, z) ≤

γ
.
βL

(B3)

We call such labels γ-good. Throughout, Λ̃ will be used to denote a generic deterministic
approximation of the true row mean parameter Λ. The relative `∞ ball of radius δ centered at
22

Λ, that is,
BΛ (δ) := {Λ̃ : kΛ̃ − Λk∞ ≤ δkΛk∞ },

(42)

will play a key role in our arguments. For sufficiently small δ and true Λ, BΛ (δ) will be the
set of δ-good row mean parameters.

6.1

Fixed label analysis

We first present the analysis assuming that all the operations are performed on the entire
adjacency matrix A. In Section 6.2, these results are extended to be applicable to sub-blocks
of A. Recall the definitions of the mean parameters an their estimates from Section 2.1. In
particular, we recall that λk∗ (y, z̃) is the mean of bi∗ (z̃) for any node i with yi = k. These
mean parameters form the kth row of Λ(y, z̃). Our first main lemma illustrate that whenever
the initial labels z̃ and ỹ are γ-good, then the parameters Λ(y, z̃) as well as the corresponding
estimates Λ̂ defined in (21) are close to the truth,Λ.
Lemma 1 (Parameter consistency).
Let Cγ = Cγ,β = β 2 γ/(1 − γ), assume that 6Cγ ω ≤ 1,

2c
3
and let hc (τ ) := 4c τ log 1 + 3 τ . Then under assumptions (A1), (A2) and (B3), we have
(a) kΛ(y, z̃) − Λk∞ ≤ Cγ kΛk∞ ,

kΛ(y, z̃)k∞ ≤ 2kΛk∞ .

(b) kΛ(ỹ, z̃) − Λ(y, z̃)k∞ ≤ 2γkΛk∞ ,

kΛ(ỹ, z̃)k∞ ≤ 4kΛk∞ .

(c) kΛ̂ − Λ(ỹ, z̃)k∞ ≤ 4τ kΛk∞ with probability at least 1 − 2p1 where
 nΛ

min h1 (τ )
p1 = p1 (τ ; n, Λmin , β) := KL exp −
,
4βK

∀τ > 0,

(43)

and Λ̂ is as defined in (21). In particular, all the estimates Λ(y, z̃), Λ(ỹ, z̃) and Λ̂ are within
relative `∞ distance of at most 4(Cγ + τ ) from Λ.
The lemma is proved in Section 9.1. Note that the lemma implies that Λ̂ ∈ BΛ (4(Cγ + τ ))
with the stated probability.
Our second key lemma shows that the LR classifiers in (33) are uniformly dominated, over
Λ̃ ∈ BΛ (δ), by a single (perturbed) classifier. To state this result, recall the block compression
b(z̃) := B(A; z̃) given in (26), and define the following:
Yikr (bi∗ , Λ̃) := Ψ(bi∗ ; λ̃r∗ | λ̃k∗ ) =

L
X
`=1

bi` log

λ̃r`
+ λ̃k` − λ̃r` ,
λ̃k`

Zik (bi∗ , Λ̃) := 1{Yikr (Λ̃) ≥ 0, for some r 6= k}.
X
1
Sk (b, Λ̃) :=
Zik (bi∗ , Λ̃),
nk (y)

(44)
(45)
(46)

i:yi =k

where Ψ is the Poisson log-likelihood ratio defined in (32). Thus, Yikr is the (pseudo) loglikelihood ratio, for k, r ∈ [K], measuring the relative likelihood of row i having label k.

23

We note that Yikr (Λ̃) < 0, ∀r 6= k implies ybi := (LR(A, Λ̃, z̃))i = k. Thus, Sk (bi∗ , Λ̃) is the
misclassification rate for the LR classifier over the kth row-class, i.e., Misk (b
y , y). Let
Jkr = LkΛk∞ /Ikr

(47)

and recalling definitions of εkr , ω and β from Section 3, set
η 0 := η 0 (δ; Λ) = 8ωδLkΛk∞ = 8ωδJkr Ikr ,

(48)

ηkr := ηkr (δ; ω, β, m, Λ)
= 21δωLkΛk∞ +

h

i 1
1
5βL2 kΛk2∞
+ log 11ω
+ 1 − log Λmin .
m
εkr − 2ω(1 + εkr )δ
2

(49)

We have the following key lemma:

Lemma 2 (Uniformity of LRC in mean parameters). Fix any row label z̃ and let b = b(z̃) be
the corresponding column compression. Let Λ0 = Λ(y, z̃) be the row mean parameter associated
with b. Assume (A1), (A2), and Λ0 ∈ BΛ (δ) with 3 ωδ < 1. Then, for all k, r ∈ [K], k 6= r,
and all i : yi = k, we have the following bounds:
(a) With η 0 defined as in (48),
P ∃Λ̃ ∈ BΛ (δ), Yikr (bi∗ , Λ̃) ≥ 0



≤ exp(−Ikr + η 0 ).

(b) If in addition εkr − 2ωδ > 0, then with ηkr defined as in (49),


P ∃Λ̃ ∈ BΛ (δ), Yikr (bi∗ , Λ̃) ≥ 0 ≤ exp −Ikr + ηkr .

(50)

(51)

The proof of Lemma 2(b) appears in Section 9.5, and that of part (a) in Appendix A.5.

Remark 4 (Typical setting). In the error exponent in Lemma 2(b), i.e. −Ikr + ηkr , the first
three terms in (49) are positive and constitute the undesirable part of the bound. Our goal is
to keep these terms dominated at the final stage of the algorithm, i.e., make them o(log Λmin ),
by making δ sufficiently small. For now, let us introduce a simple typical setting to give some
idea of the order of ηkr . In the first reading, one can consider the case where β, ω = O(1),
Ikr  I → ∞ for all k, r and some I, and assume that LkΛk∞ /I = O(1) and (A5) holds. In
this setting, Jkr = O(1) and we have ηkr = C(δ + m−1 I)I − 21 log Λmin for some constant C.
Keeping these typical orders in mind will be helpful in understanding the statements of the
subsequent results.
It is also worth noting that we always have Jkr ≥ 21 . which follows from the general bound
Ikr ≤ 2LkΛk∞ . Another important quantity is Cγ in Lemma 1, which in the typical setting
behaves as Cγ  γ when γ → 0.
Combining Lemma 2 with the Markov inequality, we can get uniform control on the misclassification rate of the LR classifier in its parameter argument (i.e., Λ̂):
Lemma 3. Fix k ∈ [K] and z̃ ∈ [L]m . Let Λ̂ ∈ RK×L
be any random matrix and set yb(z̃) :=
+
LR(A, Λ̂, z̃). Assume that (50) holds. Then, for any u ∈ R, we have
 X

Misk yb(z̃), y ≤
exp −Ikr + η 0 + u ,
r6=k


with probability at least 1 − e−u − P Λ̂ ∈
/ BΛ (δ) . The result is also true if we replace η 0 by ηrk
when (51) holds.
24

Remark 5. Edge splitting (ES) was proposed in [AS15] to generate nearly independent copies
from a single network. One might ask whether combining the edge splitting idea with Lemma 3
is enough to give us a result similar to Theorem 1. In ES, edges are randomly assigned to two
graphs G1 and G2 , with probabilities q and 1 − q. The new graphs G1 and G2 will follow a
SBM with a reduced connectivity matrix (by a factor of q and 1 − q respectively). Hence, the
corresponding parameters Λ and I are reduced by the same factor; for example I will be scaled
to qI for G1 . Let us consider the typical setting where β, K, L, ω, εkr = O(1) and Ikr  I for
all k, r and some I; assume the connectivity matrix is symmetric, i.e., Λ = Γ and I = I col . Let
z̃ and ỹ be the labels obtained by performing biclustering on G1 . Lemma 3 in the best case
scenario, with the most favorable version of ηkr —i.e., ignoring the first three positive terms
in (49)—gives a misclassification rate


X


1
max{Mis ỹ, y , Mis z̃, z } ≤ γ2 :=
exp −qIkr − log(qΛmin ) + v
2
r6=k

for some v → ∞, w.h.p.. In the second stage, given the labels z̃ and ỹ, we obtain an estimate
of the (row) mean parameters based on G2 , using the natural estimator Λ̂2 = L (G2 , ỹ, z̃). We
then obtain the second stage labels yb := LR(G2 , Λ̂, z̃). Let Λ2 = (1 − q)Λ be the row mean
parameter of G2 . By Lemma 1, Λ̂2 ∈ BΛ2 (δ) w.h.p for some δ ≥ γ2 . By Lemma 3, and the
perturbation of information (Lemma 7) we have


X

1
Mis yb, y ≤ γ3 :=
exp −(1 − q)Ikr + C(1 − q)δkΛk∞ − log Λmin + u
2
r6=k

for some u → ∞ w.h.p.. To obtain result (14) in Corollary 1, we at least hope to have
qIkr + C(1 − q)γ2 kΛk∞ = o(log Λmin ).
So we need qIkr = o(log Λmin ) and (1 − q)γ2 kΛk∞ = o(log Λmin ). Assume that we have
qIkr = o(log Λmin ). Then,
γ2 =

X
r6=k



1
−1/2−o(1) √
exp −qIkr − log(qΛmin ) + v = O(Λmin
/ q).
2

However, this is not sufficient to show (1 − q)γ2 kΛk∞ = o(log Λmin ). Therefore, applying edge
splitting and Lemma 3 does not lead to the main result of this paper.

6.2

Analysis on subblocks

We now extend the analysis of Section 6.1 to be applicable to the sub-blocks obtained by random
partitioning. Some care needs to be taken since the true (row and column) mean parameters
of the sub-blocks are changed by partitioning, due to the change in the distributions of the
labels within each sub-block among the K × L classes. The deviations of the sub-block class
proportions from the global version will be controlled by a slack parameter ξ which will be set
at the final stage of the proof (see Section 7.1.2). Throughout this section, assumptions (A1)
and (A2) will be implicit in all the stated lemmas. We will also state the result for a general
2Q × Q partitioning scheme, although Q = 4 is enough for the analysis of Algorithm 3.
25

Recall that the class priors π` (z) for the full labels are defined in (29). We will use the
same notation for sublabels z (q) , that is, π` (z (q) ) is the proportion of labels in z (q) that lie in
class `. Note that we have
π` (z) =

n` (z)
,
m

π` (z (q) ) =

n` (z (q) )
,
m/Q

hence,

π` (z (q) )
n` (z (q) )
=Q
,
π` (z)
n` (z)

(52)

since z (q) has length m/Q. We similarly we have πk (y (q) ) = nk (y (q) )/(n/(2Q)). We will work
under the assumption that the partitioning scheme satisfies:
max |πk (y (q) ) − πk (y)| ≤ ξ
k,q

max |π` (z (q) ) − π` (z)| ≤ ξ,
`,q
 1
1 
ξ ≤ min
,
.
2βK 2βL
and

(B4a)
(B4b)

When these conditions hold, we call the scheme a good partition. We note that these conditions
combined with (A2) give,
1
1 1
3β
π` (z (q) )
− 1 ≤ ξLβ ≤
=⇒
≤ π` (z (q) ) ≤
π` (z)
2
2 βL
2L

(53)

and similarly for y (q) . It follows that both z (q) and y (q) satisfy (A1) with β replaced with 2β.
Each count nk (y (q) ) follows a hypergeometric distribution with parameters (n, nk (y), n/(2Q)),
that is, the number of nodes labeled k, in a sample of size n/(2Q), from a population of size
n, with a total of nk (y) nodes labeled k. The concentration of the hypergoemtric distribution
gives the following:
Lemma 4. (B4a) holds for random partitioning, with probability at least 1 − p2 , where

p2 = 2Q(K + L) exp − min(n, m)ξ 2 /Q .

(54)

The proof of this lemma and others in this section appear in Appendix A.1.

(q)

Lemma 5. Under (B4a) and (B4b), the true local mean parameters Λ(q) = (λk` ) satisfy:
(q)

λk` −
(q)

In particular, Λmin ≥

1
2Q Λmin ,

λk`
λk`
1 λk`
≤ (ξLβ)
≤
,
Q
Q
2 Q

kΛ(q) k∞ ≤

3
2Q kΛk∞

∀q, k, `.

(55)

and Λ(q) ∈ BΛ/Q (ξLβ) for all q ∈ [Q].

Our main lemma for the sub-blocks establishes the consistency of the local mean parameter
0
estimates Λ̂(q ,q) for a good partitioning scheme. This lemma is an extension of Lemma 1. We
recall the operator L from (20):

3
Lemma 6 (Local parameter consistency). Let Cγ = β 2 γ/(1−γ) and hc (τ ) := 4c
τ log 1 + 2c
3τ
as in Lemma 1 and assume that 72 Cγ ω ≤ 1. Fix the underlying partition and fix q, q 0 ∈ [Q],
and labels z̃ and ỹ. Let
0
0
0
Λ̂(q ,q) = L (A(q ,q) , ỹ (q ) , z̃ (q) ).

26

Assume that the partition satisfies (B4a) and (B4b), and the pairs (z̃ (q) , z (q) ) and (ỹ (q) , y (q) )
satisfy the misclassification rate in (B3). Then,

0
kΛ̂(q ,q) − Λ(q) k∞ ≤ 24 Cγ + 6τ kΛ/Qk∞ , and

0
kΛ̂(q ,q) − Λ/Qk∞ ≤ 24 Cγ + 6τ + ξLβ kΛ/Qk∞
with probability at least 1 − 2p3 , where

 nΛ

min h1 (τ )
.
p3 = p3 (τ ; n, K, Λmin , Q) := KL exp −
32Q2 βK

(56)

We also have
0

(a) kΛ(q ,q) (y, z̃) − Λ(q) k∞ ≤ 4Cγ kΛ(q) k∞ .
0

0

(b) kΛ(q ,q) (ỹ, z̃) − Λ(q ,q) (y, z̃)k∞ ≤ 2γkΛ(q) k∞ .
0

0

(c) kΛ̂(q ,q) − Λ(q ,q) (ỹ, z̃)k∞ ≤ 4τ kΛ(q) k∞ , with probability at least 1 − 2p3 .
Remark 6. Similar results to those obtained above hold for the column parameters. Recall
that the dual to the row mean parameters Λ are the column mean parameters Γ. The result
of Lemma 5 can be translated to the column version by making the following substitutions
Λ → Γ, Q → 2Q and L ↔ K. For Lemma 6, in addition we need to make n → 4m. (The
reason for this is that in (115), in the proof, we need to replace n/2Q with m/Q, and Λmin /2Q
with Γmin /(4Q), and the combination of the aforementioned substitutions achieves this. We
also note for future reference that the corresponding ω inflation by a factor of 3 remains true
for column parameters.) After these substitutions, we obtain the same constant in (56), that
is, p3 has to be replaced with
p03 := p3 (τ ; 4m, L, Γmin , 2Q) = p3 (τ ; m, L, Γmin , Q).

6.3

(57)

Perturbation of information

Recall the definition of Chernoff information from (8), and let us write Ikr = Ikr (Λ) to explicitly show its dependence on the mean parameter matrix Λ. The following lemma, proved in
Appendix A.1, bounds the perturbations of Ikr (Λ) in Λ:
Lemma 7. Under (A1), for any Λ̃ ∈ BΛ (δ), we have |Ikr (Λ̃) − Ikr (Λ)| ≤ 2ωδLkΛk∞ .

6.4

Analysis of the matching step

In this section, we fill in the details of the argument sketched in Section 5.1. Specifically, we
need to give sufficient conditions so that the first and the third claims of Section 5.1 hold. We
will use the following two lemmas. Recall the notation σ ∗ (ỹ → y) introduced in Section 2.3 to
denote the optimal permutation from the set of labels ỹ to another set y.
Lemma 8. Let ỹ, y ∈ [K]n , and assume that dMis(ỹ, y) <

1
2

mink πk (y). Then,

(a) σ ∗ (ỹ → y) = id, the identity permutation, and this optimal permutation is unique, and
(b) πk (ỹ) > 21 πk (y) for all k.
27

Note that Lemma 8 implies that if dMis(σ(ỹ), y) <
then σ ∗ (ỹ → y) = σ.

1
2

mink πk (y) for some permutation σ,

Lemma 9. Consider three sets of labels y, ỹ, ỹ 0 ∈ [K]n , and assume that
max{ Mis(ỹ, y), Mis(ỹ 0 , y) } <

1
min πk (ỹ).
4 k

Let σ = σ ∗ (ỹ → y) and σ 0 = σ ∗ (ỹ 0 → y). Then, σ −1 ◦ σ 0 = σ ∗ (ỹ 0 → ỹ).
The first claim follows from Lemma 8, under the further assumption:
Mis(ỹ (q−1,q) , y (q−1,q) ) <

1
,
32βK

q ∈ [Q].

(58)

Using the permutation notations (35) of Section 5.1, we have:
Corollary 2. Under assumptions (A2), (B4a), (B4b) and (58), σq−1,q = σq−1 for all q ∈ [Q].
The third and final claim of Section 5.1 follows from Lemmas 8 and 9, by applying them
to the sub-block labels y (2) , ỹ (2) , ỹ 0(2) :
Corollary 3. Under assumptions (A2), (B4a), (B4b) and (58), σq−1 ◦ σq0 = σ ∗ (ỹ 0(q) → ỹ (q) )
for all q ∈ [Q].
The proofs of the results of this section are deferred to Appendix A.1.

7

Analysis of Algorithm 3

7.1

Proof of Theorem 2

We start with the high-level analysis of Algorithm 3 in Section 7.1.1. This analysis is parametrized
by many parameters such as ξ , τ1 , τ1col , τ2 , etc. This allows us to give the high-level idea of the
mechanics of the proof without making the arguments obscured by the expressions ultimately
chosen for these parameters. In Section 7.1.2, we make specific choices about these parameters
and finish the proof of Theorem 2.
7.1.1

Parametrized analysis of Algorithm 3

We now have all the pieces for analyzing Algorithm 3. Let ỹstep 5 and z̃step 5 be the labels from
step 5 of of Algorithm 3. As before, in all the lemmas stated, (A1) and (A2) will be implicitly
assumed. Consider the following event:
n
o
(q)
(q)
Aγ := ỹstep 5 and z̃step 5 satisfy (B3) with parameter γ, for all q ∈ [Q] .
We implicitly assume that clusters in z̃step 5 and ỹstep 5 are relabeled according to optimal
permutation relative to the truth. In other words, z̃step 5 and ỹstep 5 in the above event are not
the raw output of the algorithm, but the relabeled versions (which we do not have access to in
practice, but are well-defined and can be used in the proof.) When γ is sufficiently small, this
implies that community k in z̃step 5 is the same as community k in z̃, for all k ∈ [Q].
28

Let Π be the random partition used in Algorithm 3, and let P be the event that Π satisfies
condition (B4a). By Lemma 4, we have P(P) ≥ 1 − p2 where p2 is given in (54). For the
most part, we will work on events of the form Aγ1 ∩ P. Let us also establish some terminology.
By the probability “on an event P”, we mean the probability under the restricted measure
PP := P(· ∩ P). For example, if D = {property X holds}, we will say that “property X fails”
on P with probability at most q if P(Dc ∩ P) ≤ q. In this case, if P holds with high probability,
say ≥ 1 − p2 , and q is small, then D holds with high probability as well: P(D) ≥ 1 − q − p2 .
(q)

Let Λ̂step 6 = L (A(q−2,q) , ỹ (q−2) , z̃ (q) ), q ∈ Z/QZ, be the first local parameter estimates
obtained in step 6 of Algorithm 3 (it is easier to work with the shifted index), and let
δ1 := 24 Cγ1 + 6τ1 + ξLβ.

(59)

A better name for δ1 , and τ1 would be δ1row , and similarly τ1row contrasting with δ1col and τ1col
defined later in (62). However, for simplicity, we drop the “row” qualifier here. Recall that
ξ is a parameter controlling the tail probability related to the random partition, while τ1 will
be controlling the tail probability p3 (τ1 ) related to the local parameter estimates in Lemma 6.
These parameters will be optimized at the end of the argument (see Section 7.1.2).
Lemma 10 (First local parameters). Assume (B4b) and 72 Cγ1 ω ≤ 1, and let δ1 be as defined
in (59). Then, on event Aγ1 ∩ P,
(q)

Λ̂step 6 ∈ BΛ(q) (δ1 ),

∀q ∈ ZQ ,

fails with probability at most 2Q p3 , where p3 = p3 (τ1 ) as given in (56).
Proof. Conditioning on blocks G1 (cf. Section 5) of the (bottom) adjacency matrix Abottom —
(G1 )
(q)
denoted as Abottom
—the distribution of blocks A(q−2,q) , q ∈ ZQ used in defining Λ̂step 6 is
not changed. Under this conditioning, both initial labels ỹstep 5 and z̃step 5 are deterministic,
(q)
hence the results of Section 6.2 apply. We will apply Lemma 6 to Λ̂step 6 . Let us verify the
(q)

(q)

conditions of the lemma. On Aγ1 , for all q ∈ [Q], the sublabel pairs (z̃step 5 , z (q) ) and (ỹstep 5 , y (q) )
satisfy (B3). On P, condition (B4a) holds for the random partition and (B4b) holds by
assumption. Recall that the random partition is independent of all else, hence conditioning
on it does not change the distribution of blocks A(q−2,q) , q ∈ ZQ either. We may then apply
(G1 )
Lemma 6 to conclude that for every q ∈ ZQ , conditioned on the partition Π and Abottom
, the
(q)
event {Λ̂step 6 ∈
/ BΛ(q) (δ1 )} ∩ Aγ1 ∩ P holds with probability ≤ 2p3 . Let us write
 (q)
D = Λ̂step 6 ∈ BΛ(q) (δ1 ), ∀q ∈ ZQ

which is the desired event in this lemma. Using the union bound, and removing the conditioning, we have P(Dc ∩ Aγ1 ∩ P) ≤ 2Qp3 , unconditionally. The proof is complete.
Next we consider the first LR classifier application. Let ỹstep 7 be the row label estimates
in step 7. That is, we have
(q−2)

(q)

(q)

ỹstep 7 = LR A(q−2,q) , Λ̂step 6 , z̃step 5
for which we have the following bound on misclassification rate:
29



Lemma 11 (First LR classifier). Under the assumptions of Lemma 10, further assume that
9ωδ1 < 1. Let η step 7 := 2η 0 (δ1 ; Λ/Q) where η 0 (·) is defined in (49). Then, on event Aγ1 ∩ P,
(q)

Misk ỹstep 7 , y (q)



≤

X
r6=k

 I

kr
row
exp −
+ η step 7 + u =: γ2k
,
Q

∀q ∈ ZQ ,

(60)

fails with probability at most Q(e−u + 2Q p3 ) where p3 = p3 (τ1 ) as given in (56).
Proof. Fix q ∈ ZQ and consider ỹ (q−2) . As in the proof Lemma 10, we condition on blocks in
(q)
G1 so that z̃step 5 can be assumed deterministic. We will apply Lemma 2(a) to the subblock
A(q−2,q) . As discussed earlier, the corresponding ω is inflated to 3ω, hence we need 3(3ω)δ1 < 1
which we have assumed. We also note that Λ(q−2,q) (y, z̃) and Λ(q) play the role of Λ(y, z̃) and
Λ in Lemma 2(b), and we have the needed condition Λ(q−2,q) (y, z̃) ∈ BΛ(q) (δ1 ) from Lemma 6.
(q−2,q)
(q)
Let bi∗
be the row block compression of A(q−2,q) based on z̃step 5 . Then, Lemma 2(a) gives
P

n
o



(q−2,q)
(G1 )
(q)
∃ Λ̃ ∈ BΛ(q) (δ1 ), Yikr bi∗
, Λ̃ ≥ 0 ∩ Aγ1 ∩ P Abottom
, Π ≤ exp −Ikr + η (q)

(61)

(q)

for all rows i (in row block q − 2) with yi = k. Here Λmin is the minimum element of Λ(q) , and
(q)

Ikr := Ikr (Λ(q) )

≥

Ikr
− 2ωδ1 LkΛ/Qk∞
Q

η (q) := η 0 (δ1 ; Λ(q) ) ≤ (1 + δ1 ) η 0 (δ1 ; Λ/Q)
where η 0 (δ1 ; Λ(q) ) = 8ωδ1 LkΛ(q) k∞ as defined in (49). The first inequality uses Lemma 7 and
the second is obtained using the definition of η 0 (·) combined with Λ(q) ∈ BΛ/Q (δ1 ) (Lemma 5)
which implies kΛ(q) k∞ ≤ (1 + δ1 )kΛ/Qk∞ . By taking expectation in (61), the same bound
holds unconditionally.
(q)
By Lemma 10, on event Aγ1 ∩P, we have Λ̂step 6 ∈
/ BΛ(q) (δ1 ) with probability at most 2Q p3 .
Then, applying Lemma 3, we conclude that
X


(q−2)
(q)
Misk ỹstep 7 , y (q−2) ≤
exp −Ikr + η (q) + u
r6=k

fails on Aγ1 ∩ P with probability ≤ e−u + 2Q p3 , for each q ∈ [Q]. Note that
(q)

−Ikr + η (q) ≤ −

Ikr
+ (1 + δ1 + 9−1 ) η 0 (δ1 ; Λ/Q).
Q

Since 9ωδ1 < 1 implies δ1 < 9−1 (recall ω ≥ 1), we have 1 + δ1 + 9−1 < 2. Combining with the
previous bound and applying the union bound over q gives the result.
Note that we have called the rate in (60) γ2row for the (column) misclassification rate based
on the row information. This rate is faster than initial rate γ1 . Repeating the procedure in
steps 6 and 7 for the column labels—as prescribed in step 8 in Algorithm 3–we obtain a similar
(q)
rate for the misclassification rate of z̃step 8 relative to z (q) which we call γ2col . In deriving γ2col ,
we have to make the substitutions in Remark 6, and particular, Λ → Γ where Γ is the column
30

mean parameters defined in Section 2.1. (A minor exception is when counting the number of
blocks which will still be Q rather than 2Q.) Recall the definition of the column information
col ) from (9). Letting
matrix (I`r
δ1col := 24 Cγ1 + 6τ1col + ξKβ,

(62)

we obtain the following counterpart of Lemma 11:
Corollary 4 (First LR classifier, column version). Under the assumptions of Lemma 10, further assume that 9ωδ1col < 1. Let η step 8 := 2η 0 (δ1col ; Γ/(2Q)) where η 0 (·) is defined in (49).
Then, on event Aγ1 ∩ P,
(q)

Mis` z̃step 8 , z (q)



≤

X
r6=`

 I col

col
exp − `r + η step 8 + ucol =: γ2`
,
2Q

col

fails with probability at most Q(e−u

∀q ∈ ZQ ,

(63)

+ 2Q p03 ), where p03 = p03 (τ1col ) as given in (57).

row
col , γ row := max
Let γ2col := maxk∈[K] γ2k
`∈[L] γ2` and
2

γ2 := max{βKγ2row , βLγ2col }.

(64)

By (7), we have that (60) and (63) imply

γ2
(q)
Mis ỹstep 7 , y (q) ≤ γ2col ≤
,
βK


γ2
(q)
Mis z̃step 8 , z (q) ≤ γ2row ≤
βL

(65)

Thus, if we consider the following event:
n
o
(q)
(q)
Bγ := ỹstep 7 and z̃step 8 satisfy (B3) with parameter γ, for all q ∈ [Q] ,

after Step 8, we can work on Bγ2 ∩ P which holds with high probability: Combining Lemma 11
and Corollary 4, by union bound, P(Bcγ2 ∩ Aγ1 ∩ P) ≤ Q(2e−u + 2Q(p3 + p03 )), hence

P(Bγ2 ∩ P) ≥ P Bγ2 ∩ Aγ1 ∩ P


= P Aγ1 ∩ P − P Bcγ2 ∩ Aγ1 ∩ P
≥ 1−

P(Acγ1 )

c

− P(P ) − Q(2e

−u

+ Q(p3 +

(66)
p03 )).

(q)

Let Λ̂step 9 = L (A(q−3,q) , ỹ (q−3) , z̃ (q) ), q ∈ ZQ , be the second local parameter estimates
obtained in step 9 of Algorithm 3. Let
δ2 := 24Cγ2 + 6τ2 .

(67)

Lemma 12 (Second local parameters). Assume (B4b) and 72 Cγ2 ω ≤ 1, and let δ2 be as
defined in (67). Then, on event Bγ2 ∩ P,
(q)

Λ̂step 9 ∈ BΛ(q) (δ2 ),

∀q ∈ ZQ

fails with probability at most 2Q p3 , where p3 is given in (56).
31

Proof. Conditioning on blocks G1 ∪ G2 (cf. Section 5) of the adjacency matrix A, the distri(q)
bution of blocks A(q−3,q) used in defining Λ̂step 9 is not changed. Under this conditioning, both
initial labels ỹstep 7 and z̃step 8 are deterministic, hence the results of Section 6.2 apply. On Bγ2 ,
(q)
(q)
for all q ∈ ZQ , the sublabel pairs (z̃step 8 , z (q) ) and (ỹstep 7 , y (q) ) satisfy (B3). The rest of the
proof follows that of Lemma 10.
The key is that δ2 is much smaller than δ1 , due to γ2  γ1 (typically), i.e., the second
P (q)
parameter estimates are much more accurate. Let Λ̂step 10 =
q Λ̂step 9 be the estimate of
the global mean parameters obtained in step 10 of Algorithm 3. According to Lemma 12,
on Bγ2 ∩ P,
(q)

kΛ̂step 9 − Λ/Qk∞ ≤ δ2 kΛ/Qk∞ , ∀q

hence,

kΛ̂step 10 − Λk∞ ≤ δ2 kΛk∞

(68)

fails with probability ≤ 2Q p3 , where we have used triangle inequality. That is, on Bγ2 ∩ P,
we have Λ̂step 10 ∈ BΛ (δ2 ) with high probability.
(q)

Remark 7. Note that we could have used QΛ̂step 9 (for any q ∈ ZQ ) as our estimate Λ̂step 10 ,
leading to the same bound as in (68). The results would be the same, though in practice, we
expect the version given in the Algorithm 3 to perform better. We also note that on Bγ2 ,
(q)
the sublabels (z̃step 8 , q ∈ ZQ ) automatically define a consistent global label vector z̃step 8 , and
similarly for row labels ỹstep 7 .
Lemma 13 (Second LR classifier). Under the assumptions of Lemma 12, further assume
that

step 11
:= ηkr δ2 ; ω, β, m, Λ . Then,
δ2 defined in (67) satisfies 3ωδ2 < 1 and 6Cγ2 ω ≤ 1. Let ηkr
on event Bγ2 ∩ P,


X

step 11
Misk ybtop , ytop ≤
+ v =: γ3 ,
(69)
exp −Ikr + ηrk
r6=k

fails with probability at most e−v + 2Q p3 where p3 = p3 (τ2 ) as given in (56). The same result
step 11
holds for ηkr
= η 0 (δ2 ; Λ).

Proof. As in the proof Lemma 12, we condition on blocks in G1 ∪ G2 so that z̃step 8 can be
assumed deterministic. We will apply Lemma 2(b) to Atop . Let Top ⊂ [n] denote the row
indices of Atop . Since all the columns are present in Atop , we can directly apply Lemma 2(b) (in
contrast to the argument in Lemma 11), that is, the relevant row mean parameters are Λ(y, z̃)
and Λ—the same as those for the whole matrix A. The needed condition Λ(y, z̃) ∈ BΛ (δ2 )
11
is supplied by Lemma 1. Let bstep
= bi∗ (z̃step 8 ) be the block compression in step 11 of the
i∗
(G1 ∪G2 )
algorithm. Then, Lemma 2(b) gives (after conditioning on Abottom
and then removing the
conditioning as in (61))
n
o




11
step 11
P ∃ Λ̃ ∈ BΛ (δ1 ), Yikr bstep
,
Λ̃
≥
0
∩
B
≤
exp
−I
+
η
.
(70)
γ
kr
2
i∗
kr

for any i ∈ Top with yi = k. By (68), on Bγ2 ∩ P, we have Λ̂step 10 ∈
/ BΛ (δ2 ) with probability
at most 2Q p3 . Then, applying Lemma 3, we conclude (69) as desired. The last statement of
the theorem follows if we apply Lemma 2(a) in place of Lemma 2(b) throughout.

32

The same exact bound holds for ybbottom in step 12, with the same probability. Hence,
by union bound, the same bound on misclassification rate holds for the final row labels yb in
step 13, with probability inflated by a factor of 2; that is, Misk (b
y , y) ≤ γ3 fails on Bγ2 ∩ P,
with probability at most 2(e−v + 2Q p3 ).
To summarize, under the conditions of the lemmas, we have




P Misk (b
y , y) > γ3 ≤ P Misk (b
y , y) > γ3 ∩ Bγ2 ∩ P + P (Bγ2 ∩ P)c


≤ 2 e−v + 2Q p3 (τ2 ) + P (Bγ2 ∩ P)c


col
≤ 2 e−v + 2Q p3 (τ2 ) + P(Acγ1 ) + P(Pc ) + Q 2e−u∧u + Q(p3 (τ1 ) + p03 (τ1col )
(71)
where γ3 is the rate given in (69) and the second inequality uses (66).
7.1.2

Choosing the parameters

It remains to choose the parameters, τ1 , τ2 , ξ, etc. to simultaneously achieve the desired rate
for γ3 and ensure that the probability in (71) is o(1).
Proof of Theorem 2. First row LR classifier. Let us write τ1row = τ1 for clarity. Under
our assumptions, we will have γ2 ≤ γ1 ≤ 1/2 so that Cγi ≤ 2β 2 γi for i = 1, 2, recalling the
definition of Cγ = β 2 γ/(1 − γ). In Lemma 11, we defined (recall (59))

η step 7 = 8δ1 ωLkΛ/Qk∞ ≤ 384β 2 γ1 + 48τ1row + 8βLξ ωLkΛ/Qk∞ .
(72)

By (39), 384β 2 γ1 ωLkΛ/Qk∞ ≤ Imin /(8Q). Take
τ1row =

Imin
,
384ωLkΛk∞

ξ=

col
Imin ∧ Imin
,
64βω(K ∨ L)2 (kΛk∞ ∨ kΓk∞ )

u=

Imin
,
8Q

(73)

where u is the parameter in (60). Then from (72) we have
η step 7 ≤

Imin Imin Imin
3Imin
+
+
=
,
8Q
8Q
8Q
8Q

η step 7 + u ≤

Imin
.
2Q

Hence Lemma 11 implies that on event P,
(q)

Misk ỹstep 7 , y (q)



row
≤ γ2k
:=

X
r6=k

 I
 I 
Imin 
min
kr
exp −
+
≤ K exp −
,
Q
2Q
2Q

∀q ∈ ZQ

(74)

fails with probability at most Q(e−u + 2Q p3 (τ1row )). By (A40 ), Q log Q = o(Imin ), hence
Qe−u = o(1). By (A3),
 nΛ
row 
min h1 (τ1 )
Q2 p3 (τ1row ) = Q2 KL exp −
32Q2 βK


2
nImin
≤ Q2 KL exp −
= o(1)
256(3842 )Q2 βKL2 ω 3 kΛk∞

33

(75)

where we have used the definition (56) of p3 , h1 (τ ) ≥ τ 2 /8 for τ ≤ 1 and kΛk∞ /Λmin ≤ ω.
col )/(K ∨ L) → ∞, hence eventually (I
col
Moreover, (A40 ) implies (Imin ∧ Imin
min ∧ Imin )/(K ∨ L) ≥ 1
which gives

P(Pc ) = p2 (ξ) = 2Q(K + L) exp −(n ∧ m)ξ 2 /Q


(76)
n∧m
= o(1)
≤ 2Q(K + L) exp − 2 2 2
2
2
64 Qβ ω (K ∨ L) (kΛk∞ ∨ kΓk∞ )
where the last implication follows from (A3).
First column LR classifier.

We can apply a similar argument to z̃step 8 . Let
col
Imin
,
768ωKkΓk∞

τ1col =

ucol =

col
Imin
,
16Q

col /(4Q).
with ξ defined as in (73). By (39), and a similar argument, we obtain η step 8 + ucol ≤ Imin
By Corollary 4, on event P,
(q)

Mis` z̃step 8 , z (q)



col
≤
≤ γ2`

X
r6=`

 I col 
 I col I col 
exp − `r + min ≤ L exp − min ,
2Q
4Q
4Q

∀q ∈ ZQ ,

(77)

col

fails with probability at most Q(e−u + 2Q p03 ), where Q2 p03 = Q2 p03 (τ1col ) = o(1) by (A3) and
col
Qe−u = o(1) by (A40 ), similar to how we argued for the row labels.
Second row LR classifier. Recalling γ2 from (64) and combining with (74) and (77),

 I 
 I col 
 β(K ∨ L)2 
min
γ2 ≤ max βK 2 exp −
, βL2 exp − min
=o
col )b
2Q
4Q
(Imin ∧ Imin

(78)

col → ∞. By Lemma 13 and (49),
for any b > 0, as Imin ∧ Imin


step 11
:= ηkr δ2 ; ω, β, m, Λ
ηkr



 1
5βL2 kΛk2∞
1
+ 1 − log Λmin
+ log 11ω
m
εkr − 2ω(1 + εkr ) δ2
2
=: T1 + T2 + T3 + T4 ,
(79)

= 21δ2 ωLkΛk∞ +

where we have called the four summands above T1 , . . . , T4 in the order they appear. We have
δ2 = 24Cγ2 + 6τ2 ≤ 48β 2 γ2 + 6τ2 by (67) and the assumption γ2 ≤ 21 . Then,
T1 ≤ 21(48)β 2 ωLkΛk∞ γ2 + 21(6)ωLkΛk∞ τ2 =: T11 + T12 .
For any b > 0, by (78)
 β 3 ω(K ∨ L)3 kΛk 

∞
T11 = O β 2 ωLkΛk∞ γ2 = o
.
col
[(Imin ∧ Imin )/Q]b

34

(80)

col )/Q]a ) for some a > 0 by (A40 ). Taking
Recall that we have βω(K ∨ L)kΛk∞ = o([(Imin ∧ Imin
b = 3a in (80), we obtain T11 = o(1). Letting τ2 = (ωLkΛk∞ )−1 , we have T12 = O(1), hence,
T1 = O(1). Recalling the probability bound in Lemma 13, we have by (A3)

 nΛ
min h1 (τ2 )
Qp3 (τ2 ) = QKL exp −
32Q2 βK
(81)


n
≤ QKL exp −
=
o(1)
256Q2 βKL2 ω 3 kΛk∞

where we have used h1 (τ ) ≥ τ 2 /8 for τ ≤ 1 and kΛk∞ /Λmin ≤ ω. Using (A3) again, T2 =
5βL2 kΛk2∞ /m = O(1).
Now let us consider the third piece
 T3 in (79). Recall that Jkr = LkΛk∞ /Ikr . By Lemma 20
−1
in Section 9.3.2, εkr ≥ 2 Jkr ∧ 1 . In bounding T1 , we have shown δ2 ωLkΛk∞ = O(1),
hence 2ωδ2 = O((LkΛk∞ )−1 ). Since Ikr → ∞ and Jkr ≥ 1/2 (see Remark 4), (LkΛk∞ )−1 =
−1
o Jkr
∧ 2 . Therefore, 2ωδ2 = o(εkr ∧ 1). As a result, 2ω(1 + εkr )δ2 = o(εkr ), hence


 
1
1 
eT3 := 11ω 1 +
.
(82)
=O ω 1+
εkr − 2ω(1 + εkr )δ2
εkr
√
−v
Finally, we let v = log Λmin . Since Λmin → ∞,
√ e = o(1). Applying Lemma 13, combined
with T1 + T2 = O(1), and (82), then for ζ = 1/ log Λmin = o(1),
 X


p

1 
1
Misk ybtop , ytop = O ω
1+
exp −Ikr − log Λmin + log Λmin
εkr
2
r6=k




 X

1
1
exp −Ikr −
− ζ log Λmin
=O ω
1+
εkr
2
r6=k

col

fails w.p. ≤ 2(e−v + 2Q p3 (τ2 )) + P(Pc ) + Q(2e−u∧u + Q(p3 (τ1row ) + p03 (τ1col ))) = o(1). When
we swap Atop and Abottom and repeat the algorithm, the same misclassification rate holds. The
proof of Theorem 2 is complete.

7.2

Proof of Theorem 1

We proceed by stating a few lemmas. The proofs are deferred to Appendix A.2.
P
2
2
Lemma 14.
`∈[L] (λr` − λk` ) ≥ 2Λmin Ikr . As a consequence, Λ∧ ≥ 2Λmin Imin .

Combining Lemma 14 with Theorem 3, and noting that kΛk∞ /Λ2∧ ≤ ω/(2Imin ) as a consequence of the lemma, we obtain the following guarantee for spectral clustering in terms of the
information matrix (Ikr ):

Corollary 5. Consider the spectral algorithm given in Algorithm 4, assume that for a sufficiently small C1 > 0,
β 2 ωKL(K ∧ L) α
≤ C1 (1 + κ)−2 .
2Imin
Then the algorithm outputs estimated row labels ỹ satisfying w.h.p.
Mis(ỹ, y) ≤

(1 + κ)2 ωβL(K ∧ L)α
.
2C1 Imin
35

(83)

We next modify Corollary 5 to be applicable on sub-blocks:
Lemma 15 (Spectral clustering on subblocks). Suppose (A3) holds, and we assume for a
sufficiently small C1 > 0,
6Qβ 2 ω 2 KL(K ∧ L) α
≤ C1 (1 + κ)−2 .
Imin

(84)

Using Algorithm 4 in Step 3 of Algorithm 3, w.h.p., the misclassification rate of ỹ (q) satisfies
Mis(ỹ (q) , y (q) ) ≤

3Q(1 + κ)2 ω 2 βL(K ∧ L)α
C1 Imin

∀q ∈ [Q].

A similar result holds for misclassification rate of the spectral clustering for column labels,
with appropriate modifications.
Proof of Theorem 1. Assumption (A4) implies (84), eventually as Imin → ∞. Letting γ1row and
γ1col be bounds on the misclassification rates of the spectral clustering algorithms in steps 3
and 4, we can take, by Lemma 15 (and its column counterpart), w.h.p.
γ1row = O

 QωβL(K ∧ L)α 
Imin

,

γ1col = O

 QωβK(K ∧ L)α−1 
col
Imin

.

That is, by the end of step 4, w.h.p., Mis(ỹ (q) , y (q) ) ≤ γ1row and Mis(z̃ (q) , z (q) ) ≤ γ1col for all
q ∈ [Q]. Since the matching step increases the misclassification rate by at most a factor of
2, the same bounds hold for the overall initial labels at step 6. Taking γ1 = γ1row ∨ γ1col , we
observe that in order to satisfy condition (39) of Theorem 2, it is enough to have
 1

col
Imin
Qωβ(K ∨ L)2 (α ∨ α−1 )
Imin
=
o
∧
col
β 2 ω LkΛk∞ KkΓk∞
Imin ∧ Imin
which holds if we require the stronger condition

 1
col
Imin ∧ Imin
Qωβ(K ∨ L)2 (α ∨ α−1 )
.
=
o
col
β 2 ω (K ∨ L)(kΛk∞ ∨ kΓk∞ )
Imin ∧ Imin
But this latter condition is satisfied by assumption (A4). Thus, the assumptions of Theorem 2
hold with high probability, and so is its result. The proof is complete.

7.3

Proof of Corollary 1

Proof of Corollary 1. From the proof of Theorem 1, we have that
X 



1 
1
Misk yb, y = O
ω 1+
exp −Ikr − log Λmin + v
εkr
2

(85)

r6=k

fails with probability at most 2e−v + o(1). First, we show that

1  −1/2
χr := ω 1 +
Λ
= o(1),
εkr min
36

uniformly in r.

(86)


−1
By Lemma 20 in Section 9.3.2, εkr ≥ 2 Jkr
∧ 1 . Hence,
1+

1
3
1
≤ 1 + (Jkr ∨ 1) ≤ (Jkr ∨ 1) ≤ 3Jkr
εkr
2
2

√
using 2Jkr ≥ 1 (see Remark 4). Thus, to show (86), it is enough to show ωJkr / Λmin = o(1).
Using ω −1 kΛk∞ ≤ Λmin , we have
−1/2

ωJkr Λmin =

1/2

1/2

kΛk∞
Lω 3/2 kΛk∞
Lω 3/2 kΛk∞
−1/2
= o(1)
L ωΛmin ≤
≤
Ikr
Ikr
Imin

2 ) which is implied by (A4). Thus, we have (86),
where the last equality is by ω 3 L2 kΛk∞ = o(Imin
√
i.e., χ := maxr χr = o(1), as desired. Now, let 2v = − log χ. It follows that e−v = χ = o(1),
and we have
 X
√ X
X




exp −Ikr + v = O χ
exp −Ikr = o
exp −Ikr
Misk yb, y = O χ
r6=k

r6=k

r6=k

completing the proof.

7.4

Proof of Example 1

Proof of Example 1. Without loss of generality assume a > b so that εkr = a/b − 1. Also,
Λmin ≥ b/(βK). By (85), which holds in the general case, we have that
Misk

X b

 b 

1
yb, y = O
ω
exp −Ikr − log
+v
a
2
βK
r6=k
 (√a − √b)2

p
3/2 −1/2
βωK b
exp −
+v
=O
βK


fails with probability at most 2e−v + o(1). Assumption βω 2 K 3 = o(b) implies
p
χ := βωK 3/2 b−1/2 = o(1).

Letting 2v = − log χ, the rest of the proof follows similar to that of Corollary 1.

8

Simulations

We provide some simulation results to corroborate the theory.
model of Section 2.1 with the following connectivity matrix

1 2 3 4

α
2 3 4 5
log(mn)
√
P =C
B, B = 
3 4 5 6
mn
4 5 6 1

We generate from the SBM
5
6
1
2


6
1
.
2
3

(87)

Note that B does not have any clear assorative or dissortative structure. We let n = Kn0 and
m = Ln0 , and we vary n0 . All clusters (both row and column) will have the same number
of nodes n0 . By changing α, we can study different regimes of sparsity. In particular, when
37

C = 1.00,

C = 1.00,

= 0.75, =1.00, K=4, L=6

= 0.75, =1.00, K=4, L=6

0.8

-0.8
0.7

log miss. (overall)

-1

NMI (Overall)

0.6

Spectral

0.5

Oracle
Hard

-1.2

-1.4

-1.6

Soft

0.4

-1.8
0.3

-2

200

400

600

102

800 1000 1200 1400 1600 1800 2000

# of nodes per cluster

103

# of nodes per cluster

(a)

(b)

Figure 3: Plots of (a) the (overall) NMI and (b) the corresponding log. misclassification rate, for the
SBM model with connectivity matrix (87). The four algorithms considered are the Spectral clustering
of Algorithm 4, Soft and Hard versions of Algorithm 1 and the Oracle algorithm of Section 2.2.

α ∈ (0, 1), we are in the regime where weak recovery is possible but not exact (or strong)
recovery. We consider both the misclassification rate, and the normalized mutual information
(NMI) as measures of performance. NMI is a measure of accuracy which is between 0 and
1 (=perfect match). The NMI is quite sensitive to mismatch and tends to reveal discrepancies
between methods more clearly. Figure 3(a) shows the overall NMI versus n0 . Figure 3(b)
illustrates the corresponding log. misclassification rates.
We have considered four algorithms: (1) Spectral: the spectral clustering of Algorithm 4.
(2) Soft: Algorithm 1 with flat prior, no inner loop and no conversion to hard labels. (3) Hard:
Algorithm 1 with flat prior, no inner loop and conversion to hard labels after each label computation. (4) Oracle: The oracle classifier discussed in Section 2.2 and Remark 2: Assuming
the knowledge of z and Λ, we obtain yb by the likelihood ratio classifier, and similarly obtain
zb, assuming the knowledge of y and Γ.
Figure 3 shows the results for α = .75 (regime where no exact recovery is possible) and
C = 1. Both the soft and hard versions of Algorithm 1 are initialized with the spectral clustering
and both significantly improve over it. The soft version of Algorithm 1 also outperforms the
hard version as one would expect: soft labels carry more information between iterations. It
is also interesting to note that the slope for the log. misclassification rate of Algorithm 1
approaches that of the oracle (esp. clear for the soft version in Figure 3(b)) as predicted by
the theory. Simulation results for various other settings can be found in Appendix B, showing
qualitatively similar behavior.

9

Proofs of the main lemmas

In this section, we give the proof of the three main lemmas of Section 6.1. We we first give the
proofs of Lemma 1 and 3 in Section 9.1 and 9.2. The proof of Lemma 2(b) is more technical
and occupies the remainder of this section, including auxiliary results on the error exponents

38

and Poisson-binomial approximations, in Sections 9.3 and 9.4.
Throughout, we will use the following concentration inequality [GN15, p. 118]:
P
Proposition 2 (Prokhorov). Let S = i Xi for independent
centered variables {Xi }, each
P
bounded by c < ∞ in absolute value a.s. and let v ≥ i EXi2 , then

P S ≥ vt ≤ exp[−vhc (t)],

t ≥ 0,

where hc (t) :=

Same bound holds for P(S < −vt).

2c 
3
t log 1 + t .
4c
3

(88)

Note that hc (t)  t2 as t → 0 and hc (t)  t log t as t → ∞.

9.1

Proof of Lemma 1

Let us define the confusion matrix as R(z̃, z) ∈ [0, 1]L×L with entries
m

|j : z̃j = k, zj = `|
1 X
Rk` (z̃, z) =
.
1{z̃j = k, zj = `} =
m
m

(89)

j=1

We can similarly define Rk` (z, z̃). It is easy verify that R(z̃, z) = R(z, z̃)T . By definition (22)
of the (global) row mean parameters,
λk`0 (y, z̃) =

m X
L
X
j=1 `=1

Pk` 1{zj = `, z̃j = `0 } = mPk∗ R∗`0 (z, z̃).

(90)

To see (90), note
that since we are
 using true labels y in the first argument of λk`0 (y, z̃), the
1 P
averaging nk (y) i 1{yi = k} · · · over i, in the definition, is vacuous. That is, for any i with
P
yi = k, we have λk`0 (y, z̃) = j E[Aij ]1{z̃j = `0 }. We then further break this sum according to
column labels zj = ` to get (90).
Recall that n(z) is the vector of sizes of clusters in z and π(z) = n(z)/m is the corresponding
proportions. To simplify, let
N (z) := diag(n(z)),

Π(z) := diag(π(z)).

We have mIL = N (z)Π(z)−1 where IL is the L × L identity matrix, hence
mPk∗ R∗`0 (z, z̃) = Pk∗ N (z) Π(z)−1 R∗`0 (z, z̃) = λk∗ (y, z) Π(z)−1 R∗`0 (z, z̃)
using (2). Let use define
U (z, z̃) := Π(z)−1 R(z, z̃).
Since π(z) contains the row sums of R∗`0 (z, z̃), U (z, z̃) is the row-normalized confusion matrix,
i.e. U = (Rk` /Rk+ ). We have
λk`0 (y, z̃) = λk∗ (y, z) U∗`0 (z, z̃),

39

(91)

and its matrix version Λ(y, z̃) = Λ(y, z) U (z, z̃). We can similarly define U (ỹ, y) = Π(ỹ)−1 R(ỹ, y).
Recalling definition (22), and some algebra gives
n
1 X X
λk0 `0 (ỹ, z̃) =
λk`0 (y, z̃)1{yi = k, ỹi = k 0 }
nk0 (ỹ)
i=1 k∈[K]

=

X
1
λk`0 (y, z̃) |i : yi = k, ỹi = k 0 |,
nk0 (ỹ)
k∈[K]

P
where to get the first equality one further breaks the sums over k 1{yi = k} and use the
expression for λk`0 (y, z̃) in the comments after (90). Using the definition of the confusion
matrix in (89), adapted to row labels, and the definition of U , we have
λk0 `0 (ỹ, z̃) =

1
Rk0 ∗ (ỹ, y) λ∗`0 (y, z̃) = Uk0 ∗ (ỹ, y) λ∗`0 (y, z̃),
πk0 (ỹ)

(92)

or compactly Λ(ỹ, z̃) = U (ỹ, y)Λ(y, z̃). We also define a column-normalized confusion matrix,
V (z, z̃) := R(z, z̃)Π(z̃)−1 .
Lemma 16. (A2) and (B3) imply
max [1 − Ukk (ỹ, y)] ≤ γ,

(B3.1)

k

max [1 − U`` (z, z̃)] ≤ γ,
`

and

(B3.2)

max [1 − V`` (z, z̃)] ≤ γ.

(B3.3)

`

Proof. Without loss of generality, assume that the optimal permutation matching ỹ to y is
identity, and similarly for z̃ to z. By definition, 1 − Ukk (ỹ, y) is the misclassification rate
withing the kth community of ỹ, hence
1 − Ukk (ỹ, y) =

|i : ỹi = k, yi 6= k|/n
|i : ỹi = k, yi 6= k|/n
=
.
|i : ỹi = k|/n
|i : ỹi = k, yi 6= k|/n + |i : ỹi = yi = k|/n

Recall that we can write (see Section 2.3)
Mis(ỹ, y) =

1
1X
1X
|i : ỹi 6= yi | =
|i : ỹi = k, yi 6= k| =
|i : ỹi 6= k, yi = k|.
n
n
n
k

(93)

k

Then, (B3) and the second equality in (93) implies |i : ỹi = k, yi 6= k|/n ≤ γ/(βK), while the
third equality in (93) gives |i : ỹi = yi = k|/n ≥ πk (y) − γ/(βK). Letting f (x) = x/(x + 1),
 |i : ỹ = k, y 6= k|/n 
γ/(βK)
γ
i
i
1 − Ukk (ỹ, y) = f
≤
=
≤γ
|i : ỹi = yi = k|/n
γ/(βK) + πk (y) − γ/(βK)
πk (y)βK
where the first inequality is by monotonicity of f , and the last by (A2). This proves (B3.1).
Similarly, 1 − U`` (z, z̃) is the misclassification rate within the `th community of z, i.e.,
Mis` (z̃, z), hence
1 − U`` (z, z̃) =

|j : zj = `, z̃j 6= `|/m
γ
=
≤ γ,
π` (z)
π` (z)βL
40

proving (B3.2). The same bound holds for 1 − U`` (z̃, z) by an argument similar to that used
for Ukk (ỹ, y). To prove (B3.3), we observe
U (z̃, z) = Π(z̃)−1 R(z̃, z) = Π(z̃)−1 R(z, z̃)T = [R(z, z̃)Π(z̃)−1 ]T = V (z, z̃)T
hence 1 − V`` (z, z̃) = 1 − U`` (z̃, z) ≤ γ. All statements are true for any k ∈ [K] and ` ∈ [L].
9.1.1

Proof of Lemma 1(a)

For the lower bound, by (91) and (B3.2),
λk`0 (y, z̃) = λk∗ (y, z)U∗`0 (z, z̃) ≥ λk`0 (y, z)U`0 `0 (z, z̃)

≥ (1 − γ)λk`0 (y, z) ≥ λk`0 (y, z) − Cγ kΛk∞

where the last inequality is by γ ≤ Cγ and λk`0 (y, z) ≤ kΛk∞ . For the upper bound, we write
X
λk∗ (y, z)U∗`0 (z, z̃) = λk`0 (y, z)U`0 `0 (z, z̃) +
λk` (y, z)U``0 (z, z̃).
`6=`0

The first term obviously satisfies λk`0 (y, z) U`0 `0 (z, z̃) ≤ λk`0 (y, z), hence
X
λk`0 (y, z̃) − λk`0 (y, z) ≤
λk` (y, z)U``0 (z, z̃).

(94)

`6=`0

By (B3.3), for every `0 ∈ [L],
π`0 (z) ≥

1
|j : zj = z̃j = `0 | = π`0 (z̃)V`0 `0 (z, z̃) ≥ (1 − γ)π`0 (z̃).
m

(95)

By (A2), for every `0 and `, we have π`0 (z) ≤ β 2 π` (z), hence
U``0 (z, z̃) =

π`0 (z̃)
β2
1
R``0 (z, z̃) =
V``0 (z, z̃) ≤
V``0 (z, z̃).
π` (z)
π` (z)
1−γ

Combining with (94)
λk`0 (y, z̃) − λk`0 (y, z) ≤

β2 X
β2γ
λk` (y, z)V``0 (z, z̃) ≤
kλk∗ (y, z)k∞
1−γ
1−γ
0

(96)

`6=`

where the last inequality is by (B3.3) and that V is column normalized. This proves the upper
bound, and completes the proof of kΛ(y, z̃) − Λk∞ ≤ Cγ kΛk∞ . Since we assume Cγ ≤ 1, it
follows that kΛ(y, z̃)k∞ ≤ 2kΛk∞ .
9.1.2

Proof of Lemma 1(b)

Recalling (92), we have
λk0 `0 (ỹ, z̃) = Uk0 ∗ (ỹ, y)λ∗`0 (y, z̃) = Uk0 k0 (ỹ, y)λk0 `0 (y, z̃) +

41

X

k6=k0

Uk0 k (ỹ, y)λk`0 (y, z̃).

By (B3.1), the first term is bounded as
(1 − γ)λk0 `0 (y, z̃) ≤ Uk0 k0 (ỹ, y) λk0 `0 (y, z̃) ≤ λk0 `0 (y, z̃)
and the second term as
0 ≤

X

k6=k0

Uk0 k (ỹ, y)λk`0 (y, z̃) ≤ γkλ∗`0 (y, z̃)k∞

P
recalling that U is row normalized hence k6=k0 Uk0 k = 1 − Uk0 k0 ≤ γ, by (B3.1). Combining
the two bounds, we have

 

λk0 `0 (ỹ, z̃) − λk0 `0 (y, z̃) ∈ −γλk0 `0 (y, z̃), 0 + 0, γkλ∗`0 (y, z̃)k∞


⊆ kλ∗`0 (y, z̃)k∞ − γ, γ

showing that kΛ(ỹ, z̃) − Λ(y, z̃)k∞ ≤ γkΛ(y, z̃)k∞ . Combining with kΛ(y, z̃)k∞ ≤ 2kΛk∞ from
part (a) of the lemma, we have the first assertion of part (b). The second assertion follows
from γ ≤ 1/2 and part (a) by triangle inequality. (Note that assumption 6Cγ ω ≤ 1 in fact
implies γ ≤ 1/6 since β, ω ≥ 1 and γ ≤ Cγ .)
9.1.3

Proof of Lemma 1

Recalling definitions of λ̂k` and λk` (ỹ, z̃) from (21) and (22), we have
n X
m

 X

nk (ỹ) λ̂k` − λk` (ỹ, z̃) =
Aij − E[Aij ] 1{ỹi = k, z̃j = `}
i=1 j=1

P
which is of the form S = ij Xij with independent centered terms Xij = Aij − E[Aij ] with
P
P
P
2
|Xi | ≤ 1 and
ij EAij = nk (ỹ)λk` (ỹ, z̃). Note that the sums
ij var(Aij ) ≤
ij EXij =
in these expressions run over {(i, j) : ỹi = k, z̃j = `}. Applying the two-sided version of
Proposition 2, with v = nk (ỹ)λk` (ỹ, z̃), t = τ and c = 1, we have


P λ̂k` − λk` (ỹ, z̃) > λk` (ỹ, z̃) τ = P nk (ỹ) λ̂k` − λk` (ỹ, z̃) > nk (ỹ)λk` (ỹ, z̃) τ

≤ 2 exp −nk (ỹ)λk` (ỹ, z̃)h1 (τ ) .
Applying union bound over (k, `) ∈ [K] × [L], and using part (b) of this lemma, we have
kΛ̂ − Λ(ỹ, z̃)k∞ ≤ τ kΛ(ỹ, z̃)k∞ ≤ 4τ kΛk∞ with probability at least

1 − 2KL − min nk (ỹ) min λk` (ỹ, z̃) h1 (τ ) .
k

k,`

We have nk (ỹ) ≥ nπk (y)(1 − γ) ≥ n(βK)−1 /2 using (B3.1), (A2) and γ ≤ 1/2; see (95).
Similarly, since kΛ(ỹ, z̃) − Λk∞ ≤ 3Cγ kΛk∞ , we have
min λk` (ỹ, z̃) ≥ Λmin − 3Cγ kΛk∞ ≥ Λmin (1 − 3Cγ ω) ≥ Λmin /2.
k,`

42

9.2

Proof of Lemma 3

Let bi∗ = bi∗ (z̃). Recall (44), (45) and (46), and let
Sbk = Sk (b, Λ̂),

Zbik = Zik (bi∗ , Λ̂),

Ybikr = Zik (bi∗ , Λ̂).

For any event A and random variable X, let us write E[X; A] := E[X1A ]. Consider the
following event: A := {Λ̂ ∈ BΛ (δ)}. Pick some i ∈ [N ] with yi = k. Then,
[





E Sbk ; A = E Zbik ; A = P
Ybikr ≥ 0 ∩ A
r6=k

≤

≤
≤

X
r6=k

X
r6=k

X
r6=k

P Yikr (bi∗ (z̃), Λ̂) ≥ 0, Λ̂ ∈ BΛ (δ)




P ∃Λ̃ ∈ BΛ (δ), Yikr (bi∗ (z̃), Λ̃) ≥ 0

exp −Ikr + η 0 =: pk

where the last inequality follows from Lemma 2 with ηkr defined there. Using Markov inequality


P Sbk ≥ tpk ≤ P {Sbk ≥ tpk } ∩ A + P(Ac )
1
E[Sbk ; A]
+ P(Ac ) ≤ + P(Ac ).
≤
tpk
t

for any t > 0. The version of Markov inequality used follows from (pointwise) inequality:
1{X≥u} 1A ≤ (X1A )/u. Taking t = eu complete the proof.

9.3

Error exponents

We start by obtaining a bound on the error exponent (i.e., the negative logarithm of the
probability of error) for binary hypothesis testing in an exponential family. This result is a
generalization of the result that appears in [AS15], and is proved by the same technique. The
result (and the technique inspired by [AS15]) is interesting since it provides a bound different
than the classical Chernoff bound on the error exponent [Che52]; see also [Ver86] and [CT06,
Theorem 11.9.1]. This leads for example to a sharper control for the case of Poisson hypothesis
testing. We start with the result for a general exponential family and then in Section 9.3.2
specialize to the case of interest in this paper, the Poisson family.
9.3.1

General exponential family

Let π(t; γ) denote the density of a 1-dimensional standard exponential family w.r.t. to some
measure ν on R:

π(t; γ) = h(t) exp γt − A(γ) .
(97)
We consider distributions on RL that are products of these distributions, having density:
p(x; θ) =

L
Y

`=1

π(x` , θ` ),

x = (x` ) ∈ RL ,

θ = (θ` ) ∈ RL

with respect to µ = ν ⊗L (L-fold product measure whose coordinate measures are all ν).
43

(98)

Proposition 3. Let pr (x) := p(x; θr ), r = 0, 1 be two exponential family densities on RL
(relative to µ = ν ⊗L ) as defined in (97) and (98). Assume that ν is either the Lebesgue
measure on R or the counting measure on Z, and that θ0 6= θ1 . For s ∈ (0, 1), let


θs` = (1 − s)θ0` + sθ1` , and Is` = (1 − s)A(θ0` ) + sA(θ1` ) − A(θs` ),
(99)
as well as Is :=

PL

`=1 Is` ,

C(α) :=

T := {` : θ0` 6= θ1` } and
(
Z
e−α|t| dν(t) =

2
α
1+e−α
1−e−α

≤

2
1−e−α

ν is Lebesgue,
ν is counting.

(100)

Consider testing p0 against p1 using the likelihood ratio test based on a single observation. Let
pr be the probability of error under pr for r = 0, 1. Then, the sum of the error probabilities is
bounded as
h

i
Pe,0 + Pe,1 ≤ inf inf
e−Is kπ( · ; θs` )k∞ C min(s, 1 − s)|θ0` − θ1` | .
(101)
`∈T s∈(0,1)

Remark 8. The proof goes through for any translation invariant measure ν (e.g., a Haar
measure) with an appropriate constant C(α). It also goes through if we replace t in (97) with a
general sufficient statistic φ(t),Ras long as (1) φ is surjective from the support of the exponential
family to R and (2) C(α) = e−α|φ(t)| dν(t) < ∞ for all α > 0 and (3) φ has a measurable
inverse.
Remark 9. Let s∗ be the maximizer of s 7→ Is . Then, noting that α →
7 C(α) is decreasing,
Proposition 3 implies


Pe,0 + Pe,1 ≤ exp −Is∗ + log kπ( · ; θs∗ ` )k∞ + log C(α∗ ) , where,
(102)
α∗ = min(s∗ , 1 − s∗ ) max |θ0` − θ1` |.

(103)

`∈[L]

The bound is an improvement over the Chernoff bound if log kπ( · ; θs∗ ` )k∞ is negative and
log C(α∗ ) is controlled. This is the case for the Poisson distribution as we show in the sequel.
9.3.2

Poisson case

The Poisson case corresponds to (97) with γ = log λ, h(t) = (1/t!)1{t ≥ 0}, ν = the counting
measure and A(log λ) = λ. Letting θs` = log λs` for all s ∈ [0, 1], we have from (99)


s
λs` = λ1−s
Is` = (1 − s)λ0` + sλ1` − λs` .
0` λ1` ,

We also note that |θ0` − θ1` | = | log(λ0` /λ1` )|. Let us define
s∗ = argmax Is ,
s∈(0,1)

and,

I ∗ = max Is ,
s∈(0,1)

where

Is =

L
X

Is`

`=1

We will assume
λ0` /λ1` ∈ [1/ω, ω], ∀` ∈ [L], for some ω > 1.
The following lemma shows that s∗ stays away from the boundary:
44

(104)

1
,1 −
Lemma 17. Assuming (104), we have s∗ ∈ [ 2ω

α∗

1
2ω ].

Proof of this lemma and subsequent results appear in Section A.4. From (103), we have
= min(s∗ , 1 − s∗ ) max` | log(λ0` /λ1` )| in the Poisson case. Defining
λ
λ1` 
1
0`
ε01 := ε01 (Λ) := max
∨
− 1, α01 :=
log(1 + ε01 )
(105)
λ0`
2ω
`∈[L] λ1`

we note that α∗ = min(s∗ , 1 − s∗ ) log(1 + ε01 ), hence Lemma 17 implies α∗ ≥ α01 , that is,
C(α∗ ) ≤ C(α01 ) in (102), where C(·) has the form given in (100) for the counting measure,
i.e.,
C(α01 ) =


(a)  4
2
1 + e−α01
≤
≤
+
3
ω
1 − e−α01
1 − e−α01
ε01

(106)

where the last inequality is by the following lemma:
Lemma 18. Inequality (a) in (106) holds.
Next we bound the maximum of the density:

Lemma 19 ([HLC60]). Let π(t; log λ) = e−t (λt /t!)1{t ≥ 0} be the desnity of the Poisson
family. Then, for all λ > 0,

1  1
√
kπ( · ; log λ)k∞ ≤ 1 +
.
12λ
2πλ

In particular kπ( · ; log λ)k∞ ≤ exp − 12 log λ for λ ≥ 0.056.

Combining Lemmas 17, 18 and 19, we have the following corollary which gives the following
overall bound on the error exponent:

Corollary 6. Consider testing two Poisson vector models with mean vectors given by the rows
of Λ = [λ0∗ ; λ1∗ ] ∈ R2×L
+ , satisfying (104). Let Λmin = minr` λr` . Then, the sum of the error
probabilities for the likelihood ratio test is bounded as



 4
1
Pe,0 + Pe,1 ≤ ω
+ 3 exp −I ∗ − log Λmin .
(107)
ε01
2
We also have the following general lower bound on ε01 in terms of the information I ∗ :

Lemma 20. Let Λ = [λ0∗ ; λ1∗ ] ∈ R2×L
+ . There exists ` ∈ [L] such that
log
which implies ε01 ≥ min



2I ∗
LkΛk∞ , 2


1
8I ∗ 
λ0`
≥ log 1 +
,
λ1`
2
LkΛk∞

.

Although the bound in Lemma 20 holds without any further assumption, it is not always
tight. The difference in our two sets of results, namely (13) and (14) is due to using the sharper
bound (107) versus replacing ε01 with its universal lower bound.
Lemma 21 (Perturbation of ε01 ). Suppose Λ0 ∈ BΛ (δ), and let ε001 = ε01 (Λ0 ) and ε01 = ε01 (Λ)
as in (105), and assume (104). Then ε001 ≥ ε01 − 2ω(1 + ε01 )δ.
45

9.4

Approximation results for Lemma 2(b)

Let us collect some approximation lemmas that will be used in the proof of Lemma 2(b). The
proofs can be found in Appendix A.4. We write pmf for the probability mass functions. We
recall that a Poisson-binomial variable with parameter (p1 , . . . , pn ) is one that can be written as
P
n
i=1 Xi where Xi ∼ Ber(pi ), independent over i = 1, . . . , n. We write pmf for the probability
mass function.
Lemma 22 (Poisson-binomial approximation). Let ϕ(x; λ) be the pmf of a Poisson variable
with mean λ, and P
let ϕ(x,
e p) be the pmf of a Poisson-binomial variable with parameters p =
(p1 , . . . , pn ) where nj=1 pj = λ. Let p∗ := maxj∈[n] pj . Then,
ϕ(x;
e p)
∗
≤ exp ,
ϕ(x; λ)

∀x ∈ Z+ .

This result immediately extends to the comparison between vector versions of the two
distributions:
(`)

(`)

Corollary 7 (Poisson-binomial approximation). Let p(`) = (p1 , . . . , pn` ) ∈ [0, 1]n` be a vector
P ` (`)
of probabilities for each ` ∈ [L] and let λ(`) = ni=1
pi ∈ R+ . Let
e (p(1) , . . . , p(L) )) :=
Φ(x,

L
Y

`=1

ϕ(x
e ` ; p(`) ),

for each x = (x1 , . . . , xL ) ∈ ZL
+

be the pmf of a vector Poisson-binomial variable, and Φ(x, (λ(1) , . . . , λ(L) )) =
be the corresponding vector Poisson pmf. Then, we have
L

 X
e (p(1) , . . . , p(L) ))
Φ(x,
∗
,
≤
exp
p
x
`
Φ(x, (λ(1) , . . . , λ(L) ))
`=1

(108)

QL

`=1 ϕ(x` ; λ

(`) )

∀x ∈ ZL
+,

(`)

where p∗ = max{pi : i ∈ [n` ], ` ∈ [L]}.

Lemma 23 (Poisson likelihood approximation). Suppose max(|λ1 − λ|, |λ2 − λ|) ≤ ρ ≤ 13 λ,
then for any x ∈ Z+ , we have
 3ρx

φ(x; λ1 )
≤ exp
+ 2ρ .
φ(x; λ2 )
λ
P
Pm
Lemma 24 (Degree Truncation). Let bi+ =
`∈[L] bi` =
j=1 Aij be the degree of (row)
node i. Then,
P(bi+ > 5LkΛk∞ ) ≤ exp(−3LkΛk∞ ).

9.5

Proof of Lemma 2(b)

Fix i ∈ [n] such that yi = k, and z̃ ∈ [K]n and let bi∗ = bi∗ (z̃). Throughout, let Λ0 = (λ0k` ) :=
Λ(y, z̃) which belongs to BΛ (δ) by assumption. Denoting the kth row of Λ0 as λ0k∗ , we have
E[bi∗ ] = λk∗ . For r 6= k ∈ [K], i such that yi = k and Λ̃ ∈ BΛ (δ),
Yikr (bi∗ , Λ̃) =

L
X
`=1

bi` log

L h
i
X
λr` + ρ
λ̃r`
+ λ̃k` − λ̃r` ≤
bi` log
+ λk` − λr` + 2ρ := Y ∗
λk` − ρ
λ̃k`
`=1

46

where ρ := δkΛk∞ is the radius of BΛ (δ). Hence,
P(∃Λ̃ ∈ BΛ , Yikr (Λ̃) ≥ 0) ≤ P(Y ∗ ≥ 0) = P(bi∗ ∈ F ),
where we have defined (recalling the definition of Ψ from (32)):
n
o
F := x ∈ ZL
:
Ψ(x;
λ
+
ρ
|
λ
−
ρ)
≥
−2Lρ
.
r∗
k∗
+
Degree truncation.

Let bi+ =

P

`∈[L] bi`

=

Pm

j=1 Aij

be the degree of (row) node i, and

L
n
o
X
E = x ∈ ZL
:
x
≤
5LkΛk
∞ .
`
+

(109)

`=1

Using Lemma 24, we have P(bi∗ ∈
/ E) ≤ exp(−3LkΛk∞ ), which is faster than the rate we want
to establish. Hence, for the rest of the proof it is enough to work on {bi∗ ∈ E}. We have the
following two approximations on this event:
Poisson-binomial approximation.

Recall that P is the connectivity matrix and we have,

kP k∞ ≤

βLkΛk∞
kΛk∞
≤
mini ni (z)
m

(110)

where the first inequality follows from definition
Pm of Λ in (2), and the second from assumption (A2). We note that bi` = bi` (z̃) =
j=1 Aij 1{z̃j = `} as defined in (26), follows a
Poisson-binomial distribution. In order the describe the parameters of this distribution, let us
introduce the following notation
lab` (z̃) := (zj : j ∈ [m] such that z̃j = `},
that is, the vector of true labels associated with nodes in the `th cluster of z̃. Then, Pk,lab` (z̃) =
(Pk,zj : j ∈ [m] s.t. z̃j = `) ∈ Rn` (z̃) is the probability vector associated with the Poissonbinomial distribution of bi` . Also, let
lab(z̃) := (lab1 (z̃), . . . , labL (z̃)),

and Pk,lab(z̃) := (Pk,lab1 (z̃) , . . . , Pk,labL (z̃) ).

Then, we can say that bi∗ = bi∗ (z̃) is a product Poisson-binomial distribution with parameter
e Pk,lab(z̃) ) as defined in (108). We also note that
Pk,lab(z̃) . In particular, bi∗ has pmf Φ(x;
0
E[bi∗ (z̃)] = λk∗ (y, z̃) =: λk∗ . It follows from Corollary 7, noting that kPk,lab(z̃) k∞ ≤ kP k∞
combined with (110),
L


 5βL2 kΛk2 
e Pk,lab(z̃) )
X
Φ(x;
∞
≤
exp
kP
k
x
≤
exp
=: ζ1 ,
`
k,lab(z̃) ∞
Φ(x; λ0k∗ )
m
`=1

47

∀x ∈ E.

Poisson likelihood approximation. Recall that ρ = δkΛk∞ . Since by assumption, ωδ ≤ 13 ,
∞
1
0
0
we have ρ ≤ kΛk
3ω ≤ 3 Λmin . Recall that by assumption Λ = (λk` ) ∈ BΛ (δ). By Lemma 23,
 3ρx



Y
Φ(x; λ0k∗ )
15ρ
`
exp
LkΛk∞
≤
+ 2ρ ≤ exp 2Lρ +
Φ(x; λk∗ − ρ)
λk`
Λmin
`∈[L]


≤ exp 17ωLρ =: ζ2 , ∀x ∈ E.

e are measures as well, thus, for example,
WithPsome abuse of notation, we treat Φ and Φ
Φ(E) = x∈E Φ(x). Then, we have
e E ∩ F ; Pk∗
P(bi∗ ∈ E ∩ F ) = Φ



≤ ζ1 Φ E ∩ F ; λ0k∗




≤ ζ1 ζ2 Φ E ∩ F ; λk∗ − ρ .

(111)

Thus, it is enough to bound Φ(F ; λk∗ − ρ) which gives a further upper bound. This quantity
is closely related to testing Poisson vector distributions with mean λk∗ − ρ and λr∗ − ρ against
each other. Let us write p0 (x) := Φ(x; λk∗ − ρ) and p1 (x) := Φ(x; λr∗ + ρ) and note that
Ψ( · ; λr∗ + ρ | λk∗ − ρ) = log(p1 /p0 ). We have
X

x∈F

Φ(x; λk∗ − ρ) =

n
o
p1 (x)
≥ −2Lρ
p0 (x) 1 log
p0 (x)
L

X

x ∈ Z+

o
n e2Lρ p (x)
1
≥1
p0 (x) 1
p0 (x)
x ∈ ZL
+


X
X

≤
min e2Lρ p1 (x), p0 (x) ≤ e2Lρ
min p1 (x), p0 (x) .
=

X

x ∈ ZL
+

(112)

x ∈ ZL
+

Let us define
L
X


Is (λ0 | λ1 ) =
sλ0` + (1 − s)λ1` − λs0` λ1−s
1` ,
`=1

We can now apply Corollary 6. Since
Corollary 6 by 2ω. It follows that

kΛk∞ +ρ
Λmin −ρ

≤

ωΛmin + 13 Λmin
2
Λ
3 min

λ0 , λ1 ∈ RL
+.

(113)

≤ 2ω, we need to substitute ω in




1
min p1 (x), p0 (x) ≤ ζ3 exp −Is (λr∗ + ρ | λk∗ − ρ) − log(Λmin − ρ)
2
L

X

x ∈ Z+


1
2 
≤ ζ3 exp −Is (λk∗ | λr∗ ) + 2ωLρ −
log Λmin + log
2
3
r


3
1
ζ3 exp −Is (λk∗ | λr∗ ) + 2ωLρ − log Λmin
(114)
≤ 8
2
2

where ζ3 = ω/(εkr − 2ω(1 + εkr )δ) + ω from Lemma 21, and the second line follows from the
following elementary inequality:

sρ 
(a − ρ)1−s (b + ρ)s ≤ a1−s bs + 1−s ≤ a1−s bs + ρω
b
48

assuming a/b ≤ ω. Note that Ikr = sups∈(0,1) Is (λr∗ | λk∗ ). Putting the pieces (111), (112)
and (114) together (and taking supremum over s) we have
r


3
1
ζ2 ζ3 e2Lρ+2ωLρ exp −Ikr − log Λmin .
P(bi∗ ∈ E ∩ F ) ≤ 8
2
2
We note that
log(ζ1 ζ2 e2Lρ+2ωLρ ) ≤ 17ωLρ +

5βL2 kΛk2∞
5βL2 kΛk2∞
+ 4ωLρ ≤ 21ωLρ +
=: log ζ4
m
m

It follows that

Finally we have



1
P(bi∗ ∈ E ∩ F ) ≤ ζ3 ζ4 exp −Ikr − log Λmin .
2

P(bi∗ ∈ F ) ≤ P(bi∗ ∈ E ∩ F ) + P(bi∗ ∈ E c )
r



3
1
ζ3 ζ4 exp −Ikr − log Λmin + exp −3LkΛk∞
≤8
2
2


1
≤ 11ζ3 ζ4 exp −Ikr − log Λmin ,
2

assuming that Ikr and Λmin are sufficiently large. Noting that by the definition of ηkr in the
statement of the theorem, ηkr = log(2ζ3 ζ4 ), the proof is complete.

Acknowledgement
We thank Yunfeng Zhang for helpful discussions.

References
[Abb17]

E. Abbe. Community detection and stochastic block models: recent developments. arXiv preprint arXiv:1703.10146 (2017).

[ABH16]

E. Abbe, A. S. Bandeira, and G. Hall. Exact recovery in the stochastic block
model. IEEE Transactions on Information Theory 62.1 (2016), pp. 471–487.

[ABKK17]

N. Agarwal, A. S. Bandeira, K. Koiliaris, and A. Kolla. Multisection in the
stochastic block model using semidefinite programming. Compressed Sensing and
its Applications. Springer, 2017, pp. 125–162.

[ACBL+13]

A. A. Amini, A. Chen, P. J. Bickel, E. Levina, et al. Pseudo-likelihood methods
for community detection in large sparse networks. The Annals of Statistics 41.4
(2013), pp. 2097–2122.

[AL14]

A. A. Amini and E. Levina. On semidefinite relaxations for the block model.
arXiv preprint arXiv:1406.5647 (2014).

[AL18]

A. A. Amini and E. Levina. On semidefinite relaxations for the block model. The
Annals of Statistics 46.1 (2018), pp. 149–179.
49

[AS15]

E. Abbe and C. Sandon. Community detection in general stochastic block models: Fundamental limits and efficient algorithms for recovery. Foundations of
Computer Science (FOCS), 2015 IEEE 56th Annual Symposium on. IEEE. 2015,
pp. 670–688.

[Ban15]

A. S. Bandeira. Random Laplacian matrices and convex relaxations. Foundations
of Computational Mathematics (2015), pp. 1–35.

[BC09]

P. J. Bickel and A. Chen. A nonparametric view of network models and Newman–
Girvan and other modularities. Proceedings of the National Academy of Sciences
106.50 (2009), pp. 21068–21073.

[BCCZ13]

P. Bickel, D. Choi, X. Chang, and H. Zhang. Asymptotic normality of maximum likelihood and its variational approximation for stochastic blockmodels.
The Annals of Statistics (2013), pp. 1922–1943.

[BLM15]

C. Bordenave, M. Lelarge, and L. Massoulié. Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs. Foundations of Computer Science (FOCS), 2015 IEEE 56th Annual Symposium on.
IEEE. 2015, pp. 1347–1357.

[CC00]

Y. Cheng and G. M. Church. Biclustering of expression data. Ismb. Vol. 8. 2000.
2000, pp. 93–103.

[CDP12]

A. Celisse, J.-J. Daudin, and L. Pierre. Consistency of maximum-likelihood and
variational estimators in the stochastic block model. Electronic Journal of Statistics 6 (2012), pp. 1847–1899.

[Che52]

H. Chernoff. A measure of asymptotic efficiency for tests of a hypothesis based on
the sum of observations. The Annals of Mathematical Statistics (1952), pp. 493–
507.

[Chv79]

V. Chvátal. The tail of the hypergeometric distribution. Discrete Mathematics
25.3 (1979), pp. 285–287.

[CT06]

T. M. Cover and J. A. Thomas. Elements of information theory. John Wiley &
Sons, 2006.

[Dhi01]

I. S. Dhillon. Co-clustering documents and words using Bipartite Co-clustering
documents and words using Bipartite Spectral Graph Partitioning. Proc of 7th
ACM SIGKDD Conf (2001), pp. 269–274.

[Dhi03]

I. S. Dhillon. Information-Theoretic Co-clustering (2003), pp. 89–98.

[DKMZ11]

A. Decelle, F. Krzakala, C. Moore, and L. Zdeborová. Asymptotic analysis of the
stochastic block model for modular networks and its algorithmic applications.
Physical Review E 84.6 (2011), p. 066106.

[FH16]

S. Fortunato and D. Hric. Community detection in networks: A user guide.
Physics Reports 659 (2016), pp. 1–44.

[Fis+13]

D. E. Fishkind et al. Consistent adjacency-spectral partitioning for the stochastic
block model when the model parameters are unknown. SIAM Journal on Matrix
Analysis and Applications 34.1 (2013), pp. 23–39.

50

[GLM17]

L. Gulikers, M. Lelarge, and L. Massoulié. A spectral method for community detection in moderately sparse degree-corrected stochastic block models. Advances
in Applied Probability 49.3 (2017), pp. 686–721.

[GLMZ16]

C. Gao, Y. Lu, Z. Ma, and H. H. Zhou. Optimal estimation and completion
of matrices with biclustering structures. Journal of Machine Learning Research
17.161 (2016), pp. 1–29.

[GMZZ16]

C. Gao, Z. Ma, A. Y. Zhang, and H. H. Zhou. Community Detection in DegreeCorrected Block Models. arXiv preprint arXiv:1607.06993 (2016).

[GMZZ17]

C. Gao, Z. Ma, A. Y. Zhang, and H. H. Zhou. Achieving optimal misclassification proportion in stochastic block models. The Journal of Machine Learning
Research 18.1 (2017), pp. 1980–2024.

[GN15]

E. Giné and R. Nickl. Mathematical foundations of infinite-dimensional statistical models. Vol. 40. Cambridge University Press, 2015.

[GV16]

O. Guédon and R. Vershynin. Community detection in sparse networks via
Grothendieck’s inequality. Probability Theory and Related Fields 165.3-4 (2016),
pp. 1025–1049.

[Har72]

J. A. Hartigan. Direct clustering of a data matrix. Journal of the american
statistical association 67.337 (1972), pp. 123–129.

[HLC60]

J. L. Hodges and L. Le Cam. The Poisson approximation to the Poisson binomial
distribution. The Annals of Mathematical Statistics 31.3 (1960), pp. 737–740.

[HWX16a]

B. Hajek, Y. Wu, and J. Xu. Achieving exact cluster recovery threshold via
semidefinite programming. IEEE Transactions on Information Theory 62.5 (2016),
pp. 2788–2797.

[HWX16b]

B. Hajek, Y. Wu, and J. Xu. Achieving exact cluster recovery threshold via
semidefinite programming: Extensions. IEEE Transactions on Information Theory 62.10 (2016), pp. 5918–5937.

[Jam15]

G. J. O. Jameson. A simple proof of Stirling’s formula for the gamma function.
The Mathematical Gazette 99.544 (2015), p. 68.

[Krz+13]

F. Krzakala et al. Spectral redemption in clustering sparse networks. Proceedings
of the National Academy of Sciences 110.52 (2013), pp. 20935–20940.

[LCJ14]

D. B. Larremore, A. Clauset, and A. Z. Jacobs. Efficiently inferring community
structure in bipartite networks. Physical Review E - Statistical, Nonlinear, and
Soft Matter Physics 90.1 (2014), pp. 17–22. arXiv: 1403.2933.

[LR13]

J. Lei and A. Rinaldo. Consistency of spectral clustering in sparse stochastic
block models. arXiv preprint. arXiv preprint arXiv:1312.2050 (2013).

[Mas14]

L. Massoulié. Community detection thresholds and the weak Ramanujan property. Proceedings of the forty-sixth annual ACM symposium on Theory of computing. ACM. 2014, pp. 694–703.

[MNS15]

E. Mossel, J. Neeman, and A. Sly. Consistency thresholds for the planted bisection model. Proceedings of the forty-seventh annual ACM symposium on Theory
of computing. ACM. 2015, pp. 69–75.

51

[MS16]

A. Montanari and S. Sen. Semidefinite programs on sparse random graphs and
their application to community detection. Proceedings of the forty-eighth annual
ACM symposium on Theory of Computing. ACM. 2016, pp. 814–827.

[MTSCO10]

S. C. Madeira, M. C. Teixeira, I. Sa-Correia, and A. L. Oliveira. Identification
of regulatory modules in time series gene expression data using a linear time
biclustering algorithm. Computational Biology and Bioinformatics, IEEE/ACM
Transactions on 7.1 (2010), pp. 153–165.

[NS01]

K. Nowicki and T. A. B. Snijders. Estimation and prediction for stochastic
blockstructures. Journal of the American statistical association 96.455 (2001),
pp. 1077–1087.

[PW17]

A. Perry and A. S. Wein. A semidefinite program for unbalanced multisection in
the stochastic block model. Sampling Theory and Applications (SampTA), 2017
International Conference on. IEEE. 2017, pp. 64–67.

[PZ17]

M. Pensky and T. Zhang. Spectral clustering in the dynamic stochastic block
model. arXiv preprint arXiv:1705.01204 (2017).

[RAL17]

Z. S. Razaee, A. A. Amini, and J. J. Li. Matched bipartite block model with
covariates. Preprint (2017). arXiv: 1703.04943.

[RCY11]

K. Rohe, S. Chatterjee, and B. Yu. Spectral clustering and the high-dimensional
stochastic blockmodel. The Annals of Statistics (2011), pp. 1878–1915.

[Roh15]

K. Rohe. Co-clustering for directed graphs : the Stochastic co-Blockmodel and
spectral algorithm Di-Sim (2015), pp. 1–39. arXiv: arXiv:1204.2296v2.

[RTJM16]

F. Ricci-Tersenghi, A. Javanmard, and A. Montanari. Performance of a community detection algorithm based on semidefinite programming. Journal of Physics:
Conference Series. Vol. 699. 1. IOP Publishing. 2016, p. 012015.

[TSS02]

A. Tanay, R. Sharan, and R. Shamir. Discovering statistically significant biclusters in gene expression data. Bioinformatics 18.suppl 1 (2002), S136–S144.

[Ver86]

S. Verdú. Asymptotic error probability of binary hypothesis testing for poisson
point-process observations (corresp.) IEEE Transactions on Information Theory
32.1 (1986), pp. 113–115.

[Vu14]

V. Vu. A simple SVD algorithm for finding hidden partitions. arXiv preprint
arXiv:1404.3918 (2014).

[WFL14]

J. Wyse, N. Friel, and P. Latouche. Inferring structure in bipartite networks using
the latent block model and exact ICL. 2013 (2014), p. 23. arXiv: 1404.2911.

[YP14]

S.-Y. Yun and A. Proutiere. Accurate community detection in the stochastic
block model via spectral algorithms. arXiv preprint arXiv:1412.7335 (2014).

[ZA18]

Z. Zhou and A. A. Amini. Analysis of spectral clustering algorithms for community detection: the general bipartite setting. Preprint (Mar. 12, 2018). arXiv:
1803.04547v1 [math.ST].

[ZLZ+12]

Y. Zhao, E. Levina, J. Zhu, et al. Consistency of community detection in networks
under degree-corrected stochastic block models. The Annals of Statistics 40.4
(2012), pp. 2266–2292.

52

[ZRMZ07]

T. Zhou, J. Ren, M. Medo, and Y. C. Zhang. Bipartite network projection and
personal recommendation. Physical Review E - Statistical, Nonlinear, and Soft
Matter Physics 76.4 (2007), pp. 1–7. arXiv: 0707.0540.

[ZZ+16]

A. Y. Zhang, H. H. Zhou, et al. Minimax rates of community detection in stochastic block models. The Annals of Statistics 44.5 (2016), pp. 2252–2280.

A
A.1

Remaining proofs
Proofs of Sections 6.2, 6.3 and 6.4

Proof of Lemma 4. We have nk (y (q) ) ∼ Hypergeometric(n/4, nk (y), n). For any fixed k ∈
0
[K] and q 0 ∈ [Q], the concentration of hypergeometric distribution [Chv79] gives |πk (y (q ) ) −
πk (y)| ≤ ξ with probability at least 1 − 2 exp(−nξ 2 /Q). The same probability bound holds for
|π` (z (q) ) − π` (z)| ≤ ξ, for any fixed ` ∈ [L] and q ∈ [Q]. Taking the union bound over k, `, q, q 0
gives the desired result.
Proof of Lemma 5. Recall the definition of the true local mean parameters in (25), and the
corresponding global parameters in Section 4.1. We have

λk`
n` (z) 
(q)
λk` −
= Pk` n` (z (q) ) −
Q
Q

 λ  π (z (q) )

Pk` n` (z) π` (z (q) )
k`
`
−1 =
−1 .
=
Q
π` (z)
Q
π` (z)

Since |π` (z (q) ) − π` (z)| ≤ ξ and π` (z) ≥ 1/(βL) by assumptions (B4a) and (A2), the first
inequality in (55) follows, from which we have the second inequality by (B4a).
Proof of Lemma 6. From Lemma 5, we have kΛ(q) − Λ/Qk∞ ≤ (ξLβ)kΛ/Qk∞ and kΛ(q) k∞ ≤
3
(q) is the true (local) mean parameter matrix associated with subblock
2 kΛ/Qk∞ . Note that Λ
0
0
A(q ,q) , and this subblock has n/(2Q) rows. We will apply Lemma 1 to the submatrix A(q ,q) and
0
sublabels z (q ) and y (q) . In order to do so, we have to verify conditions (A1), (A2)) and (B3)
for the subblock. (Condition (B3) is satisfied by assumption.) By Lemma 5, we have
kΛ(q) k∞
(q)
Λmin

≤3

kΛk∞
≤ 3ω.
Λmin

By (53), the condition (A2) holds with β replaced with 2β. We also need to replace Λmin in
(q)
Lemma 5 with Λmin ≥ Λmin /(2Q), and Cγ with 4Cγ (more precisely, we are replacing Cγ,β with
Cγ,2β ). Thus, assuming 6(4Cγ )(3ω) ≤ 1, we obtain




n Λmin
0
P kΛ̂(q ,q) − Λ(q) k∞ ≤ 4(4Cγ + τ )kΛ(q) k∞ ≥ 1 − 2p1 τ ;
,
, 2β
(115)
2Q 2Q
where p1 (·) is as in (43). Since kΛ(q) k∞ ≤ 23 kΛ/Qk∞ , 4(4Cγ +τ )kΛ(q) k∞ ≤ (24Cγ +6τ )kΛ/Qk∞ .
Thus, on the event in (115), we have by triangle inequality
0

kΛ̂(q ,q) − Λ/Qk∞ ≤ 4(4Cγ + τ )kΛ(q) k∞ + (ξLβ)kΛ/Qk∞
h
i
3
≤ 4(4Cγ + τ ) + ξLβ kΛ/Qk∞ ,
2
which is the desired result.
53

2×L
Proof of Lemma 7. For the proof, it is enough to consider Λ = [λ0 ; λ1 ] ∈ R+
, where λ0 , λ1 ∈
2×L
L
R+ are the two rows of Λ. Similarly, let Λ̃ = [λ̃0 ; λ̃1 ] ∈ R+ ∈ BΛ (δ). Let us define

Is (λ0 | λ1 ) =

L
X


s
(1 − s)λ0` + sλ1` − λ1−s
0` λ1` ,
`=1

λ0 , λ1 ∈ RL
+

(116)

and α` = max{|λ0` − λ̃0` |, |λ1` − λ̃1` |}. We have
L h
i
X
1−s s
s
λ
−
λ̃
Is (λ0 | λ1 ) − Is (λ̃0 | λ̃1 ) ≤
α` + λ1−s
λ̃
.
1`
1`
0`
0`
`=1

Consider the function f (a, b) = a1−s bs for a, b > 0. Assuming max{a/b, b/a} ≤ ω, we have
k∇f (a, b)k1 ≤ (1 − s)(b/a)s + s(a/b)1−s ≤ (1 − s)ω s + sω 1−s ≤ ω
using ω ≥ 1. It follows that |a1−s bsP
− u1−s v s | ≤ ω max{|a − u|, |b − v|} for a, b, u, v > 0. Thus,
Is (λ0 | λ1 ) − Is (λ̃0 | λ̃1 ) ≤ (1 + ω) ` α` ≤ 2ωLδkΛk∞ since max` α` ≤ δkΛk∞ by assumption.
Taking the supremum over s gives part (a).
Proof of Lemma 8. Assume dMis(ỹ, y) ≤ α and let nk = |i : yi = k| and Nkk0 = |i : yi =
k, ỹi = k 0 |. Then,
XX
α
n dMis(ỹ, y) =
Nkk0 ≤ αn ≤
nk =: εnk .
π
k
0
k k 6=k

P

0
It follows that k0P
6=k Nkk0 ≤ εnk for every k. We also obtain Nkk0 ≤ εnk for all k and k such
that k 6= k 0 . Since k0 Nkk0 = nk , we have Nkk ≥ (1 − ε)nk . Thus, as long as ε < 1/2, we have
Nkk > Nkk0 for all k and k 0 such that k 6= k 0 . That is, the diagonal of the confusion matrix is
bigger than every element in the corresponding row. Now take σ 6= id. Then, there exists k
such that k 0 := σ −1 (k) 6= k.
σ
:= |i : yi = k, σ(ỹi ) = k| = |i : yi = k, ỹi = k 0 | = Nkk0 < Nkk .
Nkk

Then we have
n dMis(σ(ỹ), y) =

X
k

σ
(nk − Nkk
)>

X
(nk − Nkk ) = n dMis(ỹ, y).
k

showing that id is the unique optimal permutation and proving part (a). For part (b), we note
that |{i : ỹi = k}| ≥ Nkk ≥ (1 − ε)nk > (1/2)nk whenever ε < 1/2.
Proof of Lemma 9. Assume that Mis(ỹ, y) ≤ α and Mis(ỹ 0 , y) ≤ α where α < 41 mink πk (ỹ). By
definition of the optimal permutation, dMis(σ(ỹ), y) ≤ α and dMis(σ 0 (ỹ 0 ), y) ≤ α. Since dMis
is a metric (being the sum of discrete metrics over the coordinates), we have
dMis(σ −1 ◦ σ 0 (ỹ 0 ), ỹ) = dMis(σ 0 (ỹ 0 ), σ(ỹ)) ≤ 2α <

1
min πk (ỹ)
2 k

where the first inequality is by the triangle inequality for dMis and the second by assumption.
Applying Lemma 8 gives the desired result.
54

Proof of Corollary 2. Take q = 2 for simplicity. Assume that (58) holds with constant 8 in
place of 32, which is all we need for this lemma. We have
(n/2) dMis(σ12 (ỹ (1) ), y (1) ) ≤ n dMis(σ12 (ỹ (1,2) ), y (1,2) )
by the definition of the dMis. It then follows that
dMis(σ12 (ỹ (1) ), y (1) ) < 2

1
1
≤ min πk (y (1) )
8βK
2 k

where the second inequality holds by the counterpart of (53) for row labels. Applying Lemma (8)
we conclude that σ12 = σ ∗ (ỹ (1) → y (1) ) =: σ1 .
Proof Corollary 3. Take q = 2 for simplicity. Let ε = 1/(32βK). By assumption, we have
Mis(ỹ (1,2) , y (1,2) ) < ε and Mis(ỹ (2,3) , y (2,3) ) < ε. By Corollary 2, σ12 = σ2 and σ23 = σ20 .
Then, the argument leading to (37) implies Mis(ỹ (2) , y (2) ) < 2ε and Mis(ỹ 0(2) , y (2) ) < 2ε. By
assumption,
Mis(ỹ 0(2) , y (2) ) < 2ε =

1
1
1
≤ min πk (y (2) ) ≤ min πk (ỹ (2) )
16βK
8 k
4 k

where the second inequality holds by the counterpart of (53) for row labels, and the third
inequality follows from the second inequality and Lemma 8(b). It thus follows from Lemma 9
that σ2−1 ◦ σ20 = σ ∗ (ỹ 0(2) → ỹ (2) ) which is the desired result.

A.2

Proofs of Section 7.1.1

Proof of Lemma 14. Let us define I := Ikr and
Is =

L
X
`=1

1−s s
λr` .
(1 − s)λk` + sλr` − λk`

in this proof. For s ∈ [0, 1], s 7→ Is is a concave
I0 = I1 = 0. We have defined
 function
 and
∗
1
1
s
∗
I := Is∗ = sups∈[0,1] Is . Suppose s ≥ 2 , since 0 1 − 2s∗ + 2s∗ = 21 and 2s1∗ ≥ 12 , by concavity,

1 
1
1
I
I1/2 ≥ 1 − ∗ I0 + ∗ Is∗ ≥ Is∗ =
2s
2s
2
2

Similarly, suppose s∗ ≤ 12 , I1/2 ≥ I/2 still holds, from which it follows that
X
X p
p
p
p
(λr` − λk` )2 =
( λr` − λk` )2 ( λr` + λk` )2 ≥ (I/2)(4Λmin ) = 2Λmin I.
`∈[L]

`∈[L]

Taking the minimum over k and r completes the proof.
Proof of Lemma 15. Recall our choice of ξ in (73)—which will also be assumed in this proof—
giving P(Pc ) = o(1) as shown (76). By Lemma 5, we have Λ(q) ∈ BΛ/Q (ξLβ) for all q ∈ [Q],
which combined with Lemma 7 (applied with δ = ξLβ) gives
|Ikr (Λ(q) ) − Ikr (Λ/Q)| ≤ 2ω(ξLβ)LkΛ/Qk∞ ≤
55

2ω(Lβ)LkΛ/Qk∞
2
≤ ,
βω(K ∨ L)2 (kΛk∞ ∨ kΓk∞ )
Q

using (73). Thus
Imin (Λ(q) ) ≥ Imin (Λ/Q) −

Imin
2
≥
Q
2Q

as Imin → ∞. We are now ready to apply Corollary 5 to the Algorithm 4 operating on
subblocks in Gcol
1 . It remains to verify that assumption (84) translates to condition (83) for the
subblocks. Indeed, we have to replace Imin with Imin (Λ(q) ), ω with 3ω (by Lemma 5), β with
α
n
m
col
2β (by (53)), and α with m/4
n/2 = 2 since the subblocks in G1 are of size 2 × 4 . Therefore, by
assumption (84),
6Qβ 2 ω 2 KL(K ∧ L) α
(2β)2 (3ω)KL(K ∧ L) (α/2)
≤
≤ C1 (1 + κ)−2 ,
(q)
I
2Imin (Λ )
min

(117)

verifying condition (83) on the subblocks. Applying Corollary 5, we have the misclassification
rate of ỹ (q) satisfies
Mis(ỹ (q) , y (q) ) ≤

(1 + κ)2 (3ω)(2β)L(K ∧ L)(α/2)
2C1 (Imin /2Q)

which is the desired result.

A.3

Proofs of Section 9.3.1

Proof of Proposition 3. Step 1: Interpolation. Assume without loss of generality that θ01 6=
θ11 and fix some s ∈ (0, 1). It is enough to establish the bound for ` = 1 and this particular s.
Let Pe,+ := Pe,0 + Pe,1 the be sum of the error probabilities under the two hypothesis. Then,
Z
Z
Pe,+ = p0 1{p0 ≤ p1 }dµ + p1 1{p1 < p0 }dµ
(118)
Z
Z
(119)
= min(p0 , p1 )dµ = p01−s ps1 min(ls , ls−1 )dµ.
where l = p0 /p1 is the likelihood ratio. Let pr` := π( · ; θr` ) so that pr (x) =
Similarly, let
ps := R

p1−s
ps1
0
,
p1−s
ps1 dµ
0

and ps` := R

1−s s
p0`
p1`
1−s s
p0` p1` dν

Q

` pr` (x` ).

(120)

Q
It is easy to see that ps (x) = L
`=1 ps` (x` ) and each ps` is a probability density (w.r.t. ν). One
can also verify that
Z
s
−Is`
p1−s
, and ps` = π( · ; θs` ),
0` p1` dν = e
hence ps = p( · ; θs ) using definition (98). That is, ps defined in (120) belongs to the same
s
exponential family, with parameter θs interpolating θ0 and θ1 . We also note that p1−s
0` p1` =
1−s
−I
s
−I
e s` ps` , hence p0 p1 = e s ps . Substituting into (118), we obtain
Z
−Is
Pe,+ = e
ps min(ls , ls−1 )dµ.
(121)
56

Step 2: Reduction to the single component case (L = 1). Using pr` (t) = π(t; θr` ),
we have pr` (t)/pr` (t0 ) = exp(θr` (t − t0 )), hence

Using pr (x) =
l(x) =

Q

` pr` (x` ),



p0` (t) p1` (t0 )
= exp (θ0` − θ1` )(t − t0 )
0
p0` (t ) p1` (t)

the likelihood ratio can be written as

p0 (x) Y
=
l` (x` ),
p1 (x)

where

l` (x` ) =

`



p0` (x` )
= exp (θ0` − θ1` )x − A(θ0` ) + A(θ1` ) .
p1` (x` )

(122)

As long as θ0` 6= θ1` , l` is well defined on R and maps onto R++ . For any (x2 , . . . , xL ), let
x∗1 = x∗1 (x2 , . . . , xL ) be the solution of the following equation:
l1 (x∗1 )

L
Y

l` (x` ) = 1

`=2

which always exists in R (and not necessarily on the support of the exponential family). Then,
we have, setting δ = θ01 − θ11 ,
l(x) =
It follows that



l1 (x1 )
p01 (x1 ) p11 (x∗1 )
=
= exp (θ01 − θ11 )(x1 − x∗1 ) = exp[δ(x1 − x∗1 )].
∗
∗
l1 (x1 )
p01 (x1 ) p11 (x1 )
∗

∗

min(l(x)s , l(x)s−1 ) ≤ e− min(s,1−s)|δ(x1 −x1 )| = e−α|x1 −x1 |
Q
where we have defined α := |δ| min(s, 1 − s). Recall that ps (x) = L
`=1 ps` (x` ) which we write
QL
compactly as ps = `=1 ps` . Let us write µ = µ1 ×µ2:L as the product of underlying coordinate
measures. By Fubini theorem, we first integrate over the first coordinate in (121):
Is

e Pe,+ =

Z Y
L

ps`

`=2

hZ

i
ps1 min(ls , ls−1 )dµ1 dµ2:L

(123)

Let J = J(x2 , . . . , xL ) denote the inner integral in (123) (in brackets). We have the bound
Z
Z
∗
−α|x1 −x∗1 |
1
J ≤ ps1 (x1 )e
dµ (x1 ) ≤ kps1 k∞ e−α|x1 −x1 | dµ1 (x1 ).

Note that x∗1 is the only place where dependence on (x2 , . . . , xL ) appears in the bound. Since
µ1 is either the Lebesgue or the counting measure, and both these measures are translation
invariant, the bound is in fact independent of x∗1 . That is, we have J(x2 , . . . , xL ) ≤ C(α)kps1 k∞
for all (x2 , . . . , xL ). It follows that the same bound holds for Pe,+ by (123), that is,
Is

e Pe,+ ≤ C(α)kps1 k∞
QL

Z Y
L

`=2


ps` dµ2:L = C(α)kps1 k∞

since `=2 ps` is a probability density w.r.t µ2:L . Since the choice of the coordinate ` = 1 and
s was arbitrary, the proof is complete.

57

A.4

Proofs of Section 9.3.2

P
1−s s
Proof of Lemma 17. Let f (s) = L
`=1 (s−1)λk` −sλr` +λk` λr` , then f (s) is a concave function
∗
of s on R+ . Since f (0) = f (1) = 0, s ∈ (0, 1). First, we show the statement is true when
L = 1. In this case, s∗ satisfies


λr1
1−s∗ s∗
=0
(124)
λk1 − λr1 + λk1 λr1 log
λk1
Let x =

λr1
λk1 .

Now (124) is equivalent to
∗

1 − x + xs log x = 0.
Hence
s∗ (x) =

log((x − 1)/ log x)
.
log x

We extend the domain of s∗ (x) to 1 by defining s∗ (1) = 12 , then s∗ (x) is an continuous increasing
r1
function on (0, ∞). Since λλk1
∈ [1/ω, ω], we have s∗ ∈ [s∗ (1/ω), s∗ (ω)] ⊂ (0, 1). One can
∗
∗
observe that s (x) = 1 − s (1/x), we have s∗ ∈ [s∗ (1/ω), 1 − s∗ (1/ω)]. One can also observe
1
1
that s∗ (x) ≥ x2 for x ∈ [0, 1], so s∗ ∈ [ 2ω
, 1 − 2ω
].

1−s s
λr` , we still
Now suppose L > 1, let s∗` be the optimizer of f` (s) = (s − 1)λk` − sλr` + λk`
P
1
1
∗ ∈ [ 1 , 1 − 1 ] because
have s∗ ∈ [ 2ω
f
(s)
satisfies
s
, 1 − 2ω
]. The optimizer s∗ of f (s) = L
`=1 `
2ω
2ω
f` (s) is concave for every ` ∈ [L].

Proof of Lemma 18. We first note the following Laurent series:
1
r + 1 r2 − 1
1
=
+
+
x − O(x2 ),
1 − (1 + x)−r
rx
2r
12r

as x → 0

from which we get the inequality
1
r−1 1 + r−1
≤
+
,
1 − (1 + x)−r
x
2

for x > 0, r < 1.

Let ε = ε01 and α = α01 . Applying this inequality with r = 1/(2ω) and x = ε, and recalling
1
α = 2ω
log(1 + ε), we have
C(α) ≤

2
2ω
2
=
=2
+ (1 + 2ω).
−1/(2ω)
1 − e−α
ε
1 − (1 + ε)

Using 1 ≤ ω completes the proof.
Proof of Lemma 19. We have
eλ kπ( · ; log λ)k∞ = sup

t ∈ Z+

λt
λt
eλ
≤ sup √
=√
,
t!
2πλ(t/e)t
2πλ
t ∈ R+

where the first inequality is by Stirling’s approximation and the last equality is by plugging in
the maximizer t = λ.

58

Proof of Lemma 20. There exists ` ∈ L such that
∗

∗

1−s s
(1 − s∗ )λk` + s∗ λr` − λk`
λr` ≥

Ikr
L

Without loss of generality, we assume λk` < λr` . Let s∗ be the optimizer of Ikr . Dividing λk1
both side, we have
λr`  λr` s∗
Ikr
Ikr
1 − s∗ + s∗
−
≥
≥
λk`
λk`
Lλk`
LkΛk∞
∗

Let us definef (x) := 1 − s∗ + s∗ x − xs , then for x > 1,

Thus f (x) ≥

Ikr
LkΛk∞

1
1
f (x) ≤ (1 − s∗ )s∗ (x − 1)2 ≤ (x − 1)2
2
8
q
8Ikr
implies x ≥ 1 + LkΛk
, or equivalently, log x ≥
∞

1
2


log 1 +

8Ikr
LkΛk∞



.

Proof of Lemma 21. Without loss of generality, we assume λ01 /λ11 = 1 + ε01 . Letting ρ =
δkΛk∞ , we have kΛ0 − Λk∞ ≤ ρ by definition. Let f (x) = (λ01 − x)/(λ11 + x). Then f (x) is
convex on (0, ∞) with derivative f 0 (x) = −(λ01 + λ11 )/(λ11 + x)2 , hence
λ001
λ01 − ρ
≥
= f (ρ) ≥ f (0) + ρf 0 (0)
0
λ11
λ11 + ρ
λ01 + λ11
λ01 λ01 + λ11
ρ = 1 + ε01 −
ρ.
=
−
λ11
λ211
λ211
Combined with
ρ
δkΛk∞
=
≤ ωδ
λ11
λ11

and

λ01 + λ11
= 2 + ε01 ≤ 2(1 + ε01 )
λ11

(125)

we have λ001 /λ011 ≥ 1 + ε01 − 2ω(1 + ε01 )δ which gives the desired result.

Pn
Proof of Lemma 22. Let Xj ∼ Poi(pj ) independent over j = 1, . . . , n, so that
j=1 Xj ∼
Poi(λ). Fix x ∈ Z+ and letQS(x) = Q
{S ⊂ [n] : |S| = x}. For any subset S of [n] and vectors
α, β ∈ Rn+ , let ψ(α, β, S) = j∈S αj j ∈S
/ βj . We have
X

ϕ(x; λ) = P
Xj = x
j

≥P
=

X
j


Xj = x, Xj ∈ {0, 1}, ∀j ∈ [n]

X hY

S ∈ S(x)

P(Xj = 1)

j∈S

On the other hand ϕ(x;
e p) =

P

Y

j ∈S
/

i
X

P(Xj = 0) =
ψ (pj e−pj ), (e−pj ), S .
S ∈ S(x)

S ∈ S(x) ψ((pj ), (1

− pj ), S). Thus,


P

ψ (pj ), (1 − pj ), S
ϕ(x;
e p)
S ∈ S(x) ψ (pj ), (1 − pj ), S
 ≤ max

≤P
−pj ), (e−pj ), S
ϕ(x; λ)
S∈S(x) ψ (pj e−pj ), (e−pj ), S
S ∈ S(x) ψ (pj e

= max ψ (epj ), ((1 − pj )epj ), S
S∈S(x)

59



P
P
using ( ai )/( i bi ) ≤ max(ai /bi ) which holds assuming the sums have equal number of terms,
all of which positive. Using (1 − x)ex ≤ 1, It follows that
Y

ϕ(x;
e p)
∗
≤ max ψ (epj ), (1), S = max
epj ≤ exp .
ϕ(x; λ) S∈S(x)
S∈S(x)
j∈S

Proof of Lemma 23. We have
 2ρx

2ρ x 2ρ
φ(x; λ1 )  λ1 x λ2 −λ1  λ + ρ x 2ρ 
e
≤
e = 1+
e ≤ exp
=
+ 2ρ .
φ(x; λ2 )
λ2
λ−ρ
λ−ρ
λ−ρ
Since λ − ρ ≥ 23 λ by assumption, the result follows.

P
Proof of Lemma 24. Let row node i belong to row cluster k, and let bi+ =
`∈[L] bi` =
P
Pm
j=1 Aij be its degree, with expectation λk+ :=
`∈[L] λk` . By definition, we have λk+ ≤
LkΛk∞ . We would like to find an upper bound on the probability


P bi+ > 5LkΛk∞ ≤ P bi+ − λk+ > 4LkΛk∞

We let v = λk+ , vt = 4LkΛk∞ , so t ≥ 4. By Proposition 2, we have
h 3

 3 

2t i
P bi+ − λk+ > 4LkΛk∞ ≤ exp − vt log 1 +
≤ exp − vt ≤ exp(−3LkΛk∞ ).
4
3
4

A.5

Proof of Lemma 2(a)

Proof of Lemma 2(a). Fix z̃ and let bi∗ = bi∗ (z̃). For r 6= k ∈ [K] and i such that yi = k, and
Λ̃ ∈ BΛ (δ),
Yikr (bi∗ , Λ̃) =

L h
X
`=1

bi` log

L h
i
i
X
λ̃r`
λr` + ρ
+ λ̃k` − λ̃r` ≤
bi` log
+ λk` − λr` + 2ρ := Y ∗
λk` − ρ
λ̃k`
`=1

where ρ := δkΛk∞ is the radius of BΛ (δ). Hence P(∃Λ̃ ∈ BΛ , Yikr ≥ 0) ≤ P(Y ∗ ≥ 0). By
∗
Markov inequality, we have P(Y ∗ ≥ 0) ≤ E[esY ] for any s ≥ 0. To simplify theP
notation, let us
∗
write v` = s log[(λr` + ρ)/(λk` − ρ)] and w` = s(λk` − λr` + 2ρ), so that sY = L
`=1 bi` v` + w` .
By independence, we have
L
L
L
hY
i X
X


log E[esY∗ ] = log E
ebi` v` +w` =
log E[ebi` v` +w` ] =
log ew` Eebi` v` .
`=1

`=1

`=1

Since the mgf of a Poisson-binomial variable is bounded above by that of a Poisson variable
with the same mean,
log EesYikr ≤

L
X
`=1


w` + ψ v` , λk` (y, z̃)
60

where ψ(t, µ) = µ(et − 1) is the log-mgf of a Poi(µ) random variable. Recalling the assumption
Λ(y, z̃) ∈ BΛ , we have
L
X
`=1



L h
i
X
λr` + ρ s
− λk` (y, z̃) + s(λk` − λr` + 2ρ) .
w` + ψ v` , λk` (y, z̃) =
λk` (y, z̃)
λk` − ρ


`=1

Since Λ(y, z̃) ∈ BΛ (δ), λk` (y, z̃) ≤ λk` + δkΛk∞ = λk` + ρ. Since λk` − ρ = λk` − δkΛk∞ ≥
∞
1
2
λk` − kΛk
3ω ≥ λk` − 3 Λmin = 3 λk` ,
λk`



λr` + ρ
λk` − ρ

s

r`
λr` − ρ λλk`
+ (1 +

= λk`
≤ λk`
≤ λk`
≤ λk`

"

 λ s
r`

λk`
 λ s
r`

λk`
 λ s
r`

λk`

λr`
λk` )ρ

λk` − ρ

!s

 λ s−1  (1 + λr` )ρ 
λk`
r`
+s
λk`
λk` − ρ

#

3
+ s · 2ωρ
2
+ 3ωρ.

Moreover, λr` + ρ = λr` + δkΛk∞ = λr` + ωδΛmin ≤ λr` + 31 Λmin ≤ 43 λr` . Thus, we have
ρ



λr` + ρ
λk` − ρ

s

 4 λ s
r`
≤ ρ 23
≤ 2ωρ.
3 λk`

Taking infimum over s > 0, by Lemma 17, the maximizer s∗ of Ikr is always bounded between
0 and 1, hence we have

P ∃Λ̃ ∈ BΛ , Yikr (bi∗ , Λ̃) ≥ 0 ≤ P (Y∗ ≥ 0)


≤ exp −Ikr + 8LωδkΛk∞ = exp −(1 − η 0 )Ikr .

A.6

Proofs of Section 3

Proof of Proposition 1. The upper bound has been provided by Corollary 6. Here we will show
the lower bound, using the notation established in the proof of Proposition 3 and Section 9.3.2.
We rename λk∗ and λr∗ , and work with λ0∗ and λ1∗ instead, and we assume throughout that,
λ0` , λ1` ≥ 1 for all ` ∈ [L]. We recall from (121) that
Z
Z
Pe,+ = min(p0 , p1 )dµ = e−Is ps min(ls , ls−1 )dµ,
(126)
where ps is defined in (120) and l in (122). Since µ is the counting measure, we have
Pe,+ ≥ max min(p0 (x), p1 (x)) = max e−Is ps (x) min(ls (x), ls−1 (x)).
x∈ZL
+

x∈ZL
+

61

(127)

Finding the maximizer x over Z+ gives the lower bound. First, let us extend the Poisson
density as φ(t; λ) = λt e−λ /Γ(t + 1) to any t ∈ R+ , so that l is well-defined on RL
+ , given by
l(x) = exp

X

x log

`∈[L]


λ0`
− λ0` + λ1` , x ∈ RL
+.
λ1`

P
s
Recall that λs` = λ1−s
` [(1 − s)λ0` + sλ1` − λs` ] (cf. Section 9.3.2).
0` λ1` , λs = (λs` ) and Is =
We note that
λ 
X
dIs
0`
−λ0` + λ1` + λs` log
=
. = log(l(λs ))
ds
λ1`
`∈[L]

The function s 7→ Is is concave, smooth, nonconstant (by assumption) and we have I0 = I1 = 0.
Hence, the unique maximizer s∗ of s 7→ Is belongs to (0, 1) and satisfies dIs /ds s∗ = 0, that is,
log(l(λs∗ )) = 0, or equivalently p0 (λs∗ )/p1 (λs∗ ) = l(λs∗ ) = 1. By the definition of ps , we have
s
e−Is∗ ps∗ (λs∗ ) = p1−s
0 (λs∗ )p1 (λs∗ ) = p0 (λs∗ ) = p1 (λs∗ ).

We recall that ps is the product of Poisson densities with parameters λs` . By a version of the
Stirling’s inequality for the Gamma functions [Jam15]:
Γ(x + 1) = xΓ(x) ≤ (2π)1/2 xx+1/2 e−x e1/(12x) , ∀x > 0
hence Γ(x + 1) ≤ C0 xx+1/2 e−x for all x ≥ 1, where C0 = (2π)1/2 e1/12 . Then,
φ(λ; λ) =

λλ e−λ
≥ C0−1 λ−1/2 ,
Γ(λ + 1)

from which it follows that
ps∗ (λs∗ ) =

Y

`∈L

φ(λs∗ ` ; λs∗ ` ) ≥ C0−L

Y

−1/2

λ s∗ ` .

`∈[L]

Q −1/2
Thus, e−Is∗ C0−L ` λs∗ is a lower bound on Pe,+ whenever λs∗ ∈ ZL
+ . In general, λs∗ does
L
not have integer coordinates. Instead, pick any x ∈ Z+ satisfying kx − λs∗ k`∞ ≤ 1.
Since t 7→ φ(t; λ) is a quasi-concave function (i.e., upper-level sets are convex), we have
φ(t; λ) ≥ min{φ(a; λ), φ(b; λ)} for every t ∈ [a, b], hence, for every t ∈ [a − 1, a + 1], we obtain
using Γ(x + 1) = xΓ(x),
φ(t; λ) ≥ e−λ min

n λa−1

na
λa+1 o
e−λ λa
λ o
=
min
,
,
Γ(a) Γ(a + 2)
Γ(a + 1)
λ a+1
,

that is,
na
φ(t; λ)
λ o
≥ min
,
,
φ(a; λ)
λ a+1
Since |x` − λs∗ ` | ≤ 1,
p0` (x` ) ≥ p0` (λs∗ ` ) min

nλ

s∗ `

λ0`

,

62

t ∈ [a − 1, a + 1].

λ0` o
≥ (2ω)−1 p0` (λs∗ ` )
λs∗ ` + 1

where we have used, for any s ∈ [0, 1],
min

nλ

s`

λ0`

,

n λ λ os  1 s
1
λ0` o 
1`
0`
≥
≥
= min
,
λs`
λ0` λ1`
ω
ω

and λs∗ ` /(λs∗ ` + 1) ≥ 1/2 since λs∗ ` ≥ 1. Similarly p1` (x` ) ≥ p1` (λs∗ ` )/(2ω). Hence,
min{p0 (x), p1 (x)} ≥

min{p0 (λs∗ ), p1 (λs∗ )}
e−Is∗ ps∗ (λs∗ )
e−Is∗ Y −1/2
=
≥
λ s∗
L
L
(2ω)
(2ω)
(2C0 ω)L
`∈[L]

where we have used min{p0 (x), p1 (x)} = e−Is ps (x) min{l(x)s , l(x)s−1 } and l(λs∗ ) = 1. Thus,


L
Pe,+ ≥ exp −Is∗ − L log(2C0 ω) − log kΛk∞
2


L
≥ exp −Is∗ − L log(2C0 ω 3/2 ) − log Λmin
2

using the assumption kΛk∞ ≤ ωΛmin . The proof is complete.

63

B

Extra Simulation Results

Here we present extra simulation results under the setup of Section 8. The following figure
shows the overall NMI and log. error rate for different values of C and α:
C = 0.50,

C = 0.50,

= 0.90, =1.00, K=4, L=6

= 0.90, =1.00, K=4, L=6

0.7

-0.6
0.6

log miss. (overall)

-0.8

NMI (Overall)

0.5

Spectral

0.4

Oracle
Hard
Soft

0.3

-1

-1.2

-1.4

-1.6
0.2

-1.8
0.1
200

400

600

102

800 1000 1200 1400 1600 1800 2000

# of nodes per cluster

103

# of nodes per cluster

The next figure illustrates the results for unbalanced cluster sizes. To be specific,
π(y) =

(1, 4, 6, 9)
20

and π(z) =

(1, 3, 4, 6, 7, 9)
,
30

which implies β ≥ 3 according to (A2):
C = 1.00,

C = 1.00,

= 0.75, =3.00, K=4, L=6

= 0.75, =3.00, K=4, L=6

-0.8
0.8

-1
Spectral

-1.2

Oracle

log miss. (overall)

NMI (Overall)

0.7

Hard
Soft

0.6

0.5

-1.4
-1.6
-1.8
-2

0.4

-2.2
0.3
200

400

600

102

800 1000 1200 1400 1600 1800 2000

average # of nodes per cluster

103

average # of nodes per cluster

We also consider the setting where the number of the clusters of one side is significantly
greater that of the other. We let K = 4, L = 12 and


1 2 3 4 5 6 7 8 9 10 11 12
 4 5 6 7 8 9 10 11 12 1 2 3 

B=
 7 8 9 10 11 12 1 2 3 4 5 6  .
10 11 12 1 2 3 4 5 6 7 8 9
64

The simulation results are as follows:
C = 0.50,

C = 0.50,

= 0.75, =1.00, K=4, L=12

= 0.75, =1.00, K=4, L=12

-0.55
0.9
-0.6

log miss. (overall)

NMI (Overall)

0.8

0.7

0.6

-0.65
-0.7
-0.75
-0.8
-0.85

Spectral
Oracle

0.5

-0.9

Hard
Soft

-0.95

0.4
200

400

600

102

800 1000 1200 1400 1600 1800 2000

# of nodes per cluster

103

# of nodes per cluster

65

