arXiv:1704.07114v3 [cs.IT] 20 Jun 2017

Regular Decomposition: an information and
graph theoretic approach to stochastic block
models∗
Hannu Reittu†, Fülöp Bazsó‡, Ilkka Norros§

Abstract
A method for compression of large graphs and matrices to a block
structure is proposed. Szemerédi’s regularity lemma is used as a generic
motivation of the significance of stochastic block models. Another ingredient of the method is Rissanen’s minimum description length principle
(MDL). We propose practical algorithms and provide theoretical results
on the accuracy and consistency of the method.

Keywords: Szemerédi’s Regularity Lemma, Minimum Description Length
Principle, stochastic block model, big data

Contents
1 Introduction

2

2 Basics and definitions
2.1 Szemerédi’s Regularity Lemma . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Stochastic block models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 The Minimum Description Length (MDL) principle . . . . . . . . . . . . . . . .

3
3
5
7

3 The
3.1
3.2
3.3

family of stochastic block models in the MDL
Two-part MDL for simple graphs . . . . . . . . . . .
Two-part MDL for matrices . . . . . . . . . . . . . .
Refined MDL and asymptotic model complexity . . .

4 The Regular Decomposition approach to
stochastic block models
4.1 Block model codes . . . . . . . . . . . . . . . . . .
4.2 Algorithms for description length minimization . .
4.3 Accuracy of block structure identification by MDL
4.4 ǫ-regularity vs. stochastic block models . . . . . . .
5 Proof of Theorem 4.13

.
.
.
.

framework
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .
. . . . . . . . . . . . . . .

8
8
10
11

.
.
.
.

13
13
15
20
21

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

.
.
.
.

22

∗ This

work has been submitted to the IEEE for possible publication. Copyright may be
transferred without notice, after which this version may no longer be accessible.
† VTT Technical Research Centre of Finland Ltd, P.O. Box 1000, 02044 VTT. Email
hannu.reittu@vtt.fi
‡ Department of Theory, Institute for Particle and Nuclear Physics, Wigner Research Centre
for Physics, Hungarian Academy of Sciences, P.O. Box 49, H-1525 Budapest, Hungary. Email
bazso.fulop@wigner.mta.hu.
§ VTT Technical Research Centre of Finland Ltd, P.O. Box 1000, 02044 VTT. Email
ilkka.norros@vtt.fi

1

A Chernoff bounds and other information-theoretic preliminaries

1

30

Introduction

Szemerédi’s Regularity Lemma (SRL) [24] is a fundamental result in graph
theory. Roughly speaking, SRL states that any large enough graph can be approximated arbitrarily well by nearly regular, pseudo-random bipartite graphs,
induced by a partition of the node set into a bounded number of equal-sized
sets. For many graph problems it suffices to study the problem on a corresponding random structure, resulting in a much easier problem (for a review,
see, e.g., [10]). SRL is fundamental also in theoretical computer science, for instance, in showing the existence of polynomial time approximations for solving
dense graph problems, and in characterizing the class of so-called testable graph
properties [2].
Despite the impressive theoretical applications of SRL, it has had only few
applications to ’real-life’ problems. The main reason might be that SRL has
extremely bad worst cases in the sense that the lower bound of graph sizes, for
which the partition claim holds with reasonable accuracy and without a single
exception, is enormous. Thus, a real-world application of SRL in literal sense
is impossible. However, this drawback does not mean that regular partitions
could not appear in much smaller scales relevant to applications. On the contrary, one could conjecture that regular partitions or structures be commonplace
and worth of revealing. The goal of our work is on such realistic yet preferably
large networks that nowadays appear in almost all imaginable application areas. The regular structure granted by SRL is then replaced by a probabilistic
model, substituting the regular bipartite components with truly random bipartite graphs. This model class is well known as stochastic block models [9], and
it has recently gained much attention in research and popularity in practical
applications like community detection [16, 7, 11]. However, the fundamental
nature of SRL suggests the (heuristic) conjecture that stochastic block models
in fact present a very generic form of the separation of structure and randomness in large real-world systems. These models have good algorithmic properties
in maximum likelihood fitting, and one can use powerful tools like expectation
maximization, simulated annealing and Monte Carlo Markov Chain algorithms.
There are also some examples of graph decomposition applications that have
been explicitly inspired by SRL. The practical contexts are varied: brain cortex
analysis [12], image processing [23], peer-to-peer network[15], analysing Functional Magnetic Resonance (fMRI) data to depict functional connectivity of the
brain [14], a matrix of multiple time series [18] . In the last-mentioned work,
the method was generalized from graphs to arbitrary positive matrices using a
Poissonian construction as an intermediary step.
Interestingly, the authors of the recent work [21] define a ’practical’ variant
of SRL by relaxing algorithmic SRL in a certain way to make it more usable
in machine-learning tasks, see also [17]. Bolla [3] has developed a spectral
approach for finding regular structures of graphs and matrices. Our emphasis
is more information-theoretical by nature and continues the works [12, 15, 18].
It would be very interesting to compare the methods of [21, 3] and those of this
paper in depth.
The third main ingredient in this paper is Rissanen’s Minimum Description

2

Length (MDL) principle (see [6]), according to which the quality of a model
should be measured by the length of the bit-string that it yields for unique
encoding of the data. In our case, the data has the form of a graph or matrix, and
we present the stochastic block models as a modeling space in the sense of the
MDL theory. Within this modeling space, the model corresponding to minimum
code length coding presents the optimally regular decomposition (partition) of
the data. Therefore we call this Szemerédi-motivated and technically MDLbased approach Regular Decomposition.
By information theory, the optimal coding reveals as much redundancy in
the data as is possible from the given modeling point of view. The regular
structure has a high degree of redundancy: a regular pair is an almost structureless subgraph where almost all nodes have similar and uniform connectivity
patterns. By definition, the MDL principle should be able to discover regular
structures. Note also that this principle presents a case of ’Occam’s Razor’, a
general rule of reasoning that has proved fruitful in all areas of science.
In non-hierarchic clustering tasks, it has been a major challenge to select the
’right’ size of a partition (k). Intuitively, the optimal choice of k will strike a
balance between the simplest partition of the data using a single cluster and the
maximal partition assigning each data point to its own cluster, and selecting
something in between. A popular device has been the Akaike information criterion (AIC), which was applied also in [12]. The AIC simply adds the ’number
of model parameters’ to the maximal log-likelihood and chooses the model that
minimizes the sum. Our MDL-based approach solves the corresponding model
selection task in a better founded way.
The contributions of this paper are the following: (i) the linkage of SRL,
stochastic block models and the MDL principle, (ii) the unified handling of
graphs and matrices, (iii) effective practical algorithms for revealing regular
structures in data, and (iv) Theorem 4.13 that characterizes how accurately the
MDL principle identifies a stochastic block model.
The paper is structured as follows. Section 2 presents the definitions of
the central notions: SRL, stochastic block models and MDL. The last topic
is expanded in Section 3. Section 4 presents our main technical results, in
particular the algorithms and Theorem 4.13. The proof of Theorem 4.13 is
given in Section 5. It is structured into several propositions and makes strong
use of information-theoretic tools presented in Appendix A.

2
2.1

Basics and definitions
Szemerédi’s Regularity Lemma

Consider simple graphs G(V, E), where V is the set of nodes (vertices) and E
is the set of links (edges). The link density of a non-empty node set X ⊆ V is
defined as
d(X) =

e(X)
,
|X|
2

where e(X) = {{v, w} ∈ E : v, w ∈ X} ,

3

and |·| denotes the cardinality of a set. Similarly, the link density between two
disjoint non-empty node sets X, Y ⊆ V is defined as
d(X, Y ) :=

e(X, Y )
,
|X| |Y |

where e(X, Y ) = {{v, w} ∈ E : v ∈ X, w ∈ Y } .

Definition 2.1 Let ǫ > 0. A pair of disjoint sets A, B ⊆ V is called ǫ-regular,
if for every X ⊆ A and Y ⊆ B, such that |X| > ǫ |A| and |Y | > ǫ |B| we have
|d(X, Y ) − d(A, B)| < ǫ.
A partition ξ = {V0 , V1 , V2 , · · · , Vk } of V into k + 1 sets, where all except V0
have equal cardinalities, is called ǫ-regular, iff all except at most ǫk 2 pairs are
ǫ-regular and |V0 | < ǫ|V |.
Theorem 2.2 (Szemerédi’s Regularity Lemma, [24]) For every ǫ > 0 and for
any positive integer m, there are positive integers N (ǫ, m) and M (ǫ, m), such
that for any graph G(E, V ) with |V | ≥ N (ǫ, m) there is an ǫ-regular partition of
V into k + 1 classes with m ≤ k ≤ M (ǫ, m).
SRL states, roughly speaking, that nodes of any large enough graph can
be partitioned into a bounded number, k, of equal-sized ’clusters’ and into one
small set in such a way that most pairs of clusters look like random bipartite
graphs, whose link probability equals the link density between the pair.
The claim of SRL is significant for sufficiently dense graphs, i.e., when the
link density is higher than ǫ. However, the result can be modified to apply to
sparse graphs by multiplying ǫ at the right-hand side of the regularity definition
by the link density of the entire graph [22].
The amazing fact in SRL is that the regularity claim holds for all graphs
starting from a lower bound for size that depends only on ǫ. However, it is also
well known that this dependence on ǫ is of extremely bad kind: the known lower
bound for the graph size N (ǫ, m) behaves like a tower of powers of 2:
2
..

2

2.

,

where the height of the tower is upper-bounded by 1/ǫ5 . Such a number is too
big to be considered in any applications. Thus, all real-world networks fall into
a ’gray area’ with respect to SRL.
The magnitude of N (ǫ, m) is a problem for using algorithmic versions of the
SRL, like the one introduced by Alon et al, [1]. Although its time complexity is
only polynomial O(n2.376··· ), corresponding to the time for multiplying two n×n
binary matrices, it requires this enormous size (n ≥ N (ǫ, m)) of graph to be able
to find a regular partition. A considerable improvement was found in the recent
work [5], where the running time is only linear in the graph size. Perhaps even
more importantly, this randomized algorithm works in a more realistic fashion:
for any graph, it either finds an ǫ-regular partition or concludes that such a
partition does not exists. Thus, it does not need to work in the safe region
of graph sizes specified in SRL. Another algorithm with the same feature was
suggested by Tao [25]. In principle, such algorithms could possibly be applicable
for real-world graphs, although some extra-big numbers are hidden in some of
the constants of the algorithm.
4

2.2

Stochastic block models

The notion of an ǫ-regular partition is purely combinatorial. The stochastic
model closest to this notion is the following.
Definition 2.3 Let V be a finite set and ξ = {A1 , . . . , Ak } a partition of V .
A stochastic block model is a random graph G = (V, E) with the following
structure:
• There is a symmetric k × k matrix D = (dij )ki,j=1 of real numbers dij ∈
[0, 1] satisfying the irreducibility condition that no two rows are equal, i.e.,
for all i, j, i < j, there is qij ∈ {1, . . . , k} such that diqij 6= djqij ; (1)
• For every pair {v, w} of distinct nodes of V such that v ∈ Ai , w ∈ Aj , let
evw = ewv be a Bernoulli random variable with parameter dij , assuming
that all evw ’s are independent. The edges of G are
E = {{v, w} : v, w ∈ V, v 6= w, evw = 1} .
Note that the case of the trivial partition ξ = {V } yields the classical random
graph with edge probability d11 .
A graph sequence Gn = (Vn , En ) presenting copies of the same stochastic
block model in different sizes can, for definiteness, be constructed as follows.
Pk
Construction 2.4 Let γ1 , . . . , γk be positive real numbers such that i=1 γi =
1. Divide the interval (0, 1] into k segments
#
k−1
X
γi , 1 ,
I1 = (0, γ1 ], I2 = (γ1 , γ1 + γ2 ], . . . , Ik =
i=1

and denote Γ = {I1 , . . . , Ik }. For n = 1, 2, . . ., let the vertices of Gn be


i
: i ∈ {1, . . . , n} .
Vn =
n
For each n, let ξn be the partition of Vn into the blocks
(n)

Ai

= Ii ∩ Vn ,

i = 1, . . . , k.

For small n, we may obtain several empty copies of the empty set numbered
as
blocks.
However,
from some n0 on, all blocks are non-empty and ξn =
n
o
(n)
(n)
A1 , . . . , Ak
is a genuine partition of Vn . We can then generate stochastic
block models based on (Vn , ξn , D) according to Definition 2.3.
Remark 2.5 A slightly different kind of stochastic block model can be defined
(n)
by drawing first the sizes of blocks Ai as independent Poisson(γi n) random
variables and proceeding then with the matrix D as before. The additional level
of randomness, regarding the block sizes, is however of no interest in the present
paper.
Next, we define the notion of a Poissonian block model in complete analogy
with Definition 2.3.
5

Definition 2.6 Let V be a finite set of vertices, n = |V |, and let ξ = {A1 , . . . , Ak }
be a partition of V . The symmetric Poissonian block model is a symmetric random n × n matrix E with the following structure:
• There is a symmetric k × k matrix Λ = (λij )ki,j=1 of non-negative real
numbers satisfying the irreducibility condition that no two rows are equal,
i.e.,
for all i, j, i < j, there is qij ∈ {1, . . . , k} such that λiqij 6= λjqij ; (2)
• For every unordered pair {v, w} of distinct nodes of V such that v ∈ Ai ,
w ∈ Aj , let evw = ewv be a Poisson random variable with parameter λij ,
assuming that all evw ’s are independent. The matrix elements of E are
evw for v 6= w, and evv = 0 for the diagonal elements.
P
P
Thanks to the independence assumption, the sums u∈A v∈B euv are Poisson distributed for any A, B ∈ ξ.
Remark 2.7 The rest of the technical contents of this paper focus on the simple
binary and Poissonian models of Definitions 2.3 and 2.6. However, the following
extensions are straightforward:
• bipartite graphs: this is just a subset of simple graphs;
• m × n matrices with independent Poissonian elements: a matrix can be
seen as consisting of edge weights of a bipartite graph, where the parts are
the index sets of the rows and columns of the matrix, respectively;
• directed graphs: a directed graph can be presented as a bipartite graph
consisting of two parts of equal size, presenting the input and output ports
of each node.
The following construction, inspired by our conditionally Poissonian graph
model [13], is the key to extend the MDL methodology for data that has the
most common form of a large non-negative matrix.
Construction 2.8 Let C = (cij ) be a non-negative m × n matrix. Let N be a
(rather large) integer and denote J = {1, . . . , N }. Let
(1)

(2)

(1)
V = J1 ∪ · · · ∪ J m
,

W = J1 ∪ · · · ∪ Jn(2) ,

(ι)

Ji n= {(ι,
o i)} × J. Define a bipartite Poissonian block model PN (C) = (V ∪
(·)
(·)
, C) with blocks J· and mean matrix C.
W, J·

Although we have not studied this at the technical level, it is natural to
expect that, with large N , a partition of PN (C) with minimum description
(·)
length would with high probability keep the blocks J· unbroken. Because the
regular decomposition algorithm for Poissonian block models operates only on
means over blocks (see subsection 4.2), Construction 2.8 is a strong heuristic
argument that this algorithm applies as such to regular decomposition of nonnegative matrices.

6

2.3

The Minimum Description Length (MDL) principle

The Minimum Description Length (MDL) Principle was introduced by Jorma
Rissanen, inspired by Kolmogorov’s complexity theory, and an extensive presentation can be found in Grünwald’s monography [6], see also [20]. The basic
idea is the following: a set D of data is optimally explained by a model M,
when the combined description of the (i) model and (ii) the data as interpreted
in this model is as concise as possible. By description we mean here a code that
specifies an object uniquely.
The principle is best illustrated by our actual case, simple graphs. A graph
G = (V, E) with |V | = n can always be encoded as a binary string of length
n
2 = n(n − 1)/2, where each binary variable corresponds to a node pair and a
value 1 (resp. 0) indicates an edge (resp. absense of an edge). Thus, the MDL
of G is always at most n2 . However, G may have a structure whose disclosure
would allow a much shorter description. Our heuristic postulate is that in the
case of graphs and similar objects a good a priori class of models should be
inferred from SRL, which points to stochastic block models.
Definition 2.9 Denote by Mn/k the set of irreducible stochastic block models
(V, ξ, D) with
• |V | = n,
• |ξ| = k, and, denoting ξ = {V1 , . . . , Vk },
• for i, j ∈ {1, . . . , k},
dij =

hij
, hij ∈ N,
|Vi ||Vj |

dii =

hii
,
|Vi |
2

hii ∈ N.

The condition in the last bullet entails that each modeling space Mn/k is
finite.
Remark 2.10 Without the irreducibility condition (1), there would not be a
bijection between stochastic block models and their parameterizations.
The models in Mn/k are parameterized by Θk = (ξ, D). A good model for
a graph G is the one that gives maximal probability for G and is called the
maximum likelihood model. We denote the parameter of this model
Θ̂k (G) := arg max (P (G | Θk )),

(3)

Θk ∈Mn/k

where P (G | Θk ) denotes the probability that the probabilistic model specified
by Θk produces G.
One part of likelihood optimization is trivial: when a partition ξ is selected
for a given graph G, the optimal link probabilities are the empirical link densities:
dij =

e(Vi , Vj )
, i 6= j,
|Vi ||Vj |

dii =

e(Vi )
.
|Vi |

(4)

2

Thus, the nontrivial part is to find the optimal partition for the given graph.
This is the focus of the next sections.
7

3
3.1

The family of stochastic block models in the
MDL framework
Two-part MDL for simple graphs

Let us denote the set of all simple graphs with n nodes as
Ωn = {G : G = (V, E) is a graph, |V | = n} .
A prefix (binary) coding of a finite set Ω is an injective mapping
C : Ω → ∪s≥1 {0, 1}s

(5)

such that no code is a prefix of another code. Recall the following proposition
from information theory (see, e.g., [4]):
Theorem 3.1 (Kraft’s Inequality) For an m-element alphabet there exists a
binary prefix
with code lengths l1 , l2 , · · · , lm iff the code lengths
P coding−lscheme
i
≤
1.
2
satisfy:
i=1,··· ,m

An important application of Theorem 3.1 is the following: if letters are
drawn from an alphabet with probabilities p1 , p2 , · · · , pm , then there exists a
prefix coding with code lengths ⌈− log p1 ⌉, · · · , ⌈− log pm ⌉, and such a coding
scheme is optimal in the sense that it minimizes the expected code length (in this
section, the logarithms are in base 2). In particular, any probability distribution
P on the graph space Ωn indicates that there exists a prefix coding that assigns
codes to elements of G ∈ Ωn with lengths equal to ⌈− log P ({G})⌉.
The code length l(·) is the number of binary digits in the code of the corresponding graph. In case of a large set Ω, most such codes are long and as a result
the ceiling function can be omitted, a case we assume in sequel. A good model
results in good compression, meaning that a graph can be described by much
less bits than there are elements in the adjacency matrix. An incompressible
case corresponds to the uniform distribution on Ωn and results in code length
− log (1/ | Ω |) = n2 , equivalent to writing down all elements of the adjacency
matrix.
For every graph G from Ωn and model P we can associate an encoding with
code length distribution − log P (· | Θ̂k (G)). However, this is not all, since in
order to be able to decode we must know what particular probabilistic model
P is used. This means that also Θ̂k (G) must be prefix encoded, with some
code-length L(Θ̂k (G)). We end up with the following description length:
l(G) = ⌈− log P (G | Θ̂k (G))⌉ + L(Θ̂k (G)).

(6)

Eq. (6) presents the so-called two-part MDL, [6]. In an asymptotic regime with
n → ∞, we get an analytic expression of the refined MDL. A simple way of
estimating L(Θ̂k (G)) is just to map injectively every model in Mn/k to an
integer and then encode integers with l∗ (| Mn/k (G) |) as an upper bound of the
code-length. Here
l∗ (m) = max(0, log (m)) + max(0, log log(m)) + · · · , m ∈ N,

(7)

gives, as shown by Rissanen, the shortest length prefix coding for integers (see
[6, 19]). The size of the graph must also be encoded with l∗ (n) bits (it is assumed
that there is a way of defining an upper bound of the models with given n). In
this point, it is necessary that the modeling space is finite. This results in
8

Proposition 3.2 For any graph G ∈ Ωn , there exists a prefix coding with codelength
l(G) = ⌈− log P (G | Θ̂k (G))⌉ + m
∗

≤ mk := l (n) + l

∗

S2(n, k)



!

(k2)+k
n−k+2
+1
+ 1 + 1,
2

where S2(n, k) is the Stirling number of the second kind.
Proof The expression in (6) corresponds to a concatenation of two binary
codes. The L-part is the length of a code for maximum likelihood parameters
(in the case of a non-unique maximum, we take, say, the one with smallest
number in the enumeration of all such models). The corresponding code is called
the parametric code. The parametric code uniquely encodes the maximumlikelihood model. To create such an encoding, we just enumerate all possible
models, given in Definition 2.9, and use the integer to fix the model. The length
of a prefix code corresponding to an integer is the l∗ -function computed for that
integer, and we add 1 to handle the ceiling function.
To obtain an upper bound for the parametric code length mk , we find an
upper bound for the number of models in the modeling space. The number of
models is upper-bounded by the product of the integers, each presenting the
number of partitions of an n-element set into k non-empty sets (blocks), which
equals S2(n, k). We can view the blocks of a partition as the nodes of a ‘reduced
multi-graph’ (in a multi-graph, there can be several links between a node pair,

as well as self-loops). The range of multi-links is between zero and n−k+2
:
2
if we consider a pair of blocks (or one block internally), there can be at most
n − (k − 2) nodes in such a pair (in one set, slightly less), since there must be
at least k − 2 nodes in the other blocks of the partition. Obviously,
in such a

subgraph of n − (k − 2) nodes there can be at most n−k+2
links. Thus, the
2

number of values each multi-link can take is upper-bounded by n−k+2
+ 1.
2

Since the number of node pairs in the reduced multi-graph is k2 + k, we obtain
the second multiplier in the argument of l∗ in the proposition.
Finally, we should show that the coding of the graph is prefix. We concatenate both parts into one code that has the prescribed length and put first the
prefix code of the integer that defines the parameters of the maximum likelihood model. When we start to decode from the beginning of the entire code, we
first obtain a code of an integer, because we used a prefix coding for integers.
At this stage we are able to define the probabilistic model that was used to
create the other part of the code, corresponding to the probability distribution
P (· | Θ̂k (G)). Using this information we can decode the graph G. It remains
to show that the concatenated code itself is prefix. Assume the opposite: some
prefix of such a code is prefix to some other similar code, say, the first code is
a prefix to the second one. However, the parametric code was prefix, so both
codes must correspond to the same model. Since the first two-part code is a
prefix to the second, they both share the same parametric part, and the code for
the graph of the first is a prefix of the second one. But this is impossible, since
the encoding for graphs within the same model is prefix. This contradiction
shows that the two-part coding is prefix.


9

Finally, we call
Mn :=

[

1≤k≤n

Mn/k

(8)

the full regular decomposition modeling space of Ωn .

3.2

Two-part MDL for matrices

In this section we consider input data in the form of a n × m matrix A = (aij )
with non-negative entries. With such a matrix we associate a random bipartite
multi-graph. The set of rows and the set of columns form a bipartition. Between
row i and column j there is a random number of links that are distributed
according to Poisson distribution with mean ai,j . Such a model was introduced
in [13] and it has been used in various tasks in complex network analysis, see
[8]. The aim of this model is to back up, heuristically, a corresponding practical
algorithm for regular decomposition of matrices. Our approach is closely related
to but slightly different from the Poissonian block model. Assume that A is
used to generate random n × m matrices X with independent integer-valued
elements following Poisson(aij ) distributions. The target is to find a regular
decomposition model that minimizes the expected description length of such
random matrices. We propose the following modeling spaces:
Definition 3.3 For integers k1 , k2 from ranges 1 ≤ k1 ≤ n and 1 ≤ k2 ≤ m,
the parameters of a model Θk1 ,k2 in the modeling space Mk1 ,k2 for an integer
matrix X are partition of rows into k1 non-empty sets V = (V1 , · · · Vk1 ) and
partition of columns into k2 non-empty sets U P
= (U1 , · · · , Uk2 ) and k1 × k2
xi,j
.
block average matrix P , with elements (P )α,β := i∈Vα ,j∈Uβ |Vα ||U
β|

Thanks to the addition rule of Poisson distributions, the likelihood of X in
a model Θk1 ,k2 ∈ Mk1 ,k2 , corresponds to probabilistic models where the elements of X are independent and Poisson distributed with parameters xi,j ∼
P oisson(Pα(i),β(j) ), where i ∈ Vα(i) , j ∈ Uβ(j) in the model Θk1 ,k2 . The corresponding likelihood is denoted as P (X | Θk1 ,k2 ), the actual probability of X
is denoted as P (X | A). The maximum likelihood model is found from the
program that maximizes the expected log-likelihood:
X
Θ∗k1 ,k2 = arg max
P (X | A) log P (X | Θk1 ,k2 )
Θk1 ,k2 ∈Mk1 ,k2 X

=

arg max
Θk1 ,k2 ∈Mk1 ,k2

=

arg max
Θk1 ,k2 ∈Mk1 ,k2

X
X

!
P (X|Θk1 ,k2 )
+ P (X|A) log P (X|A)
P (X|A) log
P (X|A)

(−D(PA || PΘk1 ,k2 ) − H(PA ))

where D is the Kullback-Leibler divergence between distributions, H denotes
entropy and PA and PΘk1 ,k2 are the two families of Poisson distributions for the
matrix elements of X. Since H(PA ) is independent on Θk1 ,k2 , it does not affect
the identification of the maximum likelihood model. Thus, the final program
for finding the optimal model is
Θ∗k1 ,k2 = arg min (D(PA || PΘk1 ,k2 )).
Θk1 ,k2 ∈M

10

(9)

The description length of a model l(Θk1 ,k2 ∈ Mk1 ,k2 ) consist of the description
length l(V ) + l(U ) of the two partitions and the description length of the block
average matrix l(P (X)). For the latter we need to know only the integers
presenting the block sums of X, since the denominator is known for a fixed
partition (U, V ). The code lengths of such integers are, for large matrices,
simply the logarithms of the integers. For l(U ) + l(V ) we use the same entropy
based formula as in (14). As a result we end up with the following expression for
the description length of the random multi-graph model A using the modeling
space Mk1 ,k2 :
lk1 ,k2 (A) = D(PA || PΘ∗k ,k ) + l(V ∗ ) + l(U ∗ )
1 2
X
E(log(eα,β + 1 | PΘ∗ k
+

1 ,k2

),

1≤α≤k1 ;1≤β≤k2

where

X

eα,β =

xi,j .

i∈Vα∗ ,j∈Uβ∗

The star superscript refers to parameters corresponding to the solution of the
program (9). The expectation of logarithm is not explicitly computable. However, we assume large matrices and blocks, and then Jensen’s inequality provides
a tight upper bound that can be used in practical computations. Thus, the final
expression for the description length of A is
X
lk1 ,k2 (A) = D(PA || PΘ∗k ,k ) + l(V ∗ ) + l(U ∗ ) +
log(aα,β + 1), (10)
1

2

1≤α≤k1 ;1≤β≤k2

where
aα,β =

X

ai,j .

i∈Vα∗ ,j∈Uβ∗

The full two-part MDL would now be realized by finding the global minimum
of this expression over various (k1 , k2 ). We return to this case in the algorithm
section 4.2. Although a heuristic one, we believe that our method for matrices
is both reasonable and easy to use and implement, see [18].

3.3

Refined MDL and asymptotic model complexity

Let us next consider Rissanen’s refined MDL variant (see [6]). The idea is
to generate just one distribution on Ωn , called the normalized maximum likelihood distribution Pnml . Then a graph G ∈ Ωn has the description length
− log Pnml (G) which is at most as large as the one given by the two-part code
in (6). The function P (· | Θ̂k (·)) maps
P graphs of size n into [0, 1], and it is not
a probability distribution, because G∈Ω P (G|Θ̂k (G)) > 1. However, a related
true probability distribution can be defined as
Pnml (·) = P

P (· | Θ̂k (·))

G∈Ω P (G|Θ̂k (G))

.

(11)

The problem with this is that a computation of the normalization factor
in (11) is far too involved: finding a maximum likelihood parametrization for a
11

single graph is a ‘macroscopic’ computational task by itself and it is not possible
to solve such a problem explicitly for all graphs. Therefore the two-part variant
is a more attractive choice in a practical context. However, the refined MDL
approach is useful as an idealized target object for justifying various approximations implementations of the basic idea. It appears that in an asymptotic
sense the problem is solvable for large simple graphs. The logarithm of the normalization factor in (11) is called the parametric complexity of the model space
Mn/k :
!
X
COM P (Mn/k ) := log
P (G | Θ̂k (G)) .
(12)
G∈Ωn

In a finite modelling space case like in ours, this can be considered as a definition
of model complexity. We have now the following simple bounds:
Proposition 3.4
log (S2(n, k)) ≤ COM P (Mn/k ) ≤ mk + 1,
where we use the same notation as in Proposition 3.2.
Proof The lower bound follows from the fact that we can have at least this
number of graphs that have likelihood 1 in Mn/k . This corresponds to graphs for
which the nodes can be partitioned into k non-empty sets and inside each set we
have a full graph and no links between the sets. Thus, for every partition there
is at least one graph that has likelihood one and all such graphs are different
from each other since there is a bijection between those graphs and partitions.
For the upper bound, we notice that according to Proposition 3.2, there is
a prefix coding with code lengths that
P correspond to the two-part code. As a
result, Kraft’s inequality yields that G∈Ωn 2−lk (G) ≤ 1, or
X
X
1≥
2−⌈− log P (G|Θ̂k (G))⌉−mk ≥
2log P (G|Θ̂k (G))−1−mk ,
G∈Ωn

G∈Ωn

from which we get
X

G∈Ωn

P (G | Θ̂k (G)) ≤ 2mk +1 .

Taking logarithms, we arrive at the claimed upper bound.



When considering large-scale structures corresponding to moderate k, the
upper and lower bounds in Proposition 3.4 are asymptotically equivalent, and
we have
Corollary 3.5 Assume that k > 1 is fixed. Then
COM P (Mn/k ) ∼ n log k,

n → ∞.

Proof Denoting the lower and upper bound of parametric complexity in Proposition 3.4 respectively by bl and bu , we argue that bu ∼ bl ∼ n log k asymptotically when n → ∞. This follows from the fact that the dominant asymptotic
n
component of both bu and bl is log S2(n, k). Indeed, S2(n, k) ∼ kk! for fixed k,
the asymptotic of log S2(n, k) is linear in n, and all other terms of the asymptotics of both bounds are additive and at most logarithmic in n.

12

Remark 3.6 The speed of convergence of the upper and lower bounds in Proposition 3.4 is of type log n/n.

4

The Regular Decomposition approach to
stochastic block models

4.1

Block model codes

The previous section developed both the two-part and refined variants of the
MDL theory, as presented in [6], for the model space of stochastic block models.
In the following, we formulate a variant of two-part MDL that allows both practical implementations and a proof of consistency, i.e., that the MDL principle
identifies a block model. It was shown above that the most heavy task in the
description of a block model is identifying the partition. The same is true for
the model complexity, which is asymptotically just the logarithm of the number
of partitions. It appears that in order to prove consistency, we need quite a
delicate estimate for the description length of the partition. The asymptotic
model complexity given in Corollary 3.5 seems to be too crude for consistency.
A full resolution of this intriguing question is left for further investigations.
We call our two-part MDL construction a block model code of a graph with
respect to a partition of its nodes that allows the computation of a tight upper bound of the code length. This upper bound is also consistent with a
more generic information theoretic point of view with a semi-constructive coding scheme.
We denote by H(·) both Shannon’s entropy function of a partition and the
entropy of a binomial distribution, i.e.
H(ξ) = −

X |A|
|A|
log
,
|V |
|V |

A∈ξ

H(p) = −p log p − (1 − p) log(1 − p).

Remark 4.1 In the rest of this paper, we define also information-theoretic
functions in terms of natural logarithms, and certain notions like code lengths
should be divided by log 2 to obtain their values in bits.
Definition 4.2 A block model code of a graph G = (V, E) with respect to a
partition ξ of V is a code with the following structure:
The model part:
• first, the sizes of the blocks A ∈ ξ are given as integers;
• second, the edge density d(A) inside each block A ∈ ξ and the edge density
d(A, B) between each pair of distinct blocks A, B ∈ ξ are given as the
denominators of the rational numbers presenting the exact densities.
The aim of these two codes is to describe the parameters of two probability distributions, one for the links and the other for the membership of nodes in the
blocks of the partition.
The data part:

13

• third, the partition ξ is specified by a prefix code corresponding to membership distribution P (i ∈ A) = |A|/n, where all nodes are independent of
each other;
• fourth, the edges inside each block A ∈ ξ are specified by a prefix code
corresponding to a stochastic block model distribution of links inside each
block of ξ;
• fifth, the edges between each pair of blocks A, B ∈ ξ are specified by a prefix
code corresponding to a block model distribution of links between pairs of
blocks in ξ.
The choice of link densities as link probabilities (the second code) is natural,
since conditionally to a partition the stochastic block model is just a collection
of Bernoulli models, where the best choice is to use averages as parameters.
Note that a block model code can be given for any graph with respect to any
partition of its nodes.
From Kraft’s inequality and the above definitions, it follows that there exists
a prefix code for a graph G = (V, E) with respect to a partition ξ = {A1 , . . . , Ak }
of V with length at most (and, for large graphs, typically close to)
L(G|ξ) = L1 (G|ξ) + L2 (G|ξ) + L3 (G|ξ) + L4 (G|ξ) + L5 (G|ξ),
L1 (G|ξ) =

k
X
i=1

L2 (G|ξ) =

k
X
i=1

l∗ (|Ai |),
l∗



 X
|Ai |
l∗ (|Ai ||Aj |d(Ai , Aj )) ,
d(Ai ) +
2
i<j

(13)

L3 (G|ξ) = |V |H(ξ),

k 
X
|Ai |
H(d(Ai )),
L4 (G|ξ) =
2
i=1
X
|Ai ||Aj |H(d(Ai , Aj )),
L5 (G|ξ) =
i<j

where l∗ (m) was defined by (7). Below we shall approximate l∗ (m) by log m
without further mentioning, because their difference is insignificant in our context. Similarly, we have dropped ceiling functions systematically. Also, recall
Remark 4.1 on the use of natural logarithms.
Next, we shall extend the definition of block model codes to Poissonian block
models. Because the entries of E are integers, the matrix E can be encoded
in the same way as the simple graphs in Definition 4.2. Denote by HP (a) the
entropy of the distribution Poisson(a):
Hp (a) = −

∞
X
ak
k=0

k!

e−k (k log a − log k! − a).

Now, the Poissonian block model code length of E with respect to any partition

14

η = {B1 , . . . , Bm } is defined as
L(E|η) = L1 (E|η) + L2 (E|η) + L3 (E|η) + L4 (E|η) + L5 (E|η),
m
X
log(|Bi |),
L1 (E|η) =
i=1

m
X



 X
|Bi |
log (|Bi ||Bj |a(Bi , Bj )) ,
log
a(Bi ) +
L2 (E|η) =
2
i<j
i=1

(14)

L3 (E|η) = |V |H(η),

m 
X
|Bi |
HP (a(Bi )),
L4 (E|η) =
2
i=1
X
|Bi ||Bj |HP (a(Bi , Bj )),
L5 (E|η) =
i<j

where
a(Bi ) :=

P

v,w∈Bi

evw

|Bi |(|Bi | − 1)

,

a(Bi , Bj ) :=

P

v,w∈Bi

evw

|Bi ||Bj |

.

Remark 4.3 Definition (14) parallels (13), and it works smoothly in Theorem
4.13 below. However, the expressions of L4 and L5 correspond to the case that
each set of random variables is indeed Poisson distributed with the same parameter, which is not the case when η is not a refinement of ξn . The definition
given in Section 3.2 is based strictly on the notion of description length, and it
is also applied in the practical algorithm in Section 4.2.

4.2

Algorithms for description length minimization

In this section we present algorithms that we have used in actual computations of regular decompositions of graph and matrix data. These are written
for standard two-part MDL, where the code lengths L4 and L5 have a usual
interpretation as a minus log-likelihood of a graph corresponding to a stochastic
block model.
Thus, we use link coding lengths found in the upper bound of Proposition 3.2.
In many cases, this is all that can be computed realistically. We can obviously
describe a partition into k nonempty sets using an n × k binary matrix with all
row sums equal to one and requiring that none of the column sums equals zero.
The space of all such matrices we denote as Rk and the members of this set as
R ∈ Rk .
Definition 4.4 For a given graph G ∈ Ωn with adjacency matrix A and a
partition matrix R ∈ Rk , denote
P1 (R) := RT AR,
where ·T stands for matrix transpose, the column sums of R are denoted as
nα := (RT R)α,α ,

1 ≤ α ≤ k,

and the number of links within each block and between block pairs as
1
eα,β (R) = (1 − δα,β )(P1 (R))α,β .
2
15

Then define
(P (R))α,α := 1{nα >1}

eα,β (R)
 ,
nα

(P (R))α,β :=

2

eα,β (R)
, α 6= β.
nα nβ

Then the coding length of the graph corresponding to A using the model R
is:
Definition 4.5
lk (G(A) | R ∈ Rk ) :=

X

ni nj H((P (R))i,j ) +

1≤i<j≤k

1≤i≤k

+ lk (R),
where
lk (R) =

X

X ni 
H((P (R))i,i )
2

ni H(ni /n) +

1≤i≤k

X

l∗ (ei,j (R))

1≤i≤j≤k

is the code length of the model, according to our theory and notation.
The two-part MDL program of finding the optimal model, denoted as
Rk∗ , can now be written as:
(k ∗ , Rk∗ ) := arg min min lk (G(A) | R ∈ Rk )
1≤k≤n R∈Rk

(15)

To solve this program approximately, we can use the following greedy algorithm.
Algorithm 4.6 Greedy Two-part MDL
Input: G = G(A) ∈ Ωn a simple graph of size n.
Output: (k ∗ , Rk∗ ∈ Rk∗ ), such that the two-part code for G is shortest possible
for all models in Mn by using this pair as a model.
Start: k = 1, l∗ = ∞, R ∈ Rn = {I}, k ∗ = 1, where I is denotes the n × n unit
matrix.
1. Find
R̂k (G) := arg min(lk (G | R)
R∈Rk

using subroutine ARGMAX k (Algorithm 4.8).
2. Compute lk (G) = ⌈lk (G | R̂(G))⌉ + lk (R̂(G))
3. If lk (G) < l∗ then l∗ = lk (G), Rk∗ = R̂k (G) , k ∗ = k
4. k = k + 1
5. If k > n, Print (Rk∗ , k ∗ ) and STOP the program.
6. GoTo 1.
Definition 4.7 A mapping Φ : Rk → Rk is defined as follows:
(LogP (R))α,β := log(P (R))α,β ,
(1 − LogP (R))α,β := 1 − log(P (R))α,β ,

L(R) := −AR(LogP (R))T − (1 − A)RLog(1 − P (R)),

where we set log 0 = 0 (all log 0 values will be later multiplied by 0),
β(i, R) := inf{β : β = arg min(L(R))i,α }, 1 ≤ i ≤ n,
1≤α≤k

and finally
Φ(R)i,α = δα,β(i,R) .
16

The mapping Φ(R) moves each node to a possibly different block such that
the description length would be minimized if all other nodes stay in their current
blocks.
Algorithm 4.8 ARGMAX k
Algorithm for finding optimal regular decomposition for fixed k
Input: A: the adjacency matrix of a graph (an n × n symmetric binary matrice
with zero trace); N : an integer (the number of iterations in the search of a global
optimum); k: a positive integer.
Start: m = 1.
1. i := 0; generate a uniformly random element Ri ∈ Rk .
2. If at least one of the column sums of Ri is zero, GoTo 1. Otherwise, set
Ri+1 := Φ(Ri ).
3. If Ri+1 6= Ri , set i := i + 1 and GoTo
Pn 2.
4. R(m) := Ri ; m = m + 1; l(m) := i=1 min1≤α≤k (L(R(m)))i,α .
5. If m < N , GoTo 1.
6. M := {m : l(m) ≤ l(i), i = 1, 2, ..., N }; m∗ := inf M .
Output optimal solution: R(m∗ ).
For very large graphs, the program may not be solvable in the sense that it is
not possible and reasonable to go through all possible values of k ∈ {1, 2, · · · , n}.
One option is to limit the range of k. In case that no minimum is found, then
use as an optimal choice the model found for the largest k within this range.
Another option is to find the first minimum with smallest k and stop. When
the graph is extremely large, it makes sense to use only a randomly sampled
sub-graph as an input — indeed, when k ∗ << n, a large-scale structure can be
estimated from a sample.
Our algorithm for Poissonian block models and matrices is essentially similar,
with certain differences in formulae as detailed below. Algorithms for other
cases like directed graphs and non-quadratic matrices are written very similarly,
although two partitions are needed, one for rows and one for columns. The logic
of the solution remains the same however.
A semi-heuristic two-part MDL algorithm for finding a regular decomposition for an n × m matrix A with non-negative entries works as follows. The
decomposition takes a form of a bi-clustering: there are two partitions, one for
rows and one for columns. Such partitions are described by binary matrices
with row sums equal to one. The two-part MDL program to find an optimal
regular decomposition is written as follows:
The row partition-matrices are denoted as R ∈ Rk1 with dimensions n × k1 ,
1 ≤ k1 ≤ n, and the column partition matrices as C ∈ Ck2 with dimensions
m × k2 , 1 ≤ k2 ≤ m.
Let us formulate the cost function for the matrix case that is derived from
Eq (9). The number of matrix elements in row group α and column group β
can be written as a matrix element:
(N )α,β = (RT R)α,α (C T C)β,β := nα mβ .
Assuming that all blocks are non-empty, all Nα,β > 0, we can define an average
matrix element of block α, β. First compute the sum of all matrix elements of
17

A over such a block:
eα,β = (RT AC)α,β .
The corresponding block averages form a k1 × k2 P -matrix with elements
(P )α,β =

eα,β
.
(N )α,β

The coding length of the matrix A using a two-part MDL code with partitions
(R, C) can be written as
lk1 ,k2 (G(A) | R ∈ Rk1 , C ∈ Ck2 ))
X

=
eα,β (1 − log((P )α,β ) + l∗ ([eα,β ])
1≤α≤k1 ,1≤β≤k2

+

X

nα H

1≤α≤k1

X
mβ 
nα 
+
+ k1 × k2 c.
mβ H
n
m
1≤β≤k2

Here we assume a similar handling of log 0’s as in the binary case. nα and mβ
are the sizes of row and column blocks, and [eα,β ] denotes the integer part of
eα,β ; it is assumed that such block sums are large numbers with finite decimal
precision c > 0. The description length of such decimals is the last term and it
is small compared with other terms for large matrices and can be safely ignored.
Similarly to the binary case, the two-part MDL program is defined as:
Definition 4.9
(k1∗ , k2∗ , R, C)

:= arg min
(k1 ,k2 )

min

R∈Rk1 ,C∈Ck2


lk1 ,k2 (G(A) | R ∈ Rk1 , C ∈ Ck2 )

(16)

where
1 ≤ k1 ≤ n,

1 ≤ k2 ≤ m.

The greedy algorithm of solving this program is very similar to the case of a binary matrix. The difference is that two parametric sequences of partitions must
be searched, Ri and Cj , and in the subroutine that finds the optimal partitions
for fixed i and j. One may consider different strategies in the corresponding
search of an optimal pair. For instance, moving first along the diagonal i = j,
and finding the value where the cost function (coding length of A) has a kneepoint, and after that make an off-diagonal search near that value. Another
option could be moving along the steepest descent direction of the cost function, or alternating the directions of increments in i and j, until a saturation is
reached in one direction, and then keeping that parameter fixed and finding the
optimum on the second parameter. This question is a subject to further experimenting with real and artificial data. Therefore we write only the subroutine
that finds the optimal partitioning with fixed k1 and k2 in a greedy fashion.
First we need
Definition 4.10 Define the mappings ΦR : Rk1 × Ck2 → Rk1 and ΦC : Rk1 ×
Ck2 → Ck2 as follows. Let E be an n × m matrix with all elements equal to

18

1. Then, using the definition of P -matrix and LogP -matrix (related to P as in
binary case), define two matrices using a block-matrix notation:



CP T
L(R, C) = E A
−C(LogP )T

and

M (R, C) = E T

AT

Define





RP
−RLogP



.

β1 (i, R, C) = inf{β : β = arg min(L(R, C))i,α },

1 ≤ i ≤ n,

β2 (i, R, C) = inf{β : β = arg min(M (R, C))i,α },

1 ≤ i ≤ m.

1≤α≤k1

and
1≤α≤k2

Then,
ΦR (R, C)i,α = δα,β1 (i,R,C) ,

1 ≤ α ≤ k1 ,

1 ≤ i ≤ n,

ΦC (R, C)i,α = δα,β2 (i,R,C) ,

1 ≤ α ≤ k2 ,

1 ≤ i ≤ m.

and
The main greedy subroutine is:

Algorithm 4.11 ARGMAX (k1 , k2 )
Algorithm for finding optimal regular decomposition for fixed (k1 , k2 ).
Input: A: a real n × m matrix with non-negative entries; N : a positive integer
(the number of iterations in the search of a global optimum); (k1 , k2 ): a pair of
positive integers.
Start: m = 1.
1. i = 0; generate uniformly random elements Ri ∈ Rk1 . and Ci ∈ Ck2 .
2. If at least one of column sums of Ri or Ci is zero, GoTo 1. Otherwise, set
Ri+1 := ΦR (Ri , Ci ),

Ci+1 := ΦC (Ri , Ci ).

3. If Ri+1 6= Ri or Ci+1 6= Ci , set i := i + 1 and GoTo 2.
4. R(m) := Ri ; C(m) := Ci ; m := m + 1;

1 ≤ α ≤ k1 ,
eα,β := T T (m)AC(m) α,β ,


T
T
Nα,β := R (m)R(m) α,α C (m)C β,β ,
1 ≤ α ≤ k1 ,


X
eα,β
.
eα,β 1 − log
l(m) :=
Nα,β

1 ≤ β ≤ k2 ;
1 ≤ β ≤ k2 ;

1≤α≤k1 ,1≤β≤k2

5. If m < N , GoTo 1.
6. M := {m : l(m) ≤ l(i), i = 1, 2, ..., N }; m∗ := inf M .
OUTPUT optimal solution: (R(m∗ ), C(m∗ )).
Remark 4.12 It is also possible to find a regular decompositions in the case
of partly missing matrix elements and also in the case of mixed positive and
negative entries. In the latter case, we can use the idea of directed links already
used in[12]. In the first case, we note that the main characteristic of the regular
decomposition is the P -matrix with elements that are averages of the data matrix
over large blocks that can be estimated, in many cases, despite a portion of data
is missing.
19

4.3

Accuracy of block structure identification by MDL

The general idea of the Regular Decomposition method is to be a generic tool
for separating structure and randomness in large data sets of graph or matrix
form. A partition of a real-world data set that minimizes the (nominal) code
length given in (13) resp. (14) can often not be compared with a ‘true solution’
for the simple reason that there may not be any objective notion of a ‘true
structure’ of the data. However, it is important to analyse and understand
how the method performs when the data really originates from a stochastic
block model. This question is called the consistency of MDL. Our results on
this question are summarized in the following theorem, formulated in terms of
the asymptotic behavior of a model sequence as specified by Construction 2.4.
In such a framework, an event is said to happen with high probability, if its
probability tends to 1 when n → ∞.
Theorem 4.13 Consider a sequence of stochastic block models (Gn , ξn ) based
on a vector (γ1 , . . . , γk ) of relative block sizes and a matrix D = (dij )ki,j=1 of
link probabilities, as described in Definition 2.3 and Construction 2.4. With
high probability, the following hold:
(i) Among all partitions η of Vn such that |η| ≤ k, ξn is the single minimizer
of L(Gn |η).
(ii) For any fixed ǫ ∈ (0, mini γi ), ξn is the single minimizer of L(Gn |η) among
partitions η with minimal block size higher than nǫ.
(iii) No refinement η of ξn with |η| ≤ m improves L(Gn |ξn ) by more than
const(k, m) log n.
The corresponding claims hold for the Poissonian block model mutatis mutandis.
Proof The proof is given in the next section and structured into several propositions. Claim (i) follows from Proposition 5.1, Proposition 5.2 and Corollary
5.6. Claim (ii) is Proposition 5.7. Claim (iii) is Proposition 5.5.

The results of Section 5 offer a richer picture than what was distilled into
Theorem 4.13. For example, Proposition 5.1 shows that if |η| = k and η differs
from ξn only a little, then the remaining misplaced nodes can be immediately
identified by computing their effect to the value of L(Gn |η). On the other hand,
we have not been able to exclude the possibility that a refinement of η could
yield a slight O(log n) improvement of the code length.
Remark 4.14 It is rather obvious that with large n, the identification of the
block structure ξn is robust against independent noise. The simplest case is that
the Poissonian block model is disturbed by additive Poissonian noise to each
matrix element:
ẽij = eij + φij ,
where the φij s are i.i.d. with φij ∼ Poisson(ν), ν > 0. Then (ẽij ) is again
an irreducible Poissonian block model with the same partition. More interesting
cases are binary flips in the graph case and multiplicative noise with mean 1 in
the case of non-negative matrices. We leave these for forthcoming work.

20

4.4

ǫ-regularity vs. stochastic block models

Although the structure that a MDL-based algorithm finds typically looks like
an ǫ-regular structure, there is a principal difference. In particular cases, an
ǫ-regular graph can have a structure that allows much better compression that
that provided by the ǫ-regular partition. In this section we give an explicit
example of such a case.
An important point in SRL is that for any ǫ > 0, there is an upper bound for
the size of regular partition, M (ǫ) so that for any graph with size above some
finite threshold N (ǫ), all such graphs have a regular partition with at most
M (ǫ) sets. Based on this, we show that the ǫ-regular structure of SRL and the
structure induced by the MDL need not coincide. Let us fix an order of graph
2n, large enough so that SRL holds for some ǫ > 0, and that M (ǫ) < n1−α for
some fixed 0 < α < 1/2.
Proposition 4.15 There is a graph of order 2n such that it has a MDL structure with code length o(n2 ) and an ǫ−regular structure that allows only Θ(n2 )
code length.
Proof Take n large enough as prescribed above. Then construct a bipartite
graph (X, Y ) with |X| = |Y | = n such that n is divisible by nα with some
rational α ∈ (0, 1/2). Assume that bothP
parts of the
Pbipartition are further
partitioned into equal size blocks: X =
Xi , Y =
Yi , |Xi | = |Yi | = nα .
Define then a random graph Gp = (X, Y, E) as follows. For each pair (Xi , Yj ),
take e(Xi , Yj ) = ξi,j |Xi |Yj |, where ξi,j ∼ Ber(p) is a Bernoulli random variable
with parameter 0 < p < 1, and the variables for different pairs are independent.
Assume that there are no other edges.
We show that, with high probability, Gp is ǫ-regular with regular partition
{X,
such a structure has a coding length at least
 Y }. In the M DL-approach,
n
2
H(p
+
o(1))
=
O(n
).
This
comes
from the log-likelihood part, and the o(1)
2
corresponds to very small deviations of link densities from the expected value
p that can be made arbitrarily small by increasing n. Now we check that the
ǫ-regularity of graph Gp has a positive probability, which implies that such an ǫregular pair exists (actually, it appears that this happens with high probability).
ǫ-regularity means that for any X ′ ⊆ X, Y ′ ⊆ Y , | X ′ |, | Y ′ |> ǫn, the link
density d(X ′ , Y ′ ) deviates from the link density of the pair, d(X, Y ), no more
than by ǫ. By definition,
P
′
′
i,j ξi,j | X ∩ Xi || Y ∩ Yj |
′
′
d(X , Y ) =
,
| X ′ || Y ′ |
and as a result the expectation is
P
Ed(X ′ , Y ′ ) =

Denote
xi,j :=

i,j

Eξi,j | X ′ ∩ Xi || Y ′ ∩ Yj |
| X ′ || Y ′ |

ξi,j | X ′ ∩ Xi || Y ′ ∩ Yj |
| X ′ || Y ′ |

= p.

The range
P of xi,j is interval [0, 1] of unit length. Hoeffding’s inequality yields
for S = xi,j that
2

P (| S − ES |> t) ≤ 2e
21

− Pt 2

b
i,j

,

|X ′ ∩X ||Y ′ ∩Y |

i
j
where bi,j =
is the range of variable xi,j . The denominator of the
|X ′ ||Y ′ |
exponent in the right hand-side of the Hoeffding inequality can be bounded as
 n 2
X
1
1
1 X
n2α n2α
= 2 n2α .
| X ′ ∩ Xi |2 | Y ′ ∩ Yj |2 ≤
b2i,j ≤
2
2
(ǫn)
(ǫn)
nα
ǫ

By taking t = ǫ | X ′ || Y ′ |≥ ǫ(ǫn)2 we get for the link density:
P (| d(X ′ , Y ′ ) − p |> ǫ) ≤ 2e−2ǫ

8

n2(1−α)

.

Finally, since there are at most 4n pairs of subsets (X ′ , Y ′ ) from which to choose,
the probability that none of them violates regularity is lower bounded by
1 − 2 × 4n e−2ǫ

8

n2(1−α)

→ 1,

if the exponent has a positive power of n, and this happens when α < 1/2.
Thus, all large subsets have densities that deviate from expectation less than ǫ
with a probability tending to one. Thus, we have shown the ǫ-regularity of the
partition (X, Y ).
On the other hand, using MDL, we could reach the level of small sets Xi
and Yi , and the corresponding log-likelihood is zero. The model complexity is
o(n2 ), as can be easily seen from asymptotic formulas for the upper bound for
mk , with k = n1−α .


5

Proof of Theorem 4.13

Through this section, we consider a sequence (Gn |ξn ) of increasing versions of a
fixed stochastic block model based on a vector (γ1 , . . . , γk ) of relative block sizes
and a matrix D = (dij )ki,j=1 of link probabilities, as specified in Construction
2.4.
Consider partitions η of Vn . Denote
d(η, ξn ) =

1
max min |B \ A|.
n B∈η A∈ξn

Thus, d(η, ξn ) = 0 if and only if η is a refinement of ξn . If v ∈ Vn and B ∈ η,
denote by ηv,B the partition obtained from η by moving node v to block B (if
v ∈ B, then ηv,B = η).
Proposition 5.1 There is a number ǫ0 > 0 such that the following holds with
high probability: if |η| = k and d(η, ξn ) ≤ ǫ0 , then
if A ∈ ξn , B ∈ η,

1
|B \ A| ≤ ǫ0 and v ∈ A \ B, then L(Gn |ηv,B ) < L(Gn |η).
n

Proof Let ǫ, δ > 0 be small numbers and m a positive integer to be specified.
They can be chosen so that the following holds:
• ǫ is small so that η and ξn nearly overlap when d(η, ξn ) ≤ ǫ:
mnǫ ≤ δ min |A|;
A∈ξn

22

(17)

• all the differing link probabilities are widely separated in δ units:
δ≤

1
min {|dij1 − dij2 | : i, j1 , j2 ∈ {1, . . . , k} , dij1 6= dij2 } ;
m

(18)

• the empirical densities are close to their mean values: for any Ai , Aj ∈ ξn
(possibly i = j), we have, with high probability,
|d(Ai , Aj ) − dij | + mǫ ≤ δ.

(19)

Let η be a partition of Vn such that d(η, ξn ) ≤ ǫ. Condition (17) entails that for
each block Ai ∈ ξn there is a unique block Bi ∈ η such that |Bi \ Ai | ≤ ǫ. Let
us now assume that v ∈ Ai ∩ Bj and i 6= j, and compare the partitions η and
ηv,Bi . Denote bi = |Bi |, i = 1, . . . , k, and B̃i = Bi ∪ {v}, B̃j = Bj \ {v}. Then
L4 (Gn |η) + L5 (Gn |η) − (L4 (Gn |ηv,Bi ) + L5 (Gn |ηv,Bi ))



 

 
bi + 1
bj − 1
bj
bi
H(d(Bj )) −
H(d(B̃i )) −
H(d(Bi )) +
H(d(B̃j ))
=
2
2
2
2
Xh
+
bi bq H(d(Bi , Bq )) + bj bq H(d(Bj , Bq ))
q6=i,j

i
− (bi + 1)bq H(d(B̃i , Bq )) + (bj − 1)bq H(d(B̃j , Bq ))

+ bi bj H(d(Bi , Bj )) − (bi + 1)(bj − 1)H(d(B̃i , B̃j )).

(20)

Consider first the sum over q. Leaving out the common factor bq , each term of
the sum can be written as


bj − 1
H(d(B̃j , Bq ))
bj H(d(Bj , Bq )) −
bj


bi
H(d(Bi , Bq ))
− (bi + 1) H(d(B̃i , Bq )) −
bj + 1
 

1
bj − 1
d(B̃j , Bq ) + d({v} , Bq )
= bj H
bj
bj

bj − 1
1
−
H(d(B̃j , Bq )) − H(d({v} , Bq ))
bj
bj
 

1
bi
− (bi + 1) H
d(Bi , Bq ) +
d({v} , Bq )
bi + 1
bi + 1

1
bi
H(d(Bi , Bq )) −
H(d({v} , Bq ))
−
bi + 1
bi + 1
(note the addition and subtraction of the term H(d({v} , Bq ))). Using Lemma
A.7 and the assumptions on ǫ, δ and m, the last expression can be set to be,
with high probability, arbitrarily close to the number
I(djq : diq ) + I(diq : diq ) = I(djq : diq )
(the function I(·, ·), the Kullback-Leibler divergence of Bernoulli distributions,
is defined in (32)). Thus, the sum over q is, with high probability, close to
X
bq
I(djq , diq ).
q6=i,j

23

Let us then turn to the remaining parts of (20) that refer to two codings of
the internal links of Bi ∪ Bj . Similarly as above, we can add and subtract terms
to transform these parts into
!

 "
bj −1
bj
b
−
1
j
2
H
 d(B̃j ) + bj  d({v} , B̃j )
bj
2
2
2
#

bj −1
bj − 1
2
− bj  H(d(B̃j )) − bj  H(d({v} , B̃j ))
2

2


"
bi + 1
−
H
2

bi
−

"

bi
2
 d(Bi )
bi +1
2



2
 H(d(Bi ))
bi +1
2

+ bi bj H



−

+

bi

bi

!

 d({v} , Bi )

bi +1
2

#

 H(d({v} , Bi ))

bi +1
2


1
bj − 1
d(Bi , B̃j ) + d({v} , Bi )
bj
bj

#
1
bj − 1
H(d(Bi , B̃j )) − H(d({v} , Bi ))
−
bj
bj
" 

bi
1
− (bi + 1)(bj − 1) H
d(Bi , B̃j ) +
d({v} , B̃j )
bi + 1
bi + 1
#
bi
1
−
H(d(Bi , B̃j )) −
H(d({v} , B̃j ))
bi + 1
bi + 1
≈ (bj − 1)I(dij : djj ) − (bi + 1)I(dii : dii ) + bi I(dii : dij ) − (bj − 1)I(dij : dij )

= bj I(dij : djj ) + bi I(dii : dij ).

By the above analysis of (20), we have obtained
L4 (Gn |η) + L5 (Gn |η) − (L4 (Gn |ηv,Bi ) + L5 (Gn |ηv,Bi ))
X
≈ bq
I(djq : diq ) + bj I(dij : djj ) + bi I(dii : dij ).

(21)

q6=i,j

By the irreducibility assumption (1), there is a block Aq such that dqi 6= dqj ,
with the possibility that q ∈ {i, j}. It follows that at least one of the I(x : y)’s
in (21) is positive. Denote
κ∗ = min {I(dij1 : dij2 ) : dij1 6= dij2 } .
Thus, with high probability,
L4 (Gn |η) + L5 (Gn |η) − (L4 (Gn |ηv,Bi ) + L5 (Gn |ηv,Bi )) >

1 ∗
(κ min γi )n.
i
2

On the other hand, it is easy to compute that
L3 (Gn |η) − L3 (Gn |ηv,Bi ) = n(H(η) − H(ηv,Bi )) → log

γj
.
γi

The changes of L1 and L2 when moving from η to ηv,Bi are negligible. This
concludes the proof.

24

The proof of Proposition 5.1 showed that when d(η, ξn ) ≤ ǫ0 , moving any
node to its correct block decreases L(Gn |η) at least by (1/2)(κ∗ mini γi )n. In
particular, with high probability, ξn is the unique minimizer of L(Gn |η) among
k-partitions η satisfying d(η, ξn ) ≤ ǫ0 .
Proposition 5.2 For any ǫ ∈ (0, 1) and positive integer m, there is a constant
cǫ such that the following holds with high probability:
if |η| ≤ m and d(η, ξn ) > ǫ, then

1
(L(Gn |η) − L(Gn |ξn ∨ η)) ≥ cǫ .
n2

Proof Fix an ǫ ∈ (0, 1) and let η be a partition of Vn such that d(η, ξn ) > ǫ.
By the concavity of H, we have
L4 (Gn |η) + L5 (Gn |η)
X |B|
1
=
H(d(B)) +
2
2
B∈η

X

|B||B ′ |H(d(B, B ′ ))

B,B ′ ∈η
B6=B ′

X X |A ∩ B|
H(d(A ∩ B))
≥
2
B∈η A∈ξn

+

X1
2

B∈η

X

|A ∩ B||A′ ∩ B|H(d(A ∩ B, A′ ∩ B))

(22)

A,A′ ∈ξn
A6=A′

+

1
2

X

B,B ′ ∈η
B6=B

X

A,A′ ∈ξn

|A ∩ B||A′ ∩ B ′ |H(d(A ∩ B, A′ ∩ B ′ ))

′

= L4 (Gn |η ∨ ξn ) + L5 (Gn |η ∨ ξn ).
By assumption, there is B ∈ η such that |B \ A| > ǫn for every A ∈ ξn . It is
easy to see that there must be (at least) two distinct blocks, say Ai and Aj ,
such that
ǫ
n.
(23)
min {|Ai ∩ B|, |Aj ∩ B|} ≥
k−1

By the irreducibility assumption (1), there is a block Aq such that dqi 6= dqj ,
with the possibility that q ∈ {i, j}. Fix an arbitrary δ > 0 to be specified later.
By ǫ-regularity (claim (ii) of Lemma A.7), with high probability, every choice
of a partition η with |B \ A| > ǫn results in some blocks Ai , Aj , Aq with the
above characteristics plus the regularity properties
|d(Ai ∩ B, Aq ∩ B ′ ) − diq | ≤ δ,

|d(Aj ∩ B, Aq ∩ B ′ ) − djq | ≤ δ,

(24)

where B ′ denotes a block of η that maximizes |Aq ∩ B ′ | (note that because

25

|η| ≤ m, |Aq ∩ B ′ | ≥ |Aq |/m). By the concavity of H,
|Ai ∩ B||Aq ∩ B ′ |H(d(Ai ∩ B, Aq ∩ B ′ ))

+ |Aj ∩ B||Aq ∩ B ′ |H(d(Aj ∩ B, Aq ∩ B ′ ))

|Ai ∩ B|
H(d(Ai ∩ B, Aq ∩ B ′ ))
= |(Ai ∪ Aj ) ∩ B||Aq ∩ B ′ |
|(Ai ∪ Aj ) ∩ B|

|Aj ∩ B|
+
H(d(Aj ∩ B, Aq ∩ B ′ ))
|(Ai ∪ Aj ) ∩ B|
< |(Ai ∪ Aj ) ∩ B||Aq ∩ B ′ |H(d((Ai ∪ Aj ) ∩ B, Aq ∩ B ′ )).

In the case that q ∈ {i, j} and B = B ′ , we obtain a similar equation where
|Aq ∩ B| is partly replaced by |Aq ∩ B| − 1. Because of (24) and (23), the
difference between the sides of the equality has a positive lower bound that
holds with high probability. On the other hand, this difference is part of the
overall concavity inequality (22).

Proposition 5.3 For any refinement η of ξn , we have
(st)

L4 (Gn |ξn ) + L5 (Gn |ξn ) − (L4 (Gn |η) + L5 (Gn |η)) ≤

M(|η|)−M(k)

X

(log 2 + Yi ),

j=1

(25)

where (st) refers to stochastic order, the Yi ’s are i.i.d. Exp(1) random variables,
and
x(x + 1)
.
(26)
M (x) =
2
Proof Here we apply results presented in Appendix A. Denote by η ∩ Ai the
subset of η whose members are subsets of the block Ai of ξn . Writing the edge
code lengths of the coarser and finer partition similarly as in (22), taking the
difference and using (35), we obtain
L4 (Gn |ξn ) + L5 (Gn |ξn ) − (L4 (Gn |η) + L5 (Gn |η))
k
X |B|
X
X
1
I(d(B) : dii ) +
=
2
2
i=1
B∈η∩Ai

|B||B ′ |I(d(B, B ′ ) : dii )

B,B ′ ∈η∩Ai
B6=B ′

!


|Ai |
−
I(d(Ai ) : dii )
2
+

X
i<j

X

X

B∈η∩Ai B ′ ∈η∩Aj

′

′

!

|B||B |I(d(B, B ) : dij ) − |Ai ||Aj |I(d(Ai , Aj ) : dij ) .
(27)

Applying Proposition A.6 to each term of both outer sums now yields the claim,
because
k
X
i=1

(M (|η ∩ Ai |) − 1) +

X
i<j

(|η ∩ Ai ||η ∩ Aj | − 1) = M (|η|) − M (k).

26

Remark 5.4 It is rather surprising that the stochastic bound (25) depends only
on the number of blocks in η — not on their relative sizes, nor on the overall
model size n!
Proposition 5.5 For any positive integer m > k, the following holds with high
probability:
L(Gn |ξn ) −

min

η≥ξn , |η|≤m

L(Gn |η) ≤ (m + M (m − k)) log n,

where the relation η ≥ ξn means that η is a refinement of ξn , and M (·) was
defined in (26).
Proof Let η be a refinement of ξn . Refining the partition w.r.t. ξn yields a
gain, based on the concavity of H, in the code part L4 + L5 , but costs in the
parts L1 , L2 and L3 . We have to relate these to each other.
Consider first the value of L2 (Gn |η). In our analysis, it is important to
distinguish between ‘large’ and ‘tiny’ blocks, where the relative sizes of large
blocks exceed some pre-defined number ǫ and the rest can be arbitrarily small,
even singletons. Now, each block Ai ∈ ξn must contain at least one block
Bi ∈ η ∩ Ai such that |Bi | ≥ |Ai |/m. Define
ǫ := min |Aj |/(2mn).
j

Because no concavity gain can be obtained with an index pair {i, j} such that
dij ∈ {0, 1}, it does not restrict generality to assume that dij ∈ (0, 1) for all i, j.
Then
!
  
k
X
X
X
1
|B|
′
log(|B||B |dii )
L2 (Gn |η) =
log
dii +
2
2
i=1
B∈η∩Ai

B,B ′ ∈η∩Ai
B6=B ′

+

X
i<j

≥

k
X
i=1

+

X
i<j

X

B∈η∩Ai

log



X

B ′ ∈η∩A

2 2

ǫ n
dii
2



j

!

log(|B||B ′ |dij )

!

+ (|η ∩ Ai | − 1) log(ǫndii )



log ǫ2 n2 dij + (|η ∩ Ai | + |η ∩ Aj | − 2) log(ǫndij )

= (k|η| + k) log n + c1 (D, ǫ).

On the other hand, the bound |Ai | ≤ n yields
L2 (Gn |ξn ) ≤ (k 2 + k) log n + c2 (D).
Thus,
L2 (Gn |η) − L2 (Gn |ξn ) ≥ k(|η| − k) + c1 (D, ǫ) − c2 (D).

(28)

We obviously have also L1 (Gn |η) ≥ L1 (Gn |ξn ), but this difference is insignificant
in the present context.

27

The refinement gain in code parts L4 and L5 was bounded in Proposition 5.3
stochastically by Exp(1) random variables. The rate function (see the beginning
of Appendix A) of the distribution Exp(1) is
IE (x) = x − 1 − log x.
Denoting
Gainn (η) := L4 (Gn |ξn ) + L5 (Gn |ξn ) − (L4 (Gn |η) + L5 (Gn |η)),
Proposition 5.3 yields, for y > log 2,
P (Gainn (η) > y)


M(|η|)−M(k)
X
Yi > y − (M (|η|) − M (k)) log 2
≤P
j=1




y − (M (|η|) − M (k)) log 2
≤ exp −(M (|η|) − M (k))IE
M (|η|) − M (k)

(M(|η|)−M(k))

2e
≤ exp − y + (M (|η|) − M (k)) log y
,
M (|η|) − M (k)

where the second factor is bounded and will be henceforth neglected.
For two refinements of ξn , write η ′ ∼ η if the block sizes of η ′ in each Ai
are identical to those of η. The number of refinements η ′ of ξn with η ′ ∼ η is
upperbounded by


X
exp 
|A|H(η|ξn ) = enH(η|ξn ) .
A∈ξn

On the other hand, we have

L3 (Gn |η) − L3 (Gn |ξn ) = nH(η) − nH(ξn ) = nH(η|ξn ).
Write
∆L123 (η) = L1 (Gn |η) + L2 (Gn |η) + L3 (Gn |η)
− (L1 (Gn |ξn ) + L2 (Gn |ξn ) + L3 (Gn |ξn )),
∆L(η) = L(Gn |η|) − L(Gn |ξn |).

Denote z(η) := (|η| + M (|η| − k)) log n. Recalling (28), the union bound yields


P sup Gainn (η ′ ) > ∆L123 (η) + z(η)
η ′ ∼η

≤ exp (nH(η|ξn ) − ∆L123 (η) − z(η) + (M (|η|) − M (k)) log(∆L123 (η) + z(η)))



≤ exp − z(η) + − k(|η| − k) + M (|η|) − M (k) log n + c3 (D, ǫ) .

The number of different block size sequences ℓ1 ≥ · · · ≥ ℓ|η| is upper bounded

28

by n|η|−1 = e(|η|−1) log n . Thus, a second application of the union bound yields
!
P

sup

∆L(η) > z(η)

η≥ξn , |η|≤m


 
 
 

≤ me − k(m−k) + M(m)−M(k) − m+M(m−k) + m−1 log n+const
const
→ 0, as n → ∞.
=
n

Corollary 5.6 ξn is the unique minimizer of L(Gn |η) among partitions η with
|η| = k.
Proof
By Proposition 5.1,
min

|η|=k, d(η,ξn )≤ǫ0

L(Gn |η) > L(Gn |ξn )

with high probability. On the other hand, Proposition 5.2 yields that, with high
probability,
min

|η|=k, d(η,ξn )>ǫ0

L(Gn |η) > L(Gn |ξn ∨ η) + cǫ0 n2 .

By Proposition 5.5,
min L(Gn |ξn ∨ η) > L(Gn |ξn ) − (k + M (k 2 − k)) log n

|η|=k

with high probability. It remains to note that n2 grows faster than log n.



Proposition 5.7 Let ǫ ∈ (0, mini γi ). Consider refinements η of ξn with relative minimal block size ǫ, i.e. the set
Bǫ(n) = {η ≥ ξn : |B| ≥ nǫ ∀B ∈ η} .

(29)

With high probability,
min
(n)

η∈Bǫ

\{ξn }

L(Gn |η) > L(Gn |ξn ).

(n)

(30)

Proof The restriction η ∈ Bǫ implies |η| ≤ ⌊1/ǫ⌋ =: m. The difference to
the conditions of Proposition 5.5 is now just the magnitude of L2 (η). When
η ∈ B (n) ,
L2 (Gn |η) ≥ M (|η|) log n2 + c(D, ǫ, |η|),
so that
L2 (Gn |η) − L2 (Gn |ξn ) ≥ 2(M (|η|) − M (k)) log n + c(D, ǫ, |η|).

29

By a corresponding computation as in the proof of Proposition 5.5, we obtain
(n)
for any fixed η ∈ Bǫ that


1
′
P sup Gainn (η ) > ∆L123 (η) − log n
2
η ′ ∼η
s
≤ exp (−(M (|η|) − M (k)) log n + const ) .
Proceeding with a second union bound like in the proof of Proposition 5.5 yields
!
1
P
min
L(Gn |η) ≤ L(Gn |ξn ) + log n
(n)
2
η∈Bǫ \{ξn }

 



1
≤ max exp
− (M (q) − M (k)) + q − 1 +
log n + const
q≤m
2


1
= exp − log n + const
2
const
= √ → 0, as n → ∞,
n
where the maximum over q was obtained with the smallest value q = k + 1. 

A

Chernoff bounds and other information-theoretic
preliminaries

Consider a random variable X with moment generating function
φX (β) = E eβX ,
and denote DX = {β : φX (β) < ∞}. We restrict to distributions of X for which
DX is an open (finite or infinite) interval. The corresponding rate function is
IX (x) = − inf (log φX (β) − βx).
β∈DX

IX (x) is a strictly convex function with minimum at E X and value +∞ outside
the range of X. As a further restriction on the distribution of X, we assume
that IX (x) approaches smoothly +∞ P
both when x decreases or increases from
n
the point E X. For the mean X̄n = n1 1 Xi of i.i.d. copies of X, we have
IX̄n (x) = nIX (x).

The Chernoff bound (also known as Cramér-Lundberg bound):
Proposition A.1
P (X < x) ≤ e−IX (x)

P (X > x) ≥ e

−IX (x)

30

for x < E X,
for x > E X.

(31)

This has the following simple consequence. Let the support of X be the
closure of (x− , x+ ) and write
−
+
IX (x) = IX
(x)1{(x− ,E X]} (x) + IX
(x)1{[E X,x+ )} (x) + ∞ · 1{R\(x− ,x+ )} (x).
−
+
With the assumptions made above, the functions IX
(x) and IX
(x) are, respec−
+
tively, bijections from (x , E X] and [E X, x ) to [0, ∞).

Lemma A.2 For any z ≥ 0 we have
(st)

IX (X) ≤ log 2 + Y,
(st)

where ≤ denotes stochastic order and Y is a random variable with distribution
Exp(1).
Proof

P (IX (X) > z) = P 1{X<0} IX (X) > z or 1{X>0} IX (X) > z


= P 1{X<0} IX (X) > z + P 1{X>0} IX (X) > z




+(−1)
−(−1)
(z)
(z) + P X > IX
= P X < IX
−(−1)

≤ exp(−IX (IX
=e

−z

=e

−(z−log 2)+

+e

+(−1)

(z))) + exp(−IX (IX

(z)))

−z

,

where the inequality comes from Proposition A.1.



In the case that X has the Bernoulli(p) distribution, we have
IX (x) = I(x : p) := q log

1−q
q
+ (1 − q) log
.
p
1−p

(32)

Lemma A.3 The first and second derivatives of the functions H(x) and I(x :
p) are
H ′ (x) = log

1−x
,
x

1
,
x(1 − x)
1
I ′′ (x : p) =
.
x(1 − x)
H ′′ (x) = −

I ′ (x : p) = H ′ (p) − H ′ (x),

(33)
(34)

In particular,
H(q) − (H(p) + H ′ (p)(q − p)) = I(q : p),
and
 


1
1
1
1
lim n H (1 − )y + x − (1 − )H(y) + H(x)
n→∞
n
n
n
n


= (x − y)H ′ (y) − (H(x) − H(y))
= I(x : y).

31

(35)

Proposition A.4 Let n ≥ 2 and let X1 and X2 be independent random variables with distributions Bin(m, p) and Bin(n−m, p), respectively. Denote X12 =
X1 + X2 and X̄1 = X1 /m, X̄2 = X1 /(n − m), X̄12 = X12 /n. Then the following
identities hold:
mI(X̄1 : p) + (n − m)I(X̄2 : p) − nI(X̄12 : p)




m − X1 m
X1 m
+ (n − Xn )I
:
:
= Xn I
X12 n
n − X12 n


X12 − X1
= mI(X̄1 : X̄12 ) + (n − m)I
: X̄12 .
n−m

(36)
(37)
(38)

The identities in Proposition (A.4) are obtained by writing the full expression
of (36) and re-arranging the log terms in two other ways. Formulae (37) and (38)
are written without X2 , expressing the fact that any two of the three random
variables X1 , X2 and X12 contain same information as the full triple. Note
that (37) and (38) do not contain p. This reflects the fact that the conditional
distribution of X1 given X12 , known as the hypergeometric distribution, does
not depend on p. The identity of (36) and (38) can be interpreted so that the
two positive terms of (36) measure exactly same amount of information about
p as what is subtracted by the negative term. Moreover, (38) has the additional
interpretation of presenting the rate function of the hypergeometric distribution:
Proposition A.5 Let X have the distribution Hypergeometric(n, m, z), i.e. the
conditional distribution of X1 of Proposition A.4 given that X12 = z. The rate
function of X is


x z
z−x z
+ (n − m)I
.
(39)
:
:
IX (x) = mI
m n
n−m n
Proof

Define the bivariate moment-generating function of (X1 , X2 )
φ(α, β) = E eαX1 +βX2 .

Write
P[ X1 = m | X12 = z] =

P (X1 = m, X2 = z − m)
,
P (X12 = z)

and note that we can assume p = z/n. We can now derive the claim using
φ(α, β) in similar manner as in the well-known proof of the one-dimensional
Chernoff bound.

Proposition A.6 Let k ≥ 2 and let Xi , i ∈ {1, . . . , k}, be independent random
P
variables with distributions Bin(ni , p), respectively. Denote X1...j = ji=1 Xi ,
Pj
X̄i = Xi /ni and X̄1...j = X1...j / i=1 ni . Then
k
X
i=1

X
(st) k−1

ni I(X̄i : p) − nI(X̄12 : p) ≤

(log 2 + Yi ) ,

where Y1 , . . . , Yk−1 are independent Exp(1) random variables.
32

(40)

i=1

(41)

Proof

For k = 2, the left hand side of (40) equals


X12 − X1
: X̄12
n1 I(X̄1 : X̄12 ) + n2 I
n2

(42)

by Proposition A.4. For any N ∈ {0, . . . , n}, consider the conditional distribution of (42), given that X12 = N . By Proposition A.5, this is the distribution
of the Hypergeometric(n, n1 , N ) rate function taken at the random variable X1
with the same distribution. The claim now follows by Lemma A.2, because the
stochastic upper bound does not depend on N , i.e. on the value of X12 .
For k > 2 we proceed by induction. Assume that the claim holds for k − 1
and write
k
X
i=1

=

ni I(X̄i : p) − nI(X̄12 : p)

k−1
X
i=1

ni I(X̄i : p) − (n − nk )I(X̄1...(k−1) : p)

+ nk I(X̄k : p) + (n − nk )I(X̄1...(k−1) : p) − nI(X̄1...k : p).
By the induction hypothesis,
the first row of the second expression is stochastiPk−2
cally bounded by i=1 (log 2 + Yi ), irrespective of the value of X1...(k−1) . Similarly, the second row is stochastically bounded by log 2+Yk , where Yk ∼ Exp(1),
irrespective of the value of X1...k . It remains to note that Yk can be chosen to
be independent of (Y1 , . . . , Yk−2 ), because Xk is independent of (X1 , . . . , Xk−1 ),
and of X̄1...(k−1) in particular.

Lemma A.7 Consider the sequence (Gn , ξn ) of stochastic block models as in
Theorem 4.13. Then the following holds.
(i) For any blocks Ai and Aj such that dij 6∈ {0, 1}, it holds for an arbitrary
ǫ > 0 with high probability that
min

v∈Ai

e({v} , Aj )
1
≥ dij − n− 2 +ǫ ,
|Aj |

max

v∈Ai

e({v} , Aj )
1
≤ dij + n− 2 +ǫ .
|Aj |

(ii) For any ǫ > 0, the partition ξn is ǫ-regular with high probability.
Proof

Claim (i): By Proposition A.1 and (34),


e({v} , Aj )
P max
> dij + h
v∈Ai
|Aj |


X
e({v} , Aj )
> dij + h
≤
P
|Aj |
v∈Ai

≤ |Ai | exp (−|Aj |I(dij + h : dij ))



h2
h3 ′′′
= |Ai | exp −|Aj |
+ I (z : dij )
.
2dij (1 − dij )
6

33

1

The last expression converges to zero with the choice h = n− 2 +ǫ (recall that
|Ai | ∼ nγi and |Aj | ∼ nγj ), which proves the claim on the maximum. The case
of the minimum is symmetric.
Claim (ii): Fix ǫ > 0 and consider any i, j. Let U1 ⊆ Ai and U2 ⊆ Aj such
that |U1 | ≥ ǫ|Ai | and |U2 | ≥ ǫ|Aj |. By Proposition A.1,
P (|d(U1 , U2 ) − dij | > ǫ) ≤ e−|U1 ||U2 |I(dij +ǫ:dij ) + e−|U1 ||U2 |I(dij −ǫ:dij ) .
Let ι(ǫ) = min {I(dij + ǫ : dij ), I(dij − ǫ : dij )}. The union bound yields
P (∃U1 ⊆ Ai , U2 ⊆ Aj : |U1 | ≥ ǫ|Ai |, |U2 | ≥ ǫ|Aj |, |d(U1 , U2 ) − dij | > ǫ)

≤ 2|Ai ||Aj | exp (|Ai | + |Aj |) log 2 − ǫ2 |Ai ||Aj |ι(ǫ)

 
(γi + γj ) log 2
− γi γj ǫ2 ι(ǫ)
→ 0 as n → ∞,
≤ 2n2 exp n2
n
because ι(ǫ) > 0.



Remark A.8 Because ι(ǫ) ∝ ǫ2 as ǫ → 0, the proof of claim (ii) indicates that
with a fixed ǫ, ǫ-regularity starts to hold when n >> ǫ−4 .
Poissonian block model
Denote by HP (a) the entropy of the distribution Poisson(a)
Hp (a) = −

∞
X
ak
k=0

k!

e−k (k log a − log k! − a)

= a − a log a + e−a

∞
X
ak
k=0

k!

log k!,

(43)
(44)

and by IP (b : a) the Kullback-Leibler divergence between distributions Poisson(a)
and Poisson(b)
IP (b : a) = a − b + b log b − b log a.
(45)
We note the following facts:

Lemma A.9 The entropy function HP (α) is increasing and concave, and its
two first derivatives have the expressions
d
X +1
HP (α) = E log
,
dα
α 

1
1
d2
,
HP (α) = − + E log 1 +
dα2
α
X +1

(46)
(47)

where X denotes a random variable with distribution Poisson(α).
Lemma A.10 For any α ∈ [0, 1] and x, y > 0, denote z = αx + (1 − α)y. Then
we have
HP (z) − (αHP (x) + (1 − α)HP (y)) ≤ αIP (x : z) + (1 − α)IP (y : z)
αx
= zIBer (α) ( )
z
= IBin (z,α) (αx),
where the form (50) is applicable when z is a positive integer.
34

(48)
(49)
(50)

Proof

By (44) and (45),

HP (z) − (αHP (x) + (1 − α)HP (y)) − (αIP (x : z) + (1 − α)IP (y : z))
= f (z) − (αf (x) + (1 − α)f (y)) ≤ 0,

because the function
f (u) = e−u

∞
X
uk
k=2

is convex:
f ′′ (u) = e−u

k!

log k!

∞
X
uk
k=1

k+2
log
+ log 2
k!
k+1

!

> 0.

The expression (50) is obtained by writing the right hand side of (48) with the
substitution y = (z − αx)/(1 − α) and re-combining the log terms.

The Poissonian counterpart of Proposition A.6 is the following.
P
Proposition A.11 Let a > 0, k ≥ 2, ni ≥ 1, i = 1, . . . , k, and n = i ni .
Let Xi , i ∈ {1, . . . , k}, be independent random variables with distributions
Pj
Poisson(ni a), respectively. Denote X1...j = i=1 Xi , X̄i = Xi /ni and X̄1...j =
Pj
X1...j / i=1 ni . Then
nHP (X̄12 ) −

k
X
i=1

X
(st) k−1

ni HP (X̄i ) ≤

(log 2 + Yi ) ,

(51)

i=1

(52)

where Y1 , . . . , Yk−1 are independent Exp(1) random variables.
Proof The proof of Proposition A.6 can be imitated, using the following observations:
• Using induction, it suffices to consider the case k = 2.
• By (48), the left hand side of (40) can be upperbounded in terms of
Kullback informations.
• The latter expression can be transformed to the form (50), i.e. the value
of a the rate function of a binomial distribution taken at random variable
X1 .
• Now, the conditional distribution of X1 given X12 is the above binomial
distribution. Thus, we can apply Lemma A.2 in a similar way as in the
proof of Proposition A.6.


Acknowledgment. This work was partly supported by Academy of Finland project 294763 (Stomograph).

35

References
[1] N. Alon, R. Duke, H. Lefmann, V. Rödl, and R. Yuster. The algorithmic
aspects of the regularity lemma. Journal of Algorithms, 16:80–109, 1994.
[2] N. Alon, E. Fischer, I. Newman, and A. Shapira. A combinatorial characterization of the testable graph properties: It’s all about regularity. In
STOC’06, Seattle, U.S., May 2006.
[3] M. Bolla. Spectral Clustering and Biclustering. Wiley, 2013.
[4] T. Cover and J. Thomas. Elements of Information Theory. John Wiley
and Sons, U.S.A., 1991.
[5] E. Fischer, A. Matsliah, and A. Shapira. Approximate hypergraph partitioning and applications. Foundations of Computer Science, pages 579–589,
2007.
[6] P. Grünwald. Minimum Description Length Principle. The MIT Press,
2007.
[7] S. Heimlicher, M. Lelarge, and L. Massoulié. Community detection in the
labelled stochastic block model. In NIPS Workshop on Algorithmic and
Statistical Approaches for Large Social Networks, 2012.
[8] R. v. d. Hofstad. Random graphs and complex networks. Cambridge University Press, 2016.
[9] P. Holland, K. Laskey, and S. Leinhardt. Stochastic blockmodels: First
steps. Social Networks, 5(2):109–137, 1983.
[10] J. Komlós and M. Simonovits. Szemerédi’s regularity lemma and its applications in graph theory. In D. Miklós, V. Sós, and T. Szonyi, editors,
Combinatorics, Paul Erdös is Eighty, pages 295–352. János Bolyai Mathematical Society, Budapest, 1996.
[11] L. Massoulié. Community detection thresholds and the weak Ramanujan
property. In Proceedings of the Forty-sixth Annual ACM Symposium on
Theory of Computing, STOC ’14, pages 694–703, New York, NY, USA,
2014. ACM.
[12] T. Nepusz, L. Négyessy, G. Tusnády, and F. Bazsó. Reconstructing cortial networks: case of directed graphs with high level of reciprocity. In
B. Bollobás, R. Kozma, and D. Miklós, editors, Handbook of Large-Scale
Random Networks, number 18 in Bolyai Society of Mathematical Studies,
pages 325–368. Springer, 2008.
[13] I. Norros and H. Reittu. On a conditionally Poissonian graph process. Adv.
Appl. Prob., 38:59–75, 2006.
[14] I. Pappas and P.Pardalos.
Identifying cognitive states using regularity partitions.
PLoS ONE, 10(8):
e0137012, 2015.
https://doi.org/10.1371/journal.pone.0137012.

36

[15] V. Pehkonen and H. Reittu. Szemerédi-type clustering of peer-to-peer
streaming system. In Proc. Cnet 2011, San Francisco, U.S.A., 2011.
[16] T. P. Peixoto. Entropy of stochastic blockmodel ensembles. Physical Review
E, 85(056122), 2012.
[17] M. Pelillo, I. Elezi, and M. Fiorucci. Revealing structure in large graphs:
Szemerédis regularity lemma and its use in pattern recognition. Pattern
Recognition Letters, 87:4–11, 2017.
[18] H. Reittu, F. Bazsó, and R. Weiss. Regular decomposition of multivariate
time series and other matrices. In P. Fränti, G. Brown, M. Loog, F. Escolano, and M. Pelillo, editors, Proc. S+SSPR 2014, number 8621 in LNCS,
pages 424–433. Springer-Verlag, 2014.
[19] J. Rissanen. A universal prior for integers and estimation by minimum
description length. The Annals of Statistics, 11(2):416–431, 1983.
[20] J. Rissanen. Stochastic Complexity in Statistical Inquiry, volume 15. World
scientific, 1998.
[21] G. Sárkozy, F. Song, E. Szemerédi, and S. Trivedi. A practical regularity partitioning algorithm and its application in clustering. 2012.
arXiv:1209.6540v1 [math.CO ] 28 Sep 2012.
[22] A. Scott. Szemerédi’s regularity lemma for matrices and sparce graphs.
Combinatorics, Probability and Computing, 20(3):455–466, 2011.
[23] A. Sperotto and M. Pelillo. Szemerédi’s regularity lemma and its applications to pairwise clustering and segmentation. In Proc. EMMCVPR 2007,
Italy, 2007.
[24] E. Szemerédi. Regular partitions of graphs. In Problemés Combinatoires
et Théorie des Graphes, number 260 in Colloq. Intern. C.N.R.S., pages
399–401, Orsay, 1976.

[25] T.
Tao.
Szemerédi’s
regularity
lemma
via
random
partitions,
2009.
Blog
entry,
https://terrytao.wordpress.com/2009/04/26/szemeredis-regularity-lemma-via-random-parti

37

