# Detecting Mesoscale Structures by Surprise

Emiliano Marchese  IMT School for Advanced Studies, Piazza S.Francesco 19,
55100 Lucca (Italy)  Guido Caldarelli  ‘Ca’ Foscari’ University of Venice,
Dorsoduro 3246, 30123 Venice (Italy)  Tiziano Squartini  IMT School for
Advanced Studies, Piazza S.Francesco 19, 55100 Lucca (Italy)  Institute for
Advanced Study (IAS), University of Amsterdam, Oude Turfmarkt 145, 1012 GC
Amsterdam (The Netherlands)

###### Abstract

The importance of identifying the presence of mesoscale structures in complex
networks can be hardly overestimated. So far, much attention has been devoted
to the detection of _communities_ , _bipartite_ and _core-periphery_
structures on _binary_ networks: such an effort has led to the definition of a
unified framework based upon the score function called _surprise_ , i.e. a
p-value that can be assigned to any given partition of nodes, on both
_undirected_ and _directed_ networks. Here, we aim at making a step further,
by extending the entire framework to the _weighted_ case: after reviewing the
application of the surprise-based formalism to the detection of binary
mesoscale structures, we present a suitable generalization of it for detecting
weighted mesoscale structures, a topic that has received much less attention.
To this aim, we analyze four variants of the surprise; from a technical point
of view, this amounts at employing four variants of the hypergeometric
distribution: the _binomial_ one for the detection of binary communities, the
_multinomial_ one for the detection of binary ‘bimodular’ structures and their
_negative_ counterparts for the detection of communities and ‘bimodular’
structures on weighted networks. On top of that, we define two ‘enhanced’
variants of surprise, able to encode both binary and weighted constraints and
whose definition rests upon two suitable generalizations of the hypergeometric
distribution itself. As a result, we present a general, statistically-grounded
approach to detect mesoscale structures on networks via a unified, suprise-
based framework. To illustrate the performance of our methods, we, first, test
them on a variety of well-established, synthetic benchmarks and, then, apply
them to several real-world networks, including social, economic, financial and
ecological ones. Moreover, we attach to the paper a Python code implementing
all variants of the surprise considered in the present work.

complex networks and mesoscale structures detection and surprise minimization

###### pacs:

89.75.Fb; 02.50.Tt; 89.65.Gh

## INTRODUCTION

The importance of identifying the signature of some kind of mesoscopic
organization in complex networks (be it due to the presence of communities or
bipartite, core-periphery, bow-tie structures) can be hardly overestimated
Fortunato2016 ; Bisma2017 , the best example of complex systems whose behavior
is deeply affected by their mesoscopic structural organization (e.g.
resilience to the propagation of shocks, to the failure of nodes, etc.) being
provided by financial networks Borgatti2000 ; Craig2010 ; IntVeld2014 ;
Luu2018 .

So far, much attention has been devoted to the detection of _binary mesoscale
structures_ , i.e. _communities_ and, to a far less extent, _core-periphery
structures_ : the efforts to solve these problems have led to a number of
approaches that are briefly sketched below (for a detailed review of them, see
Fortunato2010 ; Cimini2021 ).

Community detection has been initially approached by attempting a definition
of ‘communities’ based on the concepts of _clustering_ , _cliques_ , _k-core_
; core-periphery structures have, instead, been defined in a purely top-down
fashion, by imagining a fully connected subgraph (i.e. the core) surrounded by
(peripherical) vertices exclusively linked to the first ones Borgatti2000 . As
stressed in IntVeld2014 , the deterministic character of these definitions
makes their _tout court_ application to real-world systems extremely
difficult. This is the reason why the intuitive requirements that ‘the number
of internal edges is larger than the number of external edges’ and that ‘the
core portion of a network is densely connected, while its periphery is loosely
connected’ Cimini2021  are, now, interpreted in a purely probabilistic way: a
community has, thus, become a subgraph whose vertices have a larger
probability to be inter-connected than to be connected to any other vertex in
the graph - and analogously for the core-periphery structure. In other words,
the top-down approach defining a golden standard and looking for (deviations
from) it has left the place to a bottom-up approach where structures are
supposed to emerge as the result of non-trivial (i.e. non-casual) interactions
between the nodes.

This change of perspective leads to a number of problems. The first one
concerns the definition of models stating how edges are formed and has been
solved by adopting the rich formalism defining the Exponential Random Graphs
framework Fronczak2014 ; the second, and most important, one concerns the
definition of a (statistically-sound) procedure for selecting the best model
among the ones providing competing descriptions of the data. This has led to
the identification of three broad classes of algorithms111For a different
classification of the algorithms proposed for the detection of mesoscale
structures see Pesciotto2021 .: according to our intuition, all of them
implement some kind of statistical inference, the major difference lying in
the way the corresponding test of hypothesis is implemented; from a practical
point of view, instead, all these methods are designed for _optimization_ ,
the functional form of the specific score function determining the class to
which a given algorithm belongs.

The most representative algorithms among those belonging to the first class
are the ones based on modularity. Although these methods compare the empirical
network with a benchmark, they do not provide any indication of the
_statistical significance_ of the recovered partition, the reason being that
none of them is designed as a proper statistical test. For instance, let us
consider the definition of modularity, whose generic addendum is proportional
to the term $(a_{ij}-p_{ij})$: while it embodies a comparison between the
(empirical) adjacency matrix $\mathbf{A}$ and the matrix of probability
coefficients $\mathbf{P}$ defining the benchmark, it does not implement any
proper test of hypothesis. The algorithms prescribing to maximize a plain
likelihood function belong to this group as well. Although popular, this way
of proceeding is known to be affected by _overfitting_ issues: an example is
provided by the straight maximization of the likelihood defining the
Stochastic Block Model (SBM), or its degree-corrected version, over the whole
set of possible partitions, outputting the trivial one where each vertex is a
cluster on its own Karrer2011 .

The aforementioned, major limitation is overcome by the algorithms belonging
to the second class. They implement tests of hypothesis either _‘a la Fisher_
or _à la Neyman-Pearson_ , i.e. either defining a single benchmark (the _null
hypothesis_ of the first scenario) or two, alternative ones (the _null_ and
the _alternative hypothesis_ of the second scenario): from a practical point
of view, such a result is achieved by identifying the aforementioned
benchmarks with proper probability distributions and the best partition of
nodes with the one minimizing the corresponding p-value. Surprise-based
algorithms belong to this second group: as it has been shown in DeJeude2019
\- for a particular case; such a result will be generalized in what follows -
optimizing (asymptotic) surprise amounts at carrying out a (sort of)
Likelihood Ratio Test aimed at choosing between two alternative models.

Hypothesis testing can be further refined by allowing for more than two
hypotheses to be tested at a time: results of the kind are particularly useful
for _model selection_ and, in fact, have produced a plethora of criteria (e.g.
the Akaike Information Criterion, the Bayesian Information Criterion and the
Minimum Description Length) for singling out the best statistical model out of
a basket of competing ones. Generally speaking, optimization, here, seeks for
the maximum of a ‘corrected’ likelihood function embodying the trade-off
between accuracy and parsimony of a description. An example of the algorithms
belonging to this third class is represented by Infomap; other examples are
represented by the recipes employing the SBM within a Bayesian framework (see
Pesciotto2017  and the references therein).

With this paper, we pose ourselves within the second research line and adopt a
bottom-up approach that prescribes to compare any empirical network structure
with the outcome of a properly-defined benchmark model. To this aim, we devise
a unified framework for mesoscale structures detection based upon the score
function called _surprise_ 222Notice that our function is named _surprise_
since it generalizes the function proposed in Aldecoa2013  for community
detection. However, in our case it indicates a proper probability and not its
logarithm, as in Aldecoa2013 . The same holds true for its asymptotic
expression., i.e. a p-value that can be assigned to any given partition of
nodes, on both _undirected_ and _directed_ networks: while for _binary
community detection_ this is achieved by employing the _binomial
hypergeometric distribution_ Aldecoa2013 ; Nicolini2016 , with ‘bimodular’
structures like the bipartite and the core-periphery ones, one needs to
consider its _multinomial_ variant DeJeude2019 . Here, however, we aim at
making a step further, by extending the entire framework to the _weighted_
case. As a result, we present a general, statistically-grounded approach to
the problem of detecting mesoscale structures on networks via a unified,
suprise-based framework.

Before moving forward, we would like to stress that the use of the
hypergeometric distribution to carry out tests of hypothesis on networks is
not novel: examples are provided by the papers Tumminello2011  (where the
authors introduce a method to provide a statistically-validated, monopartite
projection of a bipartite network - the considered null hypothesis encoding
the heterogeneity of the system), Bongiorno2017  (where the authors employ the
same validation procedure to detect cores of communities within each set of
nodes of a bipartite system) and Musciotto2021  (where the authors extend the
framework proposed in the aforementioned references to carry out a statistical
validation of motifs observed in hypergraphs). For a review on the use of the
hypergeometric distribution for network analyses see Micciche2019  and the
references therein.

## MESOSCALE STRUCTURES DETECTION
VIA EXACT TESTS

Surprise has recently received a lot of attention: the advantages of employing
such a score function have been extensively discussed in Aldecoa2013 ;
Nicolini2016 ; DeJeude2019 ; Jiang2014 ; Delser2016 ; Tang2019 ; Gamermann2020
where researchers have tested and compared its performance from a purely
numerical perspective. However, a characterization of the statistical
properties of surprise is still missing: for this reason, we will, first, make
an effort to ‘translate’ the problem of detecting a given mesoscale network
structure into a proper _exact significance test_ (to the best of our
knowledge, the only other attempt of the kind - specifically, to detect
communities - is the one in Sadamori2018 ) and, then, show how the rich - yet,
still underexplored - surprise-based formalism can properly answer such a
question.

The basic equation underlying exact tests reads

$\displaystyle\text{Pr}(x\geq x^{*})=\sum_{x\geq x^{*}}f(x)$ (1)

and returns the probability of observing an outcome of the random variable $X$
which is ‘more extreme’ than the realized one, i.e. $x^{*}$. In the setting
above, $f$ represents the distribution encoding the _null hypothesis_ and
$\text{Pr}(x\geq x^{*})$ \- commonly known with the name of _p-value_ \-
answers the question _is the realized value $X=x^{*}$ compatible with the
(null) hypothesis that $X$ is distributed according to $f$?_

The p-value constitutes the basic quantity for carrying out any _significance
test_ : hence, in what follows we will tackle the problem of detecting the
signature of _statistically-significant_ mesoscale organizations by
individuating specific tests (or, equivalently, suitable functional forms for
$f$).

## DETECTION OF MESOSCALE STRUCTURES
IN BINARY NETWORKS

Community detection. The topic of nodes partitioning into densely connected
groups has traditionally received a lot of attention, with applications
ranging from the analysis of social networks to the definition of novel
recommendation systems Fortunato2016 ; Bisma2017 . Within the surprise-based
framework, binary community detection is carried out via the identification

$\displaystyle f(l_{\bullet})$ $\displaystyle\equiv$
$\displaystyle\text{H}(l_{\bullet}|V,V_{\bullet},L)$ (2) $\displaystyle=$
$\displaystyle\frac{\prod_{i=\bullet,\circ}\binom{V_{i}}{l_{i}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}$

i.e. by calculating the p-value

$\mathscr{S}\equiv\sum_{l_{\bullet}\geq l_{\bullet}^{*}}f(l_{\bullet})$ (3)

of a _binomial hypergeometric distribution_ whose parameters read as above. In
this formalism, the $\bullet$ subscript will be meant to indicate quantities
that are _internal_ to communities, while the $\circ$ subscript will be meant
to indicate quantities that are _external_ to communities. More precisely, the
binomial coefficient $\binom{V_{\bullet}}{l_{\bullet}}$ enumerates the number
of ways $l_{\bullet}$ links can be redistributed _within_ communities, i.e.
over the available $V_{\bullet}$ node pairs, while the binomial coefficient
$\binom{V_{\circ}}{l_{\circ}}$ enumerates the number of ways the remaining
$l_{\circ}=L-l_{\bullet}$ links can be redistributed _between_ communities,
i.e. over the remaining $V_{\circ}=V-V_{\bullet}$ node pairs. Notice that,
although $l_{\bullet}$ is ‘naturally’ bounded by the value $V_{\bullet}$, it
cannot exceed $L$ \- whence the usual requirement
$l_{\bullet}\in[l_{\bullet}^{*},\min\\{L,V_{\bullet}\\}]$.

From a merely statistical point of view, surprise ‘considers’ a network as a
population of $V$ node pairs, $L$ of which have been drawn; out of the $L$
extracted ones, $l_{\bullet}$ node pairs have the desired feature of being
‘internal’ to communities, since they connect some of the node pairs belonging
to the set $V_{\bullet}$. Hence, for a given partition of nodes into
communities, $\mathscr{S}$ quantifies the probability of observing at least
$l_{\bullet}^{*}$ ‘successes’ (i.e. intra-cluster edges) out of $L$ draws: the
lower this probability, the ‘more surprising’ the observation of the
corresponding partition, hence the ‘better’ the partition itself.

Asymptotic results. We can gain more insight into the surprise-based formalism
above upon deriving an asymptotic expression for $\mathscr{S}$ Traag2015 . To
this aim, let us consider that it can be simplified upon Stirling-
approximating the binomial coefficients that appear within it. By exploiting
the recipe $n!\simeq\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}$, $\mathscr{S}$
can be rewritten as

$\mathscr{S}\simeq\sum_{l_{\bullet}\geq
l^{*}_{\bullet}}A(l_{\bullet})\left[\frac{\text{Ber}(V,L,p)}{\prod_{i=\bullet,\circ}\text{Ber}(V_{i},l_{i},p_{i})}\right]$
(4)

where the expression

$\text{Ber}(x,y,z)=z^{y}(1-z)^{x-y}$ (5)

defines a _Bernoulli_ probability mass function, the parameters appearing in
eq. (4) read $p=\frac{L}{V}$ and $p_{i}=\frac{l_{i}}{V_{i}}$ and the
coefficient in front of the sum is

$A(l_{\bullet})=\sqrt{\frac{\sigma^{2}}{2\pi\prod_{i=\bullet,\circ}\sigma^{2}_{i}}}$
(6)

with $\sigma^{2}=Vp(1-p)$ and $\sigma^{2}_{i}=V_{i}p_{i}(1-p_{i})$. Equation
(4) makes it explicit that employing $\mathscr{S}$ for binary community
detection ultimately amounts at comparing the description of a networked
configuration provided by the Random Graph Model (RGM), and encoded into the
expression

$\text{Ber}(V,L,p)=p^{L}(1-p)^{V-L},$ (7)

with the description of the same configuration provided by the Stochastic
Block Model (SBM) Karrer2010 , and encoded into the expression

$\prod_{i=\bullet,\circ}\text{Ber}(V_{i},l_{i},p_{i})=p_{\bullet}^{l_{\bullet}}(1-p_{\bullet})^{V_{\bullet}-l_{\bullet}}\cdot
p_{\circ}^{l_{\circ}}(1-p_{\circ})^{V_{\circ}-l_{\circ}}$ (8)

where $p_{\bullet}=\frac{l_{\bullet}}{V_{\bullet}}$ and
$p_{\circ}=\frac{l_{\circ}}{V_{\circ}}$. Naturally, the SBM takes as input the
two sets indexed by $\bullet$ and $\circ$ and distinguish the connections
found ‘within the clusters’ - contributing to the probability of the whole
configuration with the term
$\text{Ber}(V_{\bullet},l_{\bullet},p_{\bullet})=p_{\bullet}^{l_{\bullet}}(1-p_{\bullet})^{V_{\bullet}-l_{\bullet}}$
\- from the ones found ‘between the clusters’ - contributing to the
probability of the whole configuration with the term
$\text{Ber}(V_{\circ},l_{\circ},p_{\circ})=p_{\circ}^{l_{\circ}}(1-p_{\circ})^{V_{\circ}-l_{\circ}}$.

Notice also that the asymptotic expression of the surprise guarantees that the
parameters of the null models defining it are tuned according to the maximum-
of-the-likelihood principle. To see this explicitly, let us consider
$\text{Ber}(x,y,z)$ whose log-likelihood reads

$\mathscr{L}=y\ln z+(x-y)\ln(1-z);$ (9)

upon maximizing it with respect to $z$, one finds $z=\frac{y}{x}$.

The way $\mathscr{S}$ works, i.e. by comparing two different null models, is
reminiscent of more traditional likelihood ratio tests, where a null
hypothesis $H_{0}$ (in our case, a given partition is compatible with the RGM)
is tested against an alternative hypothesis $H_{1}$ (in our case, the given
partition is compatible with the SBM): as the asymptotic expression of the
surprise clarifies, minimizing it amounts at finding the partition least
likely to occour under the RGM than under the SBM.

‘Bimodular’ structures detection. The surprise-based framework can be easily
extended to detect what can be called ‘bimodular structures’, a term that will
be used to compactly indicate core-periphery Zhang2014 ; Barucca2015 ;
Kojaku2018  and bipartite structures Holme2003 ; Estrada2005 . The reason for
adopting such a terminology lies in the evidence that both kinds of structures
are defined by bimodular partitions, i.e. partitions of nodes into two
different groups.

As shown elsewhere DeJeude2019 , the issue of detecting binary ‘bimodular’
structures can be addressed by considering a _multivariate_ (or _multinomial_)
_hypergeometric distribution_ , i.e. by identifying

$\displaystyle f(l_{\bullet},l_{\circ})$ $\displaystyle\equiv$
$\displaystyle\text{MH}(l_{\bullet},l_{\circ}|V,V_{\bullet},V_{\circ},L)$ (10)
$\displaystyle=$
$\displaystyle\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}}{l_{i}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V_{\top}}{l_{\top}}}{\binom{V}{L}}$
$\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V-(V_{\bullet}+V_{\circ})}{L-(l_{\bullet}+l_{\circ})}}{\binom{V}{L}}$

where $V_{\top}\equiv V-(V_{\bullet}+V_{\circ})$ indicates the number of node
pairs between the modules $\bullet$ and $\circ$ and $l_{\top}\equiv
L-(l_{\bullet}+l_{\circ})$ indicates the number of links that must be assigned
therein. While the binomial coefficient $\binom{V_{\bullet}}{l_{\bullet}}$
enumerates the number of ways $l_{\bullet}$ links can redistributed _within_
the first module (e.g. the core portion) and the binomial coefficient
$\binom{V_{\circ}}{l_{\circ}}$ enumerates the number of ways $l_{\circ}$ links
can redistributed _within_ the second module (e.g. the periphery portion), the
third binomial coefficient
$\binom{V-(V_{\bullet}+V_{\circ})}{L-(l_{\bullet}+l_{\circ})}$ enumerates the
number of ways the remaining $L-(l_{\bullet}+l_{\circ})$ links can be
redistributed _between_ the first and the second module, i.e. over the
remaining $V-(V_{\bullet}+V_{\circ})$ node pairs.

This choice induces the definition of the _binary bimodular surprise_

$\mathscr{S}_{\sslash}\equiv\sum_{l_{\bullet}\geq
l_{\bullet}^{*}}\sum_{l_{\circ}\geq l_{\circ}^{*}}f(l_{\bullet},l_{\circ});$
(11)

analogously to the univariate case, $l_{\bullet}$ and $l_{\circ}$ are
‘naturally’ bounded by the values $V_{\bullet}$ and $V_{\circ}$ \- notice,
however, that the sum $l_{\bullet}+l_{\circ}$ cannot exceed $L$ (although it
may not reach such a value, e.g. in case $V_{\bullet}+V_{\circ}<L$).

Asymptotic results. Analogously to the univariate case, the asymptotic
expression for $\mathscr{S}_{\sslash}$ can be derived upon Stirling-
approximating the binomial coefficients appearing within it. This time, the
recipe $n!\simeq\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}$ leads to

$\mathscr{S}_{\sslash}\simeq\sum_{l_{\bullet}\geq
l^{*}_{\bullet}}\sum_{l_{\circ}\geq
l^{*}_{\circ}}B(l_{\bullet},l_{\circ})\left[\frac{\text{Ber}(V,L,p)}{\prod_{i=\bullet,\circ,\top}\text{Ber}(V_{i},l_{i},p_{i})}\right]$
(12)

where, as before, $\text{Ber}(x,y,z)=z^{y}(1-z)^{x-y}$ defines a Bernoulli
probability mass function and the parameters read $p=\frac{L}{V}$ and
$p_{i}=\frac{l_{i}}{V_{i}}$; the numerical coefficient appearing in front of
the whole expression, now, reads

$B(l_{\bullet},l_{\circ})=\frac{1}{2\pi}\sqrt{\frac{\sigma^{2}}{\prod_{i=\bullet,\circ,\top}\sigma^{2}_{i}}}$
(13)

with $\sigma^{2}=Vp(1-p)$ and $\sigma^{2}_{i}=V_{i}p_{i}(1-p_{i})$. The
quantity $\mathscr{S}_{\sslash}$ compares the description of a networked
configuration provided by the RGM, and encoded into the expression
$\text{Ber}(V,L,p)=p^{L}(1-p)^{V-L}$, with the description of the same
configuration provided by the SBM (now, defined by three - instead of two -
different blocks), ‘represented’ by the denominator of the expression defined
in eq. (12), i.e.

$\prod_{i=\bullet,\circ,\top}\text{Ber}(V_{i},l_{i},p_{i})=p_{\bullet}^{l_{\bullet}}(1-p_{\bullet})^{V_{\bullet}-l_{\bullet}}\cdot
p_{\circ}^{l_{\circ}}(1-p_{\circ})^{V_{\circ}-l_{\circ}}\cdot
p_{\top}^{l_{\top}}(1-p_{\top})^{V_{\top}-l_{\top}}.$ (14)

## DETECTION OF MESOSCALE STRUCTURES
IN WEIGHTED NETWORKS

Community detection. Within the surprise-based framework, the problem of
detecting binary communities has been rephrased as an aleatory experiment
whose random variable is the number of links within communities. Interestingly
enough, such an experiment can be easily mapped into a counting problem,
allowing us to interpret $\mathscr{S}$ as indicating the number of
configurations whose number of ‘internal’ links (i.e. within communities) is
larger than the observed one.

When dealing with weighted networks, we would like to proceed along similar
guidelines and consider the total, ‘internal’ weight as our new random
variable, to be redistributed across the available node pairs. Adopting this
approach has four major consequences: 1) weights must be considered as
composed by an integer number of binary links, 2) each node pair must be
allowed to be occupied by more than one link and 3) the total weight must be
allowed to vary even beyond the network size (when handling real-world
networks, the case $W\gg V$ is often encountered).

In this case, the proper setting to define an aleatory experiment satisfying
the requests above is provided by the so-called _stars and bars_ model, a
combinatorial technique that has been introduced to handle the counting of
configurations with multiple occupancies. Basically, the problem of counting
in how many ways $w_{\bullet}$ particles (our links) can be redistributed
among $V_{\bullet}$ boxes (our node pairs), while allowing more than one
particle to occupy each box, can be tackled by allowing _both_ the particles
_and_ the bars ‘delimiting’ the boxes to be permuted Feller1950 . Since
$V_{\bullet}$ boxes are delimited by $V_{\bullet}-1$ bars, a term like
$\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}$ is needed.

In order to better grasp the meaning of such a term, let us make a simple
example. Let us imagine to observe a network with three nodes and two links,
carrying a weight of 1 and 2, respectively. Now, were we interested in a
purely binary analysis, we may ask ourselves in how many ways we could place
the two links among the $\frac{N(N-1)}{2}=\frac{3(3-1)}{2}=3$ available pairs:
the answer is provided by the ‘binary’ binomial coefficient
$\binom{V_{\bullet}}{l_{\bullet}}=\binom{3}{2}=3$. The implicit assumption we
make is that the three links must not occupy the same node pairs - otherwise
the total number of connections wouldn’t be preserved.

This perspective changes from the purely weighted point of view. Since we are
now interested in preserving just the total weight of our network,
irrespectively of the number of connections it is placed upon, the number of
admissible configurations amounts precisely at
$\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}=\binom{2+3}{3}=10$. Such a
number is larger than before since, now, weights are ‘disaggregated’ into
binary links and multiple occupations of the latter ones are allowed.

The considerations above lead us to generalize the community detection problem
to the weighted case by identifying

$\displaystyle f(w_{\bullet})$ $\displaystyle\equiv$
$\displaystyle\text{NH}(w_{\bullet}|V+W,W,V_{\bullet})$ (15) $\displaystyle=$
$\displaystyle\frac{\prod_{i=\bullet,\circ}\binom{V_{i}+w_{i}-1}{w_{i}}}{\binom{V+W-1}{W}}=\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}}{\binom{V+W-1}{W}}$
$\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{(V-V_{\bullet})+(W-w_{\bullet})-1}{W-w_{\bullet}}}{\binom{V+W-1}{W}}$

i.e. by replacing the binomial hypergeometric distribution considered in the
purely binary case with a _negative hypergeometric distribution_ , a choice
inducing the definition of the _weighted surprise_

$\mathscr{W}\equiv\sum_{w_{\bullet}\geq w_{\bullet}^{*}}f(w_{\bullet})$ (16)

where the binomial coefficient
$\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}=\binom{V_{\bullet}+w_{\bullet}-1}{V_{\bullet}-1}$
enumerates the number of ways $w_{\bullet}$ links can be redistributed
_within_ communities, i.e. over the available $V_{\bullet}$ node pairs, and
the binomial coefficient
$\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}=\binom{V_{\circ}+w_{\circ}-1}{V_{\circ}-1}$
enumerates the number of ways the remaining $w_{\circ}=W-w_{\bullet}$ links
can be redistributed _between_ communities, i.e. over the remaining
$V_{\circ}=V-V_{\bullet}$ node pairs. Differently from the binary case, the
sum ranges up to the maximum empirical weight of the network, i.e.
$w_{\bullet}\in[w_{\bullet}^{*},W]$.

Asymptotic results. The asymptotic expression for $\mathscr{W}$ can be deduced
by following the same reasoning that has allowed us to derive the asymptotic
expression for $\mathscr{S}$. Stirling-approximating the binomial coefficients
entering into the definition of $\mathscr{W}$ leads to the writing

$\mathscr{W}\simeq\sum_{w_{\bullet}\geq
w^{*}_{\bullet}}C(w_{\bullet})\left[\frac{\text{Geo}(V,W,q)}{\prod_{i=\bullet,\circ}\text{Geo}(V_{i},w_{i},q_{i})}\right]$
(17)

where the expression

$\text{Geo}(x,y,z)=z^{y}(1-z)^{x}$ (18)

defines a _geometric_ probability mass function and the parameters appearing
in eq. (17) read $q=\frac{W}{V+W-1}$ and $q_{i}=\frac{w_{i}}{V_{i}+w_{i}-1}$.
In the weighted case, the Bernoulli probability mass function appearing in the
asymptotic expression of $\mathscr{S}$ is replaced by a geometric probability
mass function: this implies that (asymptotically) the comparison is, now,
carried out between the description of a networked configuration provided by
the Weighted Random Graph Model (WRGM), and encoded into the expression

$\text{Geo}(V,W,q)=q^{W}(1-q)^{V},$ (19)

with the description of the same configuration provided by the Weighted
Stochastic Block Model (WSBM) and encoded into the expression

$\prod_{i=\bullet,\circ}\text{Geo}(V_{i},w_{i},q_{i})=q_{\bullet}^{w_{\bullet}}(1-q_{\bullet})^{V_{\bullet}}\cdot
q_{\circ}^{w_{\circ}}(1-q_{\circ})^{V_{\circ}}$ (20)

where $q_{\bullet}=\frac{w_{\bullet}}{V_{\bullet}+w_{\bullet}}$ and
$q_{\circ}=\frac{w_{\circ}}{V_{\circ}+w_{\circ}}$. As for the binary case, the
asymptotic expression of the weighted surprise clarifies that minimizing it
amounts at finding the partition least likely to occour under the WRGM than
under the WSBM.

As in the binary case, the parameters characterizing the geometric probability
mass functions defining the asymptotic weighted surprise can be estimated via
the maximum-of-the-likelihood principle, according to which the log-likelihood
of the expression $\text{Geo}(x,y,z)$ reads

$\mathscr{L}=y\ln z+x\ln(1-z);$ (21)

upon maximizing it with respect to $z$, one finds $z=\frac{y}{x+y}$: hence,
one can pose $q\simeq\frac{W}{V+W}$ and
$q_{i}\simeq\frac{w_{i}}{V_{i}+w_{i}}$.

The numerical coefficient appearing in front of the whole expression, instead,
reads

$C(w_{\bullet})=\sqrt{\frac{\mu}{2\pi\prod_{i=\bullet,\circ}\mu_{i}}}$ (22)

with $\mu\simeq Vq$ and $\mu_{i}\simeq V_{i}q_{i}$.

‘Bimodular’ structures detection. Let us now introduce the third
generalization of the suprise-based formalism: following the same line of
reasoning that led us to approach the detection of binary ‘bimodular’
structures by considering the multinomial analogoue of the distribution
introduced for binary community detection, we focus on the _multinomial_ (or
_multivariate_) _negative hypergeometric distribution_ , i.e.

$\displaystyle f(w_{\bullet},w_{\circ})$ $\displaystyle\equiv$
$\displaystyle\text{MNH}(w_{\bullet},w_{\circ}|V+W,W,V_{\bullet},V_{\circ})$
$\displaystyle=$
$\displaystyle\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}+w_{i}-1}{w_{i}}}{\binom{V+W-1}{W}}$
$\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}\binom{V_{\top}+w_{\top}-1}{w_{\top}}}{\binom{V+W-1}{W}}$
$\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}\binom{V-(V_{\bullet}+V_{\circ})+W-(w_{\bullet}+w_{\circ})-1}{W-(w_{\bullet}+w_{\circ})}}{\binom{V+W-1}{W}};$

while the binomial coefficient
$\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}$ enumerates the number of ways
$w_{\bullet}$ links can redistributed _within_ the first module (e.g. the core
portion), the binomial coefficient $\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}$
enumerates the number of ways $w_{\circ}$ links can be redistributed _within_
the second module (e.g. the periphery portion) and the binomial coefficient
$\binom{V-(V_{\bullet}+V_{\circ})+W-(w_{\bullet}+w_{\circ})-1}{W-(w_{\bullet}+w_{\circ})}$
enumerates the number of ways the remaining $w_{\top}\equiv
W-(w_{\bullet}+w_{\circ})$ links can be redistributed _between_ the first and
the second module, i.e. over the remaining $V_{\top}\equiv
V-(V_{\bullet}+V_{\circ})$ node pairs. Such a position induces the definition
of the _weighted bimodular surprise_

$\mathscr{W}_{\sslash}\equiv\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}\sum_{w_{\circ}\geq w_{\circ}^{*}}f(w_{\bullet},w_{\circ});$
(24)

as for the (weighted) community detection, weights are understood as integer
numbers - equivalently, as composed by an integer number of binary links. For
what concerns the limits of the summations, $w_{\bullet}$ and $w_{\circ}$ are
‘naturally’ bounded by $W$; notice, however, that the sum
$w_{\bullet}+w_{\circ}$ itself cannot exceed such a value.

Asymptotic results. Let us now derive the asymptotic expression for
$\mathscr{W}_{\sslash}$. As usual, let us Stirling-approximate the binomial
coefficients entering into the definition of $\mathscr{W}_{\sslash}$; such a
simplification leads us to the expression

$\mathscr{W}_{\sslash}\simeq\sum_{w_{\bullet}\geq
w^{*}_{\bullet}}\sum_{w_{\circ}\geq
w^{*}_{\circ}}D(w_{\bullet},w_{\circ})\left[\frac{\text{Geo}(V,W,q)}{\prod_{i=\bullet,\circ,\top}\text{Geo}(V_{i},w_{i},q_{i})}\right]$
(25)

where, as before, $\text{Geo}(x,y,z)=z^{y}(1-z)^{x}$ defines a geometric
probability mass function and the parameters read $q=\frac{W}{V+W-1}$ and
$q_{i}=\frac{w_{i}}{V_{i}+w_{i}-1}$ but can be approximated, according to the
maximum-of-the-likelihood principle, as $q\simeq\frac{W}{V+W}$ and
$q_{i}\simeq\frac{w_{i}}{V_{i}+w_{i}}$; the numerical coefficient multiplying
the whole expression reads

$D(w_{\bullet},w_{\circ})=\frac{1}{2\pi}\sqrt{\frac{\mu}{\prod_{i=\bullet,\circ,\top}\mu_{i}}}$
(26)

with $\mu\simeq Vq$ and $\mu_{i}\simeq V_{i}q_{i}$. The analysis of
$\mathscr{W}_{\sslash}$ in the asymptotic regime reveals that it compares the
description of a networked configuration provided by the WRGM, and encoded
into the expression $\text{Geo}(V,W,q)=q^{W}(1-q)^{V}$, with the description
of the same configuration provided by the WSBM (now, defined by three -
instead of two - different blocks), ‘represented’ by the denominator of the
expression defined in eq. (25), i.e.

$\prod_{i=\bullet,\circ,\top}\text{Geo}(V_{i},w_{i},q_{i})=q_{\bullet}^{w_{\bullet}}(1-q_{\bullet})^{V_{\bullet}}\cdot
q_{\circ}^{w_{\circ}}(1-q_{\circ})^{V_{\circ}}\cdot
q_{\top}^{w_{\top}}(1-q_{\top})^{V_{\top}}.$ (27)

Figure 1: Graphical comparison of the three different ways of counting the
admissible configurations when dealing with the purely binary, the purely
weighted and the enhanced surprise. Let us imagine to observe a network with
three nodes and two links, carrying a weight of 1 and 2, respectively. Were we
interested in a purely binary analysis, we may ask ourselves in how many ways
we could place the two links among the $3$ available pairs: the answer is
provided by the ‘binary’ binomial coefficient
$\binom{V_{\bullet}}{l_{\bullet}}=\binom{3}{2}=3$. Were we interested in a
purely weighted analysis, we may ask ourselves in how many ways we could place
the two links among the $3$ available pairs while preserving the total weight
of our network, irrespectively of the number of connections it is placed upon;
the number of admissible configurations becomes
$\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}=\binom{2+3}{3}=10$. Such a
number is larger than before since, now, weights are ‘disaggregated’ into
binary links and multiple occupations of the latter ones are allowed. Were we
interested in the enhanced analysis, we may ask ourselves in how many ways we
could place the two links among the $3$ available pairs while preserving both
the total number of links and the total weight of the network; the number of
admissible configurations becomes
$\binom{V_{\bullet}}{l_{\bullet}}\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}=\binom{3}{2}\binom{3-1}{3-2}=3\cdot
2=6$, as can be easily verified upon explicitly listing them.

## ENHANCED DETECTION OF MESOSCALE
STRUCTURES IN NETWORKS

Community detection. The recipe to detect communities on weighted networks can
be further refined to account for the information encoded into the total
number of links, beside the one provided by the total weight. Generally
speaking, this can be realized by ‘combining’ two of the distributions
introduced above. To this aim, let us proceed in a two-step fashion: first,
let us recall that the number of ways $L$ links can be placed among $V$ node
pairs, in such a way that $l_{\bullet}$ connections are ‘internal’ to the
clusters while the remaining $L-l_{\bullet}$ ones are, instead, ‘external’ is
precisely

$\text{H}(l_{\bullet}|V,V_{\bullet},L)=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}};$
(28)

now, for each of the binary configurations listed above, $W-L$ links remain to
be assigned: while $w_{\bullet}-l_{\bullet}$ of them must be placed within the
clusters, on top of the $l_{\bullet}$ available ‘internal’ links, the
remaining $(W-L)-(w_{\bullet}-l_{\bullet})$ ones must be placed between the
clusters, on top of the $L-l_{\bullet}$ available inter-cluster connections.
Hence, the ‘conditional’ negative hypergeometric distribution reading

$\text{NH}(w_{\bullet}|W,W-L,l_{\bullet})=\frac{\binom{l_{\bullet}+(w_{\bullet}-l_{\bullet})-1}{w_{\bullet}-l_{\bullet}}\binom{(L-l_{\bullet})+(W-L)-(w_{\bullet}-l_{\bullet})-1}{(W-L)-(w_{\bullet}-l_{\bullet})}}{\binom{L+(W-L)-1}{W-L}}$
(29)

remains naturally defined; now, ‘combining’ the two distributions above,
simplifying and re-arranging, the generic term of the _enhanced hypergeometric
distribution_ can be rewritten as

$\displaystyle\text{EH}(l_{\bullet},w_{\bullet}|V,V_{\bullet},L,W)$
$\displaystyle=$
$\displaystyle\text{H}(l_{\bullet}|V,V_{\bullet},L)\cdot\text{NH}(w_{\bullet}|W,W-L,l_{\bullet})$
(30) $\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}\binom{w_{\circ}-1}{w_{\circ}-l_{\circ}}}{\binom{W-1}{W-L}}$

with a clear meaning of the symbols. An analytical characterization of it is
provided into the Appendices: for the moment, let us simply notice that the
definition provided above works for the values $0<l_{\bullet}<L$. By posing
$f(l_{\bullet},w_{\bullet})\equiv\text{EH}(l_{\bullet},w_{\bullet}|V,V_{\bullet},L,W)$,
our novel distribution induces the definition of the _enhanced surprise_ ,
i.e.

$\mathscr{E}\equiv\sum_{l_{\bullet}\geq l_{\bullet}^{*}}\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}f(l_{\bullet},w_{\bullet});$ (31)

although $l_{\bullet}$ and $w_{\bullet}-l_{\bullet}$ are ‘naturally’ bounded
by $V_{\bullet}$ and $W-L$, respectively, the former one cannot exceed $L$.

In order to better understand how the enhanced surprise works, let us consider
again the aforementioned example: given a network with three nodes and two
links, carrying a weight of 1 and 2, respectively, we observe
$\binom{V_{\bullet}}{l_{\bullet}}=\binom{3}{2}=3$ (purely binary)
configurations with exactly the same number of links and
$\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}=\binom{2+3}{3}=10$ (purely
weighted) configurations with exactly the same total weight. If we, now,
constrain both the total number of links and the total weight of the network,
the number of admissible configurations becomes
$\binom{V_{\bullet}}{l_{\bullet}}\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}=\binom{3}{2}\binom{3-1}{3-2}=3\cdot
2=6$, as it can be easily verified upon explicitly listing them. Naturally,
the configurations ‘admissible’ by the enhanced surprise are a subset of the
configurations ‘admissible’ by the weighted surprise, i.e. precisely the ones
with the desired number of links (see also fig. 1).

Asymptotic results. Analogously to the other functionals, the asymptotic
expression of $\mathscr{E}$ can be derived by Stirling-approximating the
binomial coefficients entering into its definition as well:

$\mathscr{E}\simeq\sum_{l_{\bullet}\geq l_{\bullet}^{*}}\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}E(l_{\bullet},w_{\bullet})\left[\frac{\text{BF}(V,L,W,p,r)}{\prod_{i=\bullet,\circ}\text{BF}(V_{i},l_{i},w_{i},p_{i},r_{i})}\right];$
(32)

in the formula above, the expression

$\displaystyle\text{BF}(x,y,u,z,t)$ $\displaystyle=$ $\displaystyle
z^{y}(1-z)^{x-y}\cdot t^{u-y}(1-t)^{y}$ (33) $\displaystyle=$
$\displaystyle\text{Ber}(x,y,z)\cdot\text{Geo}(y,u-y,t)$

defines a _Bose-Fermi_ probability mass function Garlaschelli2009 . While a
Bernoulli probability mass function characterizes the asymptotic behavior of
the purely binary surprise and a geometric probability mass function
characterizes the asymptotic behavior of the purely weighted surprise, the
enhanced surprise is asymptotically characterized by a distribution whose
functional form is halfway between the two previous ones: while its Bernoulli-
like portion controls for the ‘presence’ of the links, its ‘conditional’
geometric-like portion controls for the magnitude of the ‘remaining’ weights.
To stress its ‘mixed’ character, such a distribution has been named ‘Bose-
Fermi’: remarkably, it can be retrieved within the Exponential Random Graphs
framework as a consequence of Shannon entropy maximization, constrained to
simultaneously reproduce both the total number of links and the total weight
of a network Garlaschelli2009 ; Mastrandrea2014 .

The parameters appearing in eq. (32) read $p=\frac{L}{V}$,
$p_{i}=\frac{l_{i}}{V_{i}}$ and $r=\frac{W-L}{W-1}$,
$r_{i}=\frac{w_{i}-l_{i}}{w_{i}-1}$; while the first class of parameters can
be tuned according to the maximum-of-the-likelihood principle, the ones
belonging to the second class can be approximated according to the same
recipe. To see this explicitly, let us consider $\text{BF}(x,y,u,z,t)$, whose
log-likelihood reads

$\mathscr{L}=y\ln z+(x-y)\ln(1-z)+(u-y)\ln t+y\ln(1-t);$ (34)

upon maximizing it with respect to $z$, one finds $z=\frac{y}{x}$; upon
maximizing it with respect to $t$ one finds $t=\frac{u-y}{u}$ \- hence, one
can pose $r\simeq\frac{W-L}{W}$, $r_{i}\simeq\frac{w_{i}-l_{i}}{w_{i}}$.

The numerical coefficient multiplying the whole expression is defined as

$E(l_{\bullet},w_{\bullet})=\frac{1}{2\pi}\sqrt{\frac{\sigma^{2}\mu}{\prod_{i=\bullet,\circ}\sigma^{2}_{i}\mu_{i}}}$
(35)

where $\sigma^{2}=Vp(1-p)$, $\sigma^{2}_{i}=V_{i}p_{i}(1-p_{i})$, $\mu\simeq
Lr$ and $\mu_{i}\simeq l_{i}r_{i}$.

Similarly to the other cases, the asymptotic expression of the enhanced
suprise compares the description of a networked configuration provided by the
Enhanced Random Graph Model (ERGM), and encoded into the expression
$\text{BF}(V,L,W,p,r)=p^{L}(1-p)^{V-L}\cdot r^{W-L}(1-r)^{L}$, with the
description of the same configuration provided by its block-wise counterpart,
i.e. the Enhanced Stochastic Block Model (ESBM) and encoded into the
expression $\prod_{i=\bullet,\circ}\text{BF}(V_{i},l_{i},w_{i},p_{i},r_{i})$.

‘Bimodular’ structures detection. The last generalization of surprise concerns
its use for the detection of ‘bimodular’ structures within the enhanced
framework. This amounts at considering the following ‘multinomial variant’ of
the enhanced hypergeometric distribution

$\displaystyle\text{MEH}(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ}|V,V_{\bullet},V_{\circ},L,W)$
$\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V_{\top}}{l_{\top}}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}\binom{w_{\circ}-1}{w_{\circ}-l_{\circ}}\binom{w_{\top}-1}{w_{\top}-l_{\top}}}{\binom{W-1}{W-L}}$
(36) $\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V-(V_{\bullet}+V_{\circ})}{L-(l_{\bullet}+l_{\circ})}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}\binom{w_{\circ}-1}{w_{\circ}-l_{\circ}}\binom{W-(w_{\bullet}+w_{\circ})-1}{(W-L)-((w_{\bullet}+w_{\circ})-(l_{\bullet}+l_{\circ}))}}{\binom{W-1}{L-1}}$

where $V_{\top}\equiv V-(V_{\bullet}+V_{\circ})$ indicates the number of node
pairs between the modules $\bullet$ and $\circ$ and $l_{\top}\equiv
L-(l_{\bullet}+l_{\circ})$ indicates the number of links that must be assigned
therein. An analytical characterization of it is provided into the Appendices:
for the moment, let us simply notice that the definition provided above works
for the values $0<l_{\bullet},l_{\circ}<L$. The position
$f(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ})\equiv\text{MEH}(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ}|V,V_{\bullet},V_{\circ},L,W)$
induces the definition of the _enhanced bimodular surprise_

$\mathscr{E}_{\sslash}=\sum_{l_{\bullet}\geq
l_{\bullet}^{*}}\sum_{l_{\circ}\geq l_{\circ}^{*}}\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}\sum_{w_{\circ}\geq
w_{\circ}^{*}}f(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ}).$ (37)

Notice that $l_{\bullet}$ and $l_{\circ}$ are ‘naturally’ bounded by
$V_{\bullet}$ and $V_{\circ}$: still, their sum cannot exceed $L$;
analogously, $w_{\bullet}-l_{\bullet}$ and $w_{\circ}-l_{\circ}$ are
‘naturally’ bounded by $W-L$: still, their sum itself cannot exceed $W-L$.

As for its binomial counterpart, the expression of the MEH can be rearranged
in a term-by-term fashion, in such a way that the module-specific binomial
coefficients can be grouped together. Upon doing so, it becomes clearer that
the MEH counts the number of ways $w_{\bullet}-l_{\bullet}$ links can be
placed on top of the $l_{\bullet}$ binary links characterizing the connectance
of the $\bullet$ module, times the number of ways $w_{\circ}-l_{\circ}$ links
can be placed on top of the $l_{\circ}$ binary links characterizing the
connectance of the $\circ$ module, times the number of ways the remaining
$W-(w_{\bullet}+w_{\circ})-(L-(l_{\bullet}+l_{\circ}))$ links can be placed on
top of the $L-(l_{\bullet}+l_{\circ})$ binary links characterizing the
connectance of the third module.

Full expression | Community detection | ‘Bimodular’ structures detection
---|---|---
Binary case | $f=\frac{\prod_{i=\bullet,\circ}\binom{V_{i}}{l_{i}}}{\binom{V}{L}}$ | $f=\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}}{l_{i}}}{\binom{V}{L}}$
Weighted case | $f=\frac{\prod_{i=\bullet,\circ}\binom{V_{i}+w_{i}-1}{w_{i}}}{\binom{V+W-1}{W}}=\frac{\prod_{i=\bullet,\circ}\left(\kern-2.10002pt\left(\genfrac{}{}{0.0pt}{}{V_{i}}{w_{i}}\right)\kern-2.10002pt\right)}{\left(\kern-2.10002pt\left(\genfrac{}{}{0.0pt}{}{V}{W}\right)\kern-2.10002pt\right)}$ | $f=\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}+w_{i}-1}{w_{i}}}{\binom{V+W-1}{W}}=\frac{\prod_{i=\bullet,\circ,\top}\left(\kern-2.10002pt\left(\genfrac{}{}{0.0pt}{}{V_{i}}{w_{i}}\right)\kern-2.10002pt\right)}{\left(\kern-2.10002pt\left(\genfrac{}{}{0.0pt}{}{V}{W}\right)\kern-2.10002pt\right)}$
Enhanced case | $f=\frac{\prod_{i=\bullet,\circ}\binom{V_{i}}{l_{i}}\binom{w_{i}-1}{l_{i}-1}}{\binom{V}{L}\binom{W-1}{W-L}}$ | $f=\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}}{l_{i}}\binom{w_{i}-1}{l_{i}-1}}{\binom{V}{L}\binom{W-1}{W-L}}$
Asymptotic expression | Community detection | ‘Bimodular’ structures detection
Binary case | $f\propto\frac{\text{Ber}(V,L,p)}{\prod_{i=\bullet,\circ}\text{Ber}(V_{i},l_{i},p_{i})}$ | $f\propto\frac{\text{Ber}(V,L,p)}{\prod_{i=\bullet,\circ,\top}\text{Ber}(V_{i},l_{i},p_{i})}$
Weighted case | $f\propto\frac{\text{Geo}(V,W,q)}{\prod_{i=\bullet,\circ}\text{Geo}(V_{i},w_{i},q_{i})}$ | $f\propto\frac{\text{Geo}(V,W,q)}{\prod_{i=\bullet,\circ,\top}\text{Geo}(V_{i},w_{i},q_{i})}$
Enhanced case | $f\propto\frac{\text{BF}(V,L,W,p,r)}{\prod_{i=\bullet,\circ}\text{BF}(V_{i},l_{i},w_{i},p_{i},r_{i})}$ | $f\propto\frac{\text{BF}(V,L,W,p,r)}{\prod_{i=\bullet,\circ,\top}\text{BF}(V_{i},l_{i},w_{i},p_{i},r_{i})}$
Table 1: Table illustrating all the generalizations of the surprise-based
formalism proposed in the present paper and showing both the full and the
asymptotic expression of each probability mass function - with the only
exception of the numerical coefficient characterizing each asymptotic
expression. To sum up, detecting a weighted mesoscale structure implies
considering the _negative_ version of the probability mass function working in
the corresponding binary case (e.g. moving from the hypergeometric to the
negative hypergeometric one); detecting a ‘bimodular’ structure, instead,
implies considering the _multinomial_ version of the probability mass function
working in the corresponding binary case (e.g. moving from the binomial
hypergeometric to the multinomial one). Upon employing the multiset notation,
a nice formal symmetry can be recovered between the purely binary and the
purely weighted cases.

Asymptotic results. The enhanced bimodular surprise admits an asymptotic
expression as well, i.e.

$\mathscr{E}_{\sslash}\simeq\sum_{l_{\bullet}\geq
l_{\bullet}^{*}}\sum_{l_{\circ}\geq l_{\circ}^{*}}\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}\sum_{w_{\circ}\geq
w_{\circ}^{*}}F(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ})\left[\frac{\text{BF}(V,L,W,p,r)}{\prod_{i=\bullet,\circ,\top}\text{BF}(V_{i},l_{i},w_{i},p_{i},r_{i})}\right]$
(38)

that can be recovered by Stirling-approximating the binomial coefficients
entering into the definition of $\mathscr{E}_{\sslash}$. While the expression
BF denotes a Bose-Fermi probability mass function Garlaschelli2009 , the
parameters appearing in eq. (38) read $p=\frac{L}{V}$,
$p_{i}=\frac{l_{i}}{V_{i}}$ and $r=\frac{W-L}{W-1}$,
$r_{i}=\frac{w_{i}-l_{i}}{w_{i}-1}$; as for $\mathscr{E}$, the maximum-of-the-
likelihood principle determines (approximates) the first (second) class of
parameters.

The numerical coefficient multiplying the whole expression is defined as

$F(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ})=\frac{1}{(2\pi)^{2}}\sqrt{\frac{\sigma^{2}\mu}{\prod_{i=\bullet,\circ,\top}\sigma^{2}_{i}\mu_{i}}}$
(39)

where $\sigma^{2}=Vp(1-p)$, $\sigma^{2}_{i}=V_{i}p_{i}(1-p_{i})$, $\mu\simeq
Lr$ and $\mu_{i}\simeq l_{i}r_{i}$. As observed for the other functionals, the
asymptotic expression of the enhanced bimodular suprise compares the
description of a networked configuration provided by the ERGM, and encoded
into the expression $\text{BF}(V,L,W,p,r)=p^{L}(1-p)^{V-L}\cdot
r^{W-L}(1-r)^{L}$, with the description of the same configuration provided by
its block-wise counterpart, i.e. the ESBM (now, defined by three - instead of
two - different blocks), ‘represented’ by the denominator of the expression
defined in eq. (38), i.e.
$\prod_{i=\bullet,\circ,\top}\text{BF}(V_{i},l_{i},w_{i},p_{i},r_{i})$

Table 2 gathers all the variants of the surprise-based formalism, illustrating
both the full and the asymptotic expression for each of them. To sum up,
detecting a weighted mesoscale structure implies considering the _negative_
version of the probability mass function working in the corresponding binary
case (e.g. moving from the hypergeometric to the negative hypergeometric one);
detecting a ‘bimodular’ structure, instead, implies considering the
_multinomial_ version of the probability mass function working in the
corresponding binary case (e.g. moving from the binomial hypergeometric to the
multinomial one).

Figure 2: Consistency check of our surprise-based formalism in the simplest
case of community detection on binary networks. Any exact statistical test is
characterized by a parameter, usually denoted with $\alpha$ and known as the
_type I error rate_ , that quantifies the percentage of times the test
provides a _false positive_ : it can be kept below a given threshold by
adjusting the threshold of the p-value of the corresponding test, e.g. by
deeming its response as significant in case it is less than $0.05$. Over
10.000 networks sampled by the RGM with $p=0.2$, a given planted partition -
from left to right: two, three and four clusters with different dimensions -
is indeed recovered as significant $5\%$ of the times. Notably, since the
three upper panels show the CDF of the empirical values of $\mathscr{S}$, they
also provide an information about its distribution, which is uniform over the
unit interval.

## RESULTS

The previous sections have been devoted to the description of the surprise-
based formalism for detecting a number of mesoscale structures; let us now
test it on a bunch of synthetic and real-world configurations.

Consistency checks. Let us start by checking the consistency of our surprise-
based formalism. Since we have rephrased the problem of detecting any
mesoscale structure into an exact significance test, the limitations of our
formalism are the same ones affecting the tests of the kind. More precisely,
any significance test is characterized by a parameter known as _type I error
rate_ and usually denoted with $\alpha$: it quantifies the percentage of times
the considered test provides a _false positive_. In other words, any
statistical test is known to ‘fail’; in our case, this amounts at recognizing
that a significant mesoscale structure can be detected, _in case there is
none_ , a percentage $\alpha$ of the times: this number can be kept ‘small’ by
adjusting the threshold of the p-value of the corresponding test - e.g. by
deeming its response as significant in case it is less than $0.05$.

Remarkably, the aforementioned behavior can be explicitly tested for each of
the cases considered above. Whenever exercises like these are carried out, it
is of utmost importance to be as clear as possible about the null hypothesis
tested: here, we aim at testing whether _a given partition_ is significant or
not when the null hypothesis is true. For the sake of illustration, let us
consider the problem of detecting communities on binary networks: as we learnt
from the asymptotic expression of $\mathscr{S}$, it compares the description
of a network provided by the RGM with that of the same network provided by the
SBM, thus suggesting the RGM as the model playing the role of $H_{0}$. Hence,
let us generate many networks from the RGM and, for each of them, let us
impose a partition - the same for each sampled configuration, over which our
SBM is tuned; finally, let us calculate $\mathscr{S}$ on each of them.

The results of such an experiment are shown in fig. 2: over 10.000 networks
sampled from the RGM, whose only parameter has been set to $p=0.2$, the
(given) planted partition - i.e. two, three and four clusters with different
dimensions - is, indeed, recovered as significant $5\%$ of the times; notice
that such a result holds true irrespectively of the details of the partition
imposed on the sampled configurations.

We have also repeated such an experiment for the problem of detecting
communities on weighted networks: as fig. 3 shows, the same results are
recovered, i.e. over 10.000 networks sampled from the WRGM, whose only
parameter has been set to $q=0.2$, a (given) planted partition of two, three
and four clusters with different dimensions is recognized as significant $5\%$
of the times; again, such a result holds true irrespectively of the details of
the partition imposed on our sampled configurations.

Figure 3: Consistency check of our surprise-based formalism in the case of
community detection on weighted networks. As in the binary case, the type I
error rate can be kept below a given threshold by adjusting the threshold of
the p-value of the corresponding test, e.g. by deeming its response as
significant in case it is less than $0.05$. Over 10.000 networks sampled by
the WRGM with $q=0.2$, a given planted partition - from left to right: two,
three and four clusters with different dimensions - is indeed recovered as
significant $5\%$ of the times. Notably, since the three upper panels show the
CDF of the empirical values of $\mathscr{W}$, they also provide an information
about its distribution, which is uniform over the unit interval.

Our experiments tell us something deeper about the behavior of $\mathscr{S}$
and $\mathscr{W}$: since each of the three upper panels in figs. 2 and 3 shows
the cumulative distribution functions (CDF) of the empirical values of
$\mathscr{S}$ and $\mathscr{W}$, we learn that the latter ones are distributed
_uniformly_ over the unit interval, i.e. $\mathscr{S}\sim\text{U}[0,1]$ and
$\mathscr{W}\sim\text{U}[0,1]$ \- an evidence further confirming the behavior
of the p-value of any exact significance test.

Figure 4: Upper row: performance of modularity, surprise and Infomap in
correctly identifying the partitions induced by 7, different, ‘homogeneous’
rings of cliques. For each configuration, 20 cliques have been considered; the
size of the latter ones is left to vary as follows: $K_{3}$, $K_{4}$, $K_{5}$,
$K_{8}$, $K_{10}$, $K_{15}$, $K_{20}$. While surprise always recovers the
planted partition, modularity misses the partitions with $K_{3}$ and $K_{4}$
and Infomap misses the partition with $K_{3}$ \- a result that may be a
consequence of the resolution limit which is known to affect the last two
algorithms. Bottom row: comparison between the ‘ring of binary cliques’ and
the ‘ring of weighted cliques’ cases. The result according to which surprise
minimization is able to discriminate the cliques linked in a ring-like fashion
changes once the weights come into play: in fact, as the weight of the links
connecting any two cliques is risen, the algorithm reveals as ‘communities’
pairs of tightly-connected cliques.

Checking the behavior of the multivariate versions of our hypergeometric
distributions is, instead, much more difficult, the reason lying in the
evidence that the null hypothesis is, now, a multivariate one and few results
about the behavior of multivariate p-values are known. Still, we can say
something about the behavior of the _marginal p-values_ ; for the sake of
illustration, let us consider the ones induced by $\mathscr{S}_{\sslash}$ and
$\mathscr{W}_{\sslash}$ and reading

$\displaystyle\mathscr{S}_{\sslash}^{(\bullet)}$ $\displaystyle=$
$\displaystyle\sum_{l_{\bullet}\geq
l_{\bullet}^{*}}\sum_{l_{\circ}}f(l_{\bullet},l_{\circ})=\sum_{l_{\bullet}\geq
l_{\bullet}^{*}}f(l_{\bullet}),$ (40)
$\displaystyle\mathscr{S}_{\sslash}^{(\circ)}$ $\displaystyle=$
$\displaystyle\sum_{l_{\circ}\geq
l_{\circ}^{*}}\sum_{l_{\bullet}}f(l_{\bullet},l_{\circ})=\sum_{l_{\circ}\geq
l_{\circ}^{*}}f(l_{\circ})$ (41)

and

$\displaystyle\mathscr{W}_{\sslash}^{(\bullet)}$ $\displaystyle=$
$\displaystyle\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}\sum_{w_{\circ}}f(w_{\bullet},w_{\circ})=\sum_{w_{\bullet}\geq
w_{\bullet}^{*}}f(w_{\bullet}),$ (42)
$\displaystyle\mathscr{W}_{\sslash}^{(\circ)}$ $\displaystyle=$
$\displaystyle\sum_{w_{\circ}\geq
w_{\circ}^{*}}\sum_{w_{\bullet}}f(w_{\bullet},w_{\circ})=\sum_{w_{\circ}\geq
w_{\circ}^{*}}f(w_{\circ})$ (43)

respectively. As evident from the formulas above, marginal distributions of
multivariate hypergeometric distributions are hypergeometric themselves:
hence, the behavior of our marginal p-values is expected to match the one
observed in the univariate cases, leading to recover
$\mathscr{S}_{\sslash}^{(\bullet)}\sim\text{U}[0,1]$,
$\mathscr{S}_{\sslash}^{(\circ)}\sim\text{U}[0,1]$,
$\mathscr{W}_{\sslash}^{(\bullet)}\sim\text{U}[0,1]$ and
$\mathscr{W}_{\sslash}^{(\circ)}\sim\text{U}[0,1]$. Analogously, for the
enhanced surprise.

Figure 5: Examples of the benchmarks used in the present paper to compare
modularity, surprise and Infomap: Lancichinetti-Fortunato-Radicchi benchmark
(LFR, left panel), constituted by 1000 nodes arranged in ‘big’ communities and
with mixing coefficient $\mu_{t}=0.1$; Aldecoa’s ‘relaxed-caveman’ benchmark
(RC, right panel), constituted by 16 communities whose size is distributed
according to a power-law and with ‘degradation’ coefficient $p=0.1$.

Our findings can be also described from a different perspective, i.e. by
answering the question _what is the expected value of the surprise on a
network sampled from the RGM ensemble?_ To answer this question let us focus
on the simplest case of detecting communities on binary networks and consider
that, irrespectively from the partition that is planted on the sampled
configuration, the expected number of links will be $\langle L\rangle=pV$
while the expected number of ‘internal’ links will be $\langle
l_{\bullet}^{*}\rangle=pV_{\bullet}$, where $p$ is the parameter defining the
RGM; hence, the surprise becomes

$\mathscr{S}=\sum_{l_{\bullet}\geq\langle
l_{\bullet}^{*}\rangle}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{\langle
L\rangle-l_{\bullet}}}{\binom{V}{\langle L\rangle}}=\sum_{l_{\bullet}\geq
pV_{\bullet}}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{pV-
l_{\bullet}}}{\binom{V}{pV}}$ (44)

(to be fully consistent we should have rounded both the value of $\langle
L\rangle$ and that of $\langle l_{\bullet}^{*}\rangle$ to the nearest integer
but our conclusions are still valid). In order to evaluate the expression
above, let us consider that the expected value of the hypergeometric
distribution defining the surprise reads

$\langle l_{\bullet}\rangle=\langle
L\rangle\frac{V_{\bullet}}{V}=pV_{\bullet}$ (45)

the validity of the first passage resting upon the evidence that the network
is generated via the RGM. Hence, the surprise becomes a sum over the values
$l_{\bullet}\geq\langle l_{\bullet}\rangle$, i.e. those that are more extreme
that the average. Since the hypergeometric distribution is peaked around its
average, the result above suggests that summing over larger values encodes
half of the probability mass, i.e. $\simeq\frac{1}{2}$: as a consequence, the
expected value of surprise, on a network generated via the RGM, amounts at
$\simeq\frac{1}{2}$, irrespectively from the partition that is planted on the
sampled configuration; its fluctuations are, instead, compatible with those of
a uniform distribution whose support coincides with the unit interval.

Comparison with the binary modularity. For the sake of comparison, let us
consider the modularity function for the problem of community detection on
binary, undirected networks. It is defined as

${\color[rgb]{0,0,0}Q=\frac{1}{2L}\sum_{i\neq j}(a_{ij}-\langle
a_{ij}\rangle)\delta_{c_{i},c_{j}}}$ (46)

where, for consistency, we employ the RGM as a benchmark, i.e. $\langle
a_{ij}\rangle=p=\frac{L}{V}=\frac{2L}{N(N-1)}$, $\forall\>i<j$. From the
definition of modularity, it follows that

$\displaystyle Q$ $\displaystyle=\frac{1}{2L}\sum_{i\neq
j}a_{ij}\delta_{c_{i},c_{j}}-\frac{1}{2L}\sum_{i\neq j}p\delta_{c_{i},c_{j}}$
(47) $\displaystyle=\frac{1}{2L}\sum_{c}[2l_{c}^{*}-pN_{c}(N_{c}-1)]$ (48)
$\displaystyle=\frac{1}{2L}\sum_{c}\left[2l_{c}^{*}-2L\frac{N_{c}(N_{c}-1)}{N(N-1)}\right]$
(49)
$\displaystyle=\left(\frac{l_{\bullet}^{*}}{L}-\frac{V_{\bullet}}{V}\right)$
(50)

where $N_{c}$ and $l_{c}^{*}$ indicate the size of community $c$ and the
number of links _within_ it, respectively - naturally,
$l_{\bullet}^{*}=\sum_{c}l_{c}^{*}$ and $V_{\bullet}=\sum_{c}N_{c}(N_{c}-1)$.
This calculation clarifies that a positive value of $Q$ implies that the
‘internal’ probability of connection $p_{\bullet}$ is larger than the one
‘predicted’ by the RGM, i.e.

$Q\geq 0\Longrightarrow
p_{\bullet}=\frac{l_{\bullet}^{*}}{V_{\bullet}}\geq\frac{L}{V}=p$ (51)

while a null modularity value implies that $p_{\bullet}$ matches the only
parameter defining the RGM. The result above can be also restated as follows:
a null modularity value points out that the mesoscale structure of the
configuration at hand _has not been generated_ by a SBM (or, equivalently,
_does not need_ a SBM to be explained). Notice, however, that the modularity
provides no indication about the statistical significance of the recovered
partition.

Figure 6: Comparison of three different algorithms for community detection on
the LFR benchmark, i.e. modularity maximization (left column), Infomap
(central column) and surprise minimization (right column): the chosen
configurations are binary undirected networks (first row), binary directed
networks (second row), weighted undirected networks (third row) and weighted
directed networks (fourth row). Four curves are shown, corresponding to two
different network sizes for binary configurations (blue and red: 1000 nodes;
green and black: 5000 nodes) and, for a given size, to two different ranges
for the size of communities, indicated by the letters S (‘small’, communities
have between 10 and 50 nodes) and B (‘big’, communities have between 20 and
100 nodes). For weighted configurations, instead, all networks have the same
size (5000 nodes): what differentiates the two pairs of curves is the value of
the binary mixing parameter. When focusing on binary, undirected networks, we
have considered $\tau_{1}=-2$ and $\tau_{2}=-1$ (the average degree is 20 and
the maximum degree is 50). When binary, directed configurations are
considered, $\mu_{t}$ refers to in-degrees, which are distributed according to
a power-law while the out-degrees are kept constant for all nodes; the other
input parameters, instead, are the same used for undirected configurations.
When weighted, undirected configurations are considered, an additional mixing
parameters is needed, i.e. $\mu_{w}$, accounting for the percentage of a node
strength to be distributed on the links that connect it to the nodes outside
its own community; the exponent of the strength distribution has been set to
1.5 for all realizations considered here. When weighted, directed
configurations are considered, $\mu_{w}$ refers to in-strengths. The trend of
the NMI, plotted as a function of the mixing parameters, reveals Infomap to be
a strong performer even if its performance decreases abruptly as the mixing
parameters exceed a threshold value that depends on the particular setting of
the benchmark; the performance of modularity, instead, ‘degrades’ less
rapidly. The performance of surprise constitutes a good compromise between the
robustness of modularity and the steadily high accuracy of Infomap.

Comparison with the weighted modularity. A similar conclusion can be reached
upon considering the modularity function for community detection on weighted,
undirected networks, defined as

${\color[rgb]{0,0,0}Q=\frac{1}{2W}\sum_{i\neq j}(w_{ij}-\langle
w_{ij}\rangle)\delta_{c_{i},c_{j}}}$ (52)

where we have employed the Weighted Random Graph Model (WRGM) as a benchmark,
according to which $\langle
w_{ij}\rangle_{\text{WRGM}}=\frac{W}{V}=\frac{2W}{N(N-1)}$, $\forall\>i<j$
Garlaschelli2009b . From the definition of weighted modularity provided above,
it follows that

$\displaystyle Q$ $\displaystyle=\frac{1}{2W}\sum_{i\neq
j}w_{ij}\delta_{c_{i},c_{j}}-\frac{1}{2W}\sum_{i\neq j}\langle
w_{ij}\rangle_{\text{WRGM}}\delta_{c_{i},c_{j}}$ (53)
$\displaystyle=\frac{1}{2W}\sum_{c}[2w_{c}^{*}-\langle
w_{ij}\rangle_{\text{WRGM}}N_{c}(N_{c}-1)]$ (54)
$\displaystyle=\frac{1}{2W}\sum_{c}\left[2w_{c}^{*}-2W\frac{N_{c}(N_{c}-1)}{N(N-1)}\right]$
(55)
$\displaystyle=\left(\frac{w_{\bullet}^{*}}{W}-\frac{V_{\bullet}}{V}\right)$
(56)

where $N_{c}$ and $w_{c}^{*}$ indicate the number of nodes of community $c$
and the weight of links _within_ community $c$, respectively - naturally,
$w_{\bullet}^{*}=\sum_{c}w_{c}^{*}$. As in the binary case, the calculation
above clarifies that a positive value of $Q$ implies that the expected weight
of any ‘internal’ link is larger than the one predicted by the WRGM, i.e.

$Q\geq 0\Longrightarrow\langle
w_{ij}\rangle_{\text{WSBM}}=\frac{w_{\bullet}^{*}}{V_{\bullet}}\geq\frac{W}{V}=\langle
w_{ij}\rangle_{\text{WRGM}};$ (57)

a null modularity value, on the other hand, implies that $\langle
w_{ij}\rangle_{\text{WSBM}}$ coincides with the expectation coming from the
WRGM. Analogously to the binary case, the result above can be also restated as
follows: a null modularity value points out that the mesoscale structure of
the configuration at hand _has not been generated_ by a WSBM (or,
equivalently, _does not need_ a WSBM to be explained). Again, however, the
modularity provides no indication about the statistical significance of the
recovered partition.

Figure 7: Comparison of three different algorithms for community detection on
the RC benchmark, i.e. modularity maximization, Infomap and surprise
minimization. The chosen configurations are binary, undirected (first row) and
directed (second row) networks, consisting of 512 nodes, grouped in 16
communities arranged in a ring-like fashion and whose sizes obey a power-law
whose exponent has been set to 1.8; the smallest community is composed by 3
nodes. Such initial configurations are progressively ‘degraded’ according to
the following mechanism: first, a percentage $p$ of links is selected randomly
and removed; afterwards, a percentage $p$ of links is selected randomly and
rewired. The trend of the NMI, plotted as a function of the (single)
‘degradation’ parameter $p$, driving the evolution of the initial ring-of-
cliques towards a progressively less-defined clustered configuration, reveals
surprise to outperform both modularity and Infomap. From a more general
perspective, these results confirm what has been already observed elsewhere,
i.e. that the best-performing algorithms on the LFR benchmarks often perform
poorly on the RC benchmarks and vice versa.

Figure 8: First, second and third row: benchmark for testing surprise on the
recovery of core-periphery structures (left column: binary, undirected
networks; right column: binary, directed networks). The one we have defined
here progressively ‘degrades’ an initial configuration, defined by 1) a
completely connected core, 2) an empty periphery, 3) an intermediate part
whose link density amounts at $p_{cp}=0.5$. Such a configuration is ‘degraded’
by progressively filling the periphery and emptying the core. This is achieved
by 1) considering all peripherical node pairs and link them with probability
$q$; 2) considering all core node pairs and keep them linked with probability
$1-q$: varying $q$ in the interval $[0,p_{cp}]$ allows us to span a range of
configurations starting with Borgatti-Everett and ending with Erdös-Rényi. As
expected, the performance of the surprise worsens as the ‘degradation’
parameter becomes closer to $p_{cp}=0.5$; however, both the NMI and the ARI
indices steadily remain very close to 1. Bottom row: the presence of weights
affects the detection of ‘bimodular’ mesoscale structures as well. In fact,
rising the weight of any two links connecting the core with the periphery of a
toy network allows the two nodes originally part of the periphery to be
detected as belonging to the core; analogously, if a bipartite topology is
modified by adding weights between some of the nodes belonging to the same
layer, a core-periphery structure will now be detected as significant.

Comparing mesoscale structures detection methods. The previous subsections
have been devoted to check the consistency of our surprise-based formalism and
reveal the ‘statistical flaws’ of other approaches. Let us now consider two
popular algorithms for mesoscale structures detection, i.e. modularity
maximization (now, $Q$ has been considered in its ‘full’ definition, e.g.
$\langle a_{ij}\rangle=p_{ij}=\frac{k_{i}k_{j}}{2L}$, $\forall\>i<j$ for
binary, undirected configurations) and Infomap Rosvall2008 , and carry out
more systematical comparisons between the former and the surprise. Upon doing
so, we are able to compare one algorithm per class, i.e. modularity for the
first class, surprise for the second class and Infomap for the third class.

To this aim, we have focused on different kinds of benchmarks, i.e. classes of
synthetic networks with well-defined planted partitions, the aim being that of
inspecting the goodness of a given algorithm in recovering the imposed
partition. As an indicator of the goodness of the partition retrieved by each
algorithm, we have followed Bongiorno2017  and employed three different
indices. The first one is the _normalized mutual information_ (NMI), defined
as

$\overline{I}(X,Y)=\frac{2I(X,Y)}{H(X)+H(Y)}$ (58)

where partitions $X$ and $Y$ are compared, $H(X)=-\sum_{x}f_{x}\ln f_{x}$ and
$f_{x}=\frac{n_{x}}{n}$ is the fraction of nodes assigned to the cluster
labeled with $x$; analogously, for the partition $Y$. The term
$I(X,Y)=\sum_{x}\sum_{y}f_{xy}\ln\left(\frac{f_{xy}}{f_{x}f_{y}}\right)$ is
the ‘proper’ _mutual information_ and $f_{xy}=\frac{n_{xy}}{n}$ is the
fraction of nodes assigned to cluster $x$ in partition $X$ and to cluster $y$
in partition $Y$. Naturally, $\overline{I}(X,Y)$ equals 1 if the partitions
are identical and 0 if the partitions are independent.

The second index we have considered is the _adjusted Rand index_ (ARI),
defined as

$\text{ARI}=\frac{TP+TN-\langle TP+TN\rangle}{TP+FP+TN+FN-\langle
TP+TN\rangle}$ (59)

and representing a sort of accuracy333The number of true positives (TP) is the
number of pairs of nodes being in the same community both in the considered
and in the reference partition; the number of false positives (FP) is the
number of pairs of nodes being in the same community in the considered
partition but in different communities in the reference partition; the number
of true negatives (TN) is the number of pairs of nodes being in the same
community neither in the considered nor in the reference partition; the number
of false negatives (FN) is the number of pairs of nodes being in the same
community in the reference partition but not in the considered partition.
‘corrected’ by a term that quantifies the agreement between the reference
partition and a random partition - the term ‘random’ referring to the
Permutation Model; equivalently, the closer the ARI to 0, the more ‘random’
the provided partition.

Figure 9: Application of our framework for the detection of binary and
weighted communities on a bunch of social networks: (top left) binary network
of co-occurrences between Star Wars characters repository0 , (top right)
weighted network of co-occurrences between Star Wars characters repository0 ,
(bottom left) weighted ‘friendship’ network among the terrorists involved in
the train bombing of Madrid in 2004 repository , (bottom right) weighted
‘friendship’ network among the residents living in an Australian University
Campus repository . Overall, link weights refine the picture provided by just
considering the presence of links. A first example is provided by the two
‘friendship’ networks where sparse subgraphs can be considered as communities
whenever their connections are ‘heavy’ enough; a second example is provided by
the weighted communities of Star Wars characters: the three major clusters, in
this case, are the one induced by the characters of Episodes I-III, the one
induced by the characters of Episodes IV-IX and the one contaning the villains
of Episodes VII-IX.

The third index we have considered is the _adjusted Wallace index_ (AWI),
defined as

$\text{AWI}=\frac{TP-\langle TP\rangle}{TP+FP-\langle TP\rangle}$ (60)

and representing a sort of ‘corrected’ positive predicted value. Again, the
closer the AWI to 0, the more ‘random’ the provided partition.

First, let us inspect the performance of modularity, surprise and Infomap to
detect cliques arranged in a ring. Specifically, we have considered 7,
different ring-like configurations, each one linking 20 binary cliques (i.e.
$K_{3}$, $K_{4}$, $K_{5}$, $K_{8}$, $K_{10}$, $K_{15}$, $K_{20}$). As fig. 4
reveals, surprise always recovers the planted partition; on the other hand,
modularity maximization leads to miss the partitions with $K_{3}$ and $K_{4}$
and Infomap misses the partition with $K_{3}$, a result that may be a
consequence of the resolution limit, affecting both the aforementioned
algorithms.

Let us now ask us if the presence of weights affects the detection of
mesoscale structures. Generally speaking, the answer is yes, as the comparison
between the ‘ring of binary cliques’ and the ‘ring of weighted cliques’ cases
shows. In particular, the result according to which surprise minimization is
able to discriminate the inter-linked cliques changes once the weight of the
links connecting any two cliques is risen: in fact, this leads the algorithm
to reveal as ‘communities’ two tightly-connected pairs of cliques, now. We
explicitly notice that the results shown in fig. 4 also depend on the relative
magnitude of the weights of the intra-cliques and of the inter-cliques links:
however, as long as the inter-cliques weight is up to two orders of magnitude
larger than the intra-cliques one, it holds true.

Figure 10: Top panels: results of the application of our framework for the
detection of ‘bimodular’ structures on the network of co-occurrences of the
‘Les Miserables’ characters repository  (first panel: binary version; second
panel: weighted version). Link weights refine the picture provided by just
considering the presence of links: the core (in black) is, in fact,
constituted by the nodes connected by the ‘heavier’ links, irrespectively from
the link density of the former one. Bottom panels: results of the application
of our framework for the detection of weighted ‘bimodular’ structures on the
electronic Italian Interbank Money Market (e-MID) Iori2006 , for the
maintenance period approximately corresponding to May 2009 (third panel:
binary version; fourth panel: weighted version). As the picture shows, the
(purely binary) core links are also the ‘heavier’ ones - an evidence
confirming an ubiquitous tendency in economic and financial systems, i.e.
binary and weighted quantities are closely related.

In order to expand the set of comparisons, we have focused on two different
kinds of well-established benchmarks, i.e. the Lancichinetti-Fortunato-
Radicchi (LFR) one and the Aldecoa’s ‘relaxed-caveman’ (RC) one Aldecoa2021;
Lancichinetti2009 .

The LFR benchmark is a special case of the planted-partition model, in which
groups have different sizes and nodes have different degrees - hence
constituting a refinement of the GN benchmark, where groups have equal size
and nodes have the same expected degree Lancichinetti2009 . When binary,
undirected configurations are considered, the degrees of nodes are distributed
according to a power-law with exponent $\tau_{1}$ while the sizes of
communities are distributed according to a power-law with exponent $\tau_{2}$.
Once the sizes of communities have been drawn each node ‘receives’ its own
degree, say $k_{i}$: the percentage of these links connecting node $i$ with
other internal nodes is, then, chosen to be $(1-\mu_{t})k_{i}$, with $\mu_{t}$
playing the role of mixing parameter that controls for the sharpness of the
planted partition. For binary, directed configurations, $\mu_{t}$ refers to
in-degrees, which are distributed according to a power-law while the out-
degrees are kept constant for all nodes; the other input parameters, instead,
are the same used for undirected configurations.

Results on specific implementations of the LFR benchmark are shown in fig. 6.
Infomap is, generally speaking, a strong performer; as evident upon looking at
the first row, however, its performance decreases abruptly as the mixing
parameter exceeds a threshold value that depends on the particular setting of
the LFR benchmark. Modularity, instead, seems to be more robust (i.e. its
performance ‘degrades’ less rapidly as $\mu_{t}$ increases) although the
resolution limit manifests itself when configurations with small communities
are considered. Overall, the performance of surprise seems to constitute a
good compromise between the robustness of modularity and the steadily high
accuracy of Infomap. We would also like to stress that surprise competes with
modularity although it employes much less information than the latter: in
fact, while the benchmark employed by modularity coincides with the (sparse
version of the) Configuration Model - hence, encodes the information on the
entire degree sequence - surprise compares the RGM with the SBM, hence
employing the information on the link density, both in a global and in a
block-wise fashion.

Surprise becomes the best performer when binary, directed configurations are
considered - see the second row of fig. 6: while the performance of modularity
starts decreasing as soon as the value of $\mu_{t}$ is risen and
$\text{NMI}_{Infomap}\simeq 0$ when $\mu_{t}$ crosses the value of $0.6$, the
performance of surprise ‘degrades’ much more slowly - in fact, for some
instances of the LFR benchmark, it achieves a large value of NMI even for
values $\mu_{t}\geq 0.8$.

Let us now comment on the performance of our algorithms when weighted
configurations are considered. The LFR benchmark, in fact, can be extended to
account for weights as well. In this case, there are two mixing parameters:
while the first one is the usual, topological mixing parameter, the second
one, that will be indicated with $\mu_{w}$, accounts for the heterogeneity of
weights - more precisely, controlling for the percentage of a node strength to
be distributed on the links that connect it to the nodes outside its own
community. The third parameter to be determined is the exponent of the
strength distribution - that has been set to 1.5 for all realizations
considered here. The results are, again, shown in fig. 6 (to be noticed that
we have kept one of the two parameters fixed and studied the dependence of NMI
and ARI on the other: specifically, we have frozen the topological mixing
parameter and studied the dependence of the results on $\mu_{w}$, thus
inspecting the performance of our algorithms as the weights are redistributed
on a fixed topology).

Infomap is, again, a strong performer although its performance keeps
decreasing abruptly as $\mu_{w}$ exceeds a threshold value depending on the
particular setting of the LFR benchmark; modularity, instead, performs worse
than in the binary case although it is still more robust than Infomap.
Although ‘degrading’ less sharply than Infomap, the performance of the purely
weighted surprise seems to be the worst, here; on the other hand, the enhanced
surprise outperforms the competing algorithms for intermediate values of the
topological mixing parameter, irrespectively from the size of the communities:
in fact, $\text{NMI}_{surprise}=1$ even for $\mu_{w}=0.8$. Similar
considerations hold true when weighted, directed configurations are considered
(with the only difference that, now, modularity steadily performs worse than
the other algorithms, except for the largest values of the mixing parameter).
As for the binary cases, surprise competes with modularity although it
employes much less information than the latter: in fact, while the benchmark
employed by modularity now coincides with the (sparse version of the) Weighted
Configuration Model - hence, encodes the information on the entire strength
sequence - surprise compares the WRGM with the WSBM, hence employing the
information on the magnitude of the total weight, both in a global and in a
block-wise fashion.

Figure 11: Result of the application of our framework for the detection of
weighted ‘bimodular’ structures on the WTW in the year 2000 Gleditsch2002 .
Top panels: upon running $\mathscr{W}_{\sslash}$ in a hierarchical fashion, a
‘core-within-the-core’ is detected, revealing that the richest world countries
(i.e. Canada, USA, the richest European countries, China and Russia - right
panel, in dark green) constitute an even more tightly-linked cluster of
nations among the ones with the largest trading activity (left panel, in
black). Bottom panels: $\mathscr{E}_{\sslash}$ penalizes the countries with
both a small degree and a small strength; in fact, the second run exlcudes the
Russia, a result seemingly indicating that while its strength is ‘large
enough’ to be a member of the core, its degree is not.

Let us now consider the RC benchmark. It consists of 512 nodes, grouped in 16
communities arranged in a ring-like fashion and whose sizes obey a power-law
whose exponent has been set to 1.8; the smallest community is composed by 3
nodes. Such a configuration is progressively ‘degraded’ according to the
following mechanism: first, a percentage $p$ of links is selected randomly and
removed; afterwards, a percentage $p$ of links is selected randomly and
rewired. In other words, a single ‘degradation’ parameter $p$ drives the
evolution of the initial ring-of-cliques towards a progressively less-defined
clustered configuration.

Results on specific implementations of the RC benchmark are shown in fig. 7:
surprise outperforms both competing algorithms across the entire domain of the
‘degradation’ parameter $p$. More specifically, while modularity ‘degrades’
slowly as the value of $p$ is risen, Infomap ‘degrades’ abruptly as $p\geq
0.4$. Hence, for small values of such a parameter, Infomap outperforms
modularity; on the other hand, for large values of $p$, modularity outperforms
Infomap (although both $\text{NMI}_{modularity}$ and $\text{ARI}_{modularity}$
achieve a value which is around $0.6$, i.e. already far from the maximum).
Interestingly, for small values of $p$, the performance of Infomap and that of
surprise overlap, both achieving NMI and ARI values which are very close to 1:
as $p$ crosses the value of $0.4$, however, the two trends become increasingly
different with Infomap being outperformed by modularity which is, in turn,
outperformed by surprise. From a more general perspective, these results
confirm what has been already observed elsewhere Aldecoa2013 , i.e. that the
best-performing algorithms on the LFR benchmarks often perform poorly on the
RC benchmarks and vice versa.

Let us now inspect the performance of surprise in recovering binary
‘bimodular’ structures. To this aim, we have defined a novel benchmark
mimicking the philosophy of the RC one, i.e. progressively ‘degrading’ an
initial, well-defined configuration:

  * •

let us consider $N_{c}$ core nodes and $N_{p}$ periphery nodes. The core is
completely connected (i.e. the link density of the $N_{c}\times N_{c}$ block
is 1) and the periphery is empty (i.e. the link density of the $N_{p}\times
N_{p}$ block is 0). So far, our benchmark is reminiscent of a core-periphery
structure _à la Borgatti-Everett_ ;

  * •

let us now focus on the topology of the $N_{c}\times N_{p}$ bipartite network
embodying the connections between the core and the periphery: in particular,
let us consider each entry of such an adjacency matrix and pose $a_{cp}=1$
with probability $p_{cp}$. Upon doing so, such a subgraph will have a link
density amounting precisely at $p_{cp}$;

  * •

let us now ‘degrade’ such an initial configuration, by progressively filling
the periphery and emptying the core. This can be achieved by 1) considering
all peripherical node pairs and link them with probability $q$; 2) considering
all core node pairs and keep them linked with probability $1-q$ (or,
equivalently, disconnect them with probability $q$). Upon doing so, we end up
with a core whose link density is precisely $1-q$ and with a periphery whose
link density is precisely $q$. Now, varying $q$ in the interval $[0,p_{cp}]$
allows us to span a range of configurations starting with Borgatti-Everett and
ending with Erdös-Rényi.

Binary networks |  $\\#$ communities |
---|---|---
| Modularity | Infomap | Surprise | $\text{NMI}_{\text{SM}}$ | $\text{NMI}_{\text{SI}}$ | $\text{ARI}_{\text{SM}}$ | $\text{ARI}_{\text{SI}}$ | $\text{AWI}_{\text{SM}}$ | $\text{AWI}_{\text{SI}}$
Madrid bombing terrorists | 6 | 5 | 26 | 0.50 | 0.44 | 0.33 | 0.27 | 0.83 | 0.73
Star Wars characters | 9 | 5 | 35 | 0.37 | 0.53 | 0.15 | 0.28 | 0.39 | 0.94
Australian University Campus | 13 | 6 | 30 | 0.47 | 0.72 | 0.20 | 0.46 | 0.76 | 0.93
‘Les Miserables’ characters | 8 | 6 | 33 | 0.56 | 0.56 | 0.46 | 0.50 | 0.92 | 0.83
Weighted networks |  $\\#$ communities |
| Modularity | Infomap | Surprise | $\text{NMI}_{\text{SM}}$ | $\text{NMI}_{\text{SI}}$ | $\text{ARI}_{\text{SM}}$ | $\text{ARI}_{\text{SI}}$ | $\text{AWI}_{\text{SM}}$ | $\text{AWI}_{\text{SI}}$
Madrid bombing terrorists | 8 | 5 | 22 | 0.50 | 0.67 | 0.37 | 0.62 | 0.65 | 0.94
Star Wars characters | 10 | 6 | 52 | 0.36 | 0.43 | 0.18 | 0.24 | 0.79 | 0.79
Australian University Campus | 14 | 6 | 25 | 0.47 | 0.74 | 0.22 | 0.46 | 0.72 | 0.91
‘Les Miserables’ characters | 10 | 6 | 32 | 0.51 | 0.61 | 0.37 | 0.49 | 0.72 | 0.74
Table 2: Table comparing the performance of modularity, Infomap and surprise
on the bunch of real-world networks considered here. While Infomap is the
algorithm producing the smallest number of clusters, surprise is the one
producing the largest number of clusters; as a consequence, the three
partitions output by our algorithms have an overall small overlap.
Interestingly enough, the values of the AWI are quite large: since it just
focuses on the percentage of true positives, a good performance under such an
index indicates that the tested algorithms ‘agree’ on the nodes to be
clustered together: hence, the discrepancy between the ARI and the AWI may be
explained by the presence of statistical ‘noise’ (i.e. pairs of nodes on whose
classification the algorithms ‘disagree’) around the bulk of nodes to be put
together - on which the algorithms ‘agree’.

Specifically, here we have considered $N_{c}=100$, $N_{p}=300$ and
$p_{cp}=0.5$. The result of our exercise is shown in fig. 8: as expected, the
performance of the surprise worsens as the ‘degradation’ parameter becomes
closer to $p_{cp}=0.5$; however, both the NMI and the ARI indices steadily
remain very close to 1 - a result meaning that surprise optimization is not
only able to correctly classify true positives (i.e. to keep the nodes
originally in the same communities together) but also the other, possible
kinds of node pairs. As with the exercise on community detection, let us now
ask us if the presence of weights affects the detection of mesoscale
structures. Generally speaking, the answer is, again, yes. Let us consider a
toy core-periphery network: rising the weight of any two links connecting the
core with the periphery allows the two nodes originally part of the periphery
to be detected as belonging to the core. Analogously, if a bipartite topology
is modified by adding weights between some of the nodes belonging to the same
layer, $\mathscr{W}_{\sslash}$ will detect a core-periphery structure as
significant, the core nodes being the ones linked by the ‘heaviest’
connections.

Testing surprise on real-world networks. Let us now apply our formalism to the
detection of mesoscale structures in real-world networks.

When coming to study real systems, particularly insightful examples are
provided by social networks. To this aim, let us consider the one induced by
the co-occurrences of characters within the Star Wars saga (i.e. the three
trilogies) repository0 . As shown in fig. 9 we have both considered the binary
and the weighted version of it. In both cases two major clusters are visible.
For what concerns the binary version of such a network, the optimization of
$\mathscr{S}$ reveals the presence of two major clusters: remarkably, those
clusters are induced by the characters of Episodes I-III (e.g. Yoda, Qui-Gon,
Obi-Wan, Anakin, Padme, the Emperor, Count Dooku, etc.) and by the characters
of Episodes IV-IX (e.g. C-3PO, Leia, Han, Lando, Poe, Finn); a third cluster,
instead, concerns the villains of Episodes VII-IX (i.e. Snoke, Kylo Ren,
Phasma, Hux). Interestingly, Rey, BB-8, Maz-Kanata and other characters living
on Jakku are clustered together. Quite remarkably, the interactions between
the characters of Episodes IV-VI and those of Episodes VII-IX causes the
former ones and the latter ones to be recovered within the same cluster. This
picture is further refined once weights are taken into account: in fact, two
of the aforementioned clusters are now merged, giving origin to the cluster of
heroes of Episodes IV-IX (e.g. C-3PO, Leia, Han, Lando, Poe, Finn, Rey, Maz-
Kanata).

Let us now inspect the effectiveness of our framework in revealing weighted
communities by considering the ‘friendship’ networks among the terrorists
involved in the train bombing of Madrid in 2004 repository  and the one among
the residents living in an Australian University Campus repository  (see fig.
9). As the optimization of $\mathscr{W}$ reveals, while fully connected
subsets of nodes are considered as communities in case links have unitary
weights, sparser subgraphs can be considered as communities as well whenever
their inner connections are ‘heavy’ enough. On the other hand, both bottom
panels of fig. 9 seem to confirm that one of the main limitations of surprise-
like functionals is that of recovering a large number of small cluster of
nodes.

Let us now compare the performance of $\mathscr{S}_{\sslash}$ and
$\mathscr{W}_{\sslash}$ in order to see if, and how, the presence of weights
affects the ‘bimodular’ mesoscale organization of networks. To this aim, let
us focus on the network of co-occurrences of the characters of the novel ‘Les
Miserables’ repository . As shown in fig. 10, link weights indeed modify the
picture provided by just considering the simple presence of links (see also
DeJeude2019 ): the core of the weighted network is, in fact, constituted by
the nodes connected by the ‘heavier’ links, irrespectively from the link
density of the former one.

After having applied our framework to the analysis of social networks, let us
move to consider financial networks. One of the most popular examples of the
kind is provided by the electronic Italian Interbank Money Market (e-MID)
Iori2006 , depicted in fig. 10. Notice that, for such a network, the vast
majority of core links are also the ‘heavier’ ones, an evidence confirming a
tendency that is ubiquitous in financial and economic systems, i.e. binary and
weighted quantities - even at the mesoscale - are closely related.

Remarkably, the surprise-based formalism presented in this paper can be
employed in a hierarchical fashion to highlight either nested communities or
nested ‘bimodular’ structures. To clarify this point, let us consider the
World Trade Web (WTW) in the year 2000 as a case-study Gleditsch2002 . First,
let us run $\mathscr{W}_{\sslash}$ to highlight the core portion of the
weighted version of such a network; as fig. 11 shows, the bipartition
distinguishes countries with a large strength from those whose trade volume is
low (basically, a bunch of African, Asian and South-American countries).
Repeating our analysis _within_ the core portion of the network allows us to
discover the presence of a (statistically-significant) nested core: in fact,
the secound-round optimization reveals that the ‘core-inside-the-core’ is
composed by countries such as Canada, USA, the richest European countries,
China and Russia.

Let us now compare it with $\mathscr{E}_{\sslash}$, run in a hierarchical
fashion as well. The results of our exercise are shown in fig. 11. As evident
from looking at it, the enhanced surprise is more restrictive than the purely
weighted one, as a consequence of constraining the degrees beside the
strengths. Hence, while the first run excludes the countries with both a small
degree and a small strength, the second run exlcudes the Russia, a result
seemingly indicating that while its strength is ‘large enough’ to allow it to
be a member of the core, its degree is not.

In a sense, the optimization of $\mathscr{E}_{\sslash}$ corrects the picture
provided by the optimization of $\mathscr{W}_{\sslash}$ as the core becomes
less populated by ‘low degree’ nodes - an effect which is likely to become
more evident on systems that are neither financial nor economic in nature.

Let us now run, and compare, modularity, Infomap and surprise on the bunch of
real-world networks above. Table 2 sums up the results. A first observation
concerns the number of detected communities: while Infomap is the algorithm
producing the smallest number of clusters, surprise is the one producing the
largest number of clusters - more precisely, surprise outputs _more_ and
_smaller_ clusters than the other two methods. As a consequence, our three
algorithms produce partitions with an overall small overlap, as indicated by
the NMI; the ARI confirms such an observation - although indicating that the
pictures provided by Infomap and surprise are (overall) more similar than
those provided by modularity and surprise (and, as a consequence, by
modularity and Infomap). Interestingly enough, the values of the AWI are quite
large - and larger than the corresponding NMI and ARI values: since it just
focuses on the percentage of true positives, a good performance under such an
index indicates that the two tested algorithms ‘agree’ on the nodes to be
clustered together (although they may not - and, in general, will not - agree
on the number of communities). Hence, the discrepancy between the ARI and the
AWI may be explained by the presence of statistical ‘noise’ (i.e.
‘misclassified’ pairs of nodes, although the word may not be correct as the
information about the ‘true’ partition is not available) around the bulk of
nodes to be put together.

Let us stress once more that, whenever real-world networks are considered,
information about the existence of a ‘true’ partition is rarely available; for
this reason, exercises as the one we have carried out here may be useful to
gain insight on the system under study: instead of trusting just one
algorithm, combining pairs of them - e.g. by considering as communities the
subsets of nodes output by both - may be the right solution to overcome the
limitations affecting each single method.

The ‘SurpriseMeMore’ Python package. As an additional result, we release a
comprehensive package, coded in Python, that implements all aforementioned
surprise-like variants. Its name is ‘SurpriseMeMore’ - the name recalls the
package released in 2014 Aldecoa2014  \- and is freely downloadable at the
following URL: https://github.com/EmilianoMarchese/SurpriseMeMore.

## CONCLUSIONS

The hypergeometric distribution - together with its many variants - has
recently revived the interest of researchers who have employed it to define
novel network ensembles Casiraghi2018 , recipes for projecting bipartite
networks Cimini2021b , etc.

Remarkably, the distributions related to it allow for a wide variety of
benchmarks to be defined, each one embodying a different set of constraints.
In the present paper we have explored the ‘power’ of the hypergeometric
formalism to carry out the detection of mesoscale structures: remarkably, it
allows proper statistical tests to be definable for revealing the presence of
modular, core-periphery and bipartite structures on any kind of network, be it
binary or weighted, undirected or directed. According to the classification
proposed in the introductory section of the paper, we believe surprise to
belong to the second class of algorithms - its asymptotic expression embodying
a sort of LRT aimed at choosing between alternative hypotheses (the RGM and
the SBM, the WRG and the WSBM, etc.).

More in general, our approach reveals the superiority of the algorithms for
mesoscale structures detection belonging to the second and to the third class
with respect to those belonging to the first one: still, the two classes of
statistically-grounded approaches compete on some, specific benchmarks, as the
comparisons carried out on the LFR and the RC ones clearly show (specifically,
methods performing well on LFR benchmarks do not on RC benchmarks and vice
versa). This also suggests a strategy to handle mesoscale structures detection
on real-world networks: as the information on the presence of a possible,
‘true’ partition is rarely available, a good strategy may be that of running
different algorithms on the same empirical networks, check the consistency of
their output and combine them, e.g. by taking the overlap - in a way that is
reminiscent of multi-model inference.

Although the surprise-based approach is powerful and versatile, its downside
is represented by the specific kinds of tests that are induced by the
optimization of surprise-like score functions: comparisons between benchmarks
‘ignoring’ the local structure of nodes (i.e. their degree) are, in fact,
carried out. While this seems to be perfectly reasonable when considering
core-periphery structures - see also the contribution Kojaku2018  whose
authors claim that a core-periphery structure is always compatible with a
network degree sequence - this is no longer true for the community detection
task Karrer2010  \- indeed, as it has been noticed elsewhere, ignoring degrees
may be at the origin of the large number of singletons output by surprise as
assigning nodes with few neighbors to larger clusters may be disfavored from a
statistical point of view.

The observations above call for the ‘extension’ of the hypergeometric
formalism we have studied here to include more refined benchmarks as the ones
constraining the entire degree sequence.

## DATA AVAILABILITY

Network data that support the findings of this study are freely downloadable
at the following URLs: https://github.com/evelinag/StarWars-social-network,
http://konect.cc/networks, http://ksgleditsch.com/exptradegdp.html.

## CODE AVAILABILITY

The code employed for the analyses carried out in this study is freely
downloadable at the following URL:
https://github.com/EmilianoMarchese/SurpriseMeMore.

## COMPETING INTERESTS

The authors declare no competing financial interests.

## ACKNOWLEDGEMENTS

The authors acknowledge support from the EU project SoBigData-PlusPlus (grant
no. 871042).

## References

  * (1) S. Fortunato, D. Hric, Community detection in networks: a user guide, _Phys. Rep._ 659, 1-44 (2016).
  * (2) B. S. Khan, M. A. Niazi, Network community detection: a review and visual survey, arXiv:1708.00977 (2017).
  * (3) S. P. Borgatti, M. G. Everett, Models of core/periphery structures, _Soc. Networks_ 21 (4), 375-395 (2000).
  * (4) B. Craig, G. von Peter, Interbank tiering and money center banks, _J. Finan. Intermediation_ 23 (3), 322-347 (2014).
  * (5) D. in ’t Veld, I. van Lelyveld, Finding the core: Network structure in interbank markets, _J. Bank. Finance_ 49, 27-40 (2014).
  * (6) D. T. Luu et al., Collateral unchained: Rehypothecation networks, concentration and systemic effects, _Sciences Po OFCE working paper n. 07, 2018/01/31_ (2018).
  * (7) S. Fortunato, Community detection in graphs, _Phys. Rep._ 486, 75-174 (2010).
  * (8) G. Cimini, R. Mastrandrea, T. Squartini, Reconstructing networks, in Elements in Structure and Dynamics of Complex Networks, Cambridge University Press (2021).
  * (9) A. Fronczak, Exponential random graph models, in Encyclopedia of Social Network Analysis and Mining, R. Alhajj, J. Rokne (Eds.), Springer-Verlag (2014).
  * (10) T. Peixoto, Descriptive vs. inferential community detection: pitfalls, myths and half-truths, arXiv:2112.00183 (2021).
  * (11) B. Karrer, M. E. J. Newman, Stochastic blockmodels and community structure in networks, _Phys. Rev. E_ 83, 016107 (2011).
  * (12) J. van Lidth de Jeude, G. Caldarelli, T. Squartini, Detecting core-periphery structures by surprise, _Europhys. Lett._ 125 (68001) (2018).
  * (13) T. Peixoto, Bayesian stochastic blockmodeling, arXiv:1705.10225 (2017).
  * (14) R. Aldecoa, I. Marin, Surprise maximization reveals the community structure of complex networks, _Sci. Rep._ 3 (1060) (2013).
  * (15) C. Nicolini, A. Bifone, Modular structure of brain functional networks: breaking the resolution limit by Surprise, _Sci. Rep._ 6 (19250) (2016).
  * (16) M. Tumminello, S. Micciché, F. Lillo, J. Piilo, R. N. Mantegna, Statistically validated networks in bipartite complex systems, _PLoS ONE_ 6 (3), e17994 (2011).
  * (17) C. Bongiorno, A. London, S. Micciché, R. N. Mantegna, Core of communities in bipartite networks, _Phys. Rev. E_ 96 (022321) (2017).
  * (18) F. Musciotto, F. Battiston, R. N. Mantegna, Detecting informative higher-order interactions in statistically-validated hypergraphs, arXiv:2103.16484 (2021).
  * (19) S. Micciché. R. N. Mantegna, A primer on statistically-validated networks, _Proceedings of the International School of Physics ‘Enrico Fermi’_ 203 (2019).
  * (20) Y. Jiang, C. Jia, J. Yu, An efficient community detection algorithm using greedy surprise maximization, _J. Phys. A: Math. Theor._ 47 (165101) (2014).
  * (21) J. Del Ser, J. L. Lobo, E. Villar-Rodriguez, M. N. Bilbao and C. Perfecto, Community detection in graphs based on surprise maximization using firefly heuristics, 2016 IEEE Congress on Evolutionary Computation (CEC), 2233-2239 (2016).
  * (22) Y. -N. Tang et al., An effective algorithm for optimizing surprise in network community detection, in IEEE Access 7, 148814-148827 (2019).
  * (23) D. Gamermann, J. A. Pellizaro, An algorithm for network community structure determination by surprise, arXiv:2012.13780 (2020).
  * (24) S. Kojaku, N. Masuda, A generalised significance test for individual communities in networks, _Sci. Rep._ 8 (7351) (2018).
  * (25) V. A. Traag, R. Aldecoa, J.-C. Delvenne, Detecting communities using asymptotical surprise, _Phys. Rev. E_ 92 (022816) (2015).
  * (26) B. Karrer, M. E. J. Newman, Stochastic blockmodels and community structure in networks, _Phys. Rev. E_ 83 (016107) (2011).
  * (27) X. Zhang, T. Martin, M. E. J. Newman, Identification of core-periphery structure in networks, _Phys. Rev. E_ 91 (032803) (2014).
  * (28) P. Barucca, F. Lillo, Disentangling bipartite and core-periphery structure in financial networks, arXiv:1511.08830v1 (2015).
  * (29) S. Kojaku, N. Masuda, Core-periphery structure requires something else in the network, _New J. Physics_ 4 (20), 359-395 (2018).
  * (30) P. Holme et al., Network bipartivity, _Phys. Rev. E_ 68 (5) (2003).
  * (31) E. Estrada, J. A. Rodríguez-Velázquez, Spectral measures of bipartivity in complex networks, _Phys. Rev. E_ 72 (4), 1-16 (2005).
  * (32) W. Feller, An introduction to probability theory and its applications, 1 (2nd ed.), Wiley (1950).
  * (33) D. Garlaschelli, M. I. Loffredo, Generalized Bose-Fermi statistics and structural correlations in weighted networks, _Phys. Rev. Lett._ 102 (038701) (2009).
  * (34) R. Mastrandrea, T. Squartini, G. Fagiolo, D. Garlaschelli, Enhanced reconstruction of weighted networks from strengths and degrees, _New J. Phys._ 16 (043022) (2014).
  * (35) D. Garlaschelli, The weighted random graph model, _New J. Phys._ 11 (073005) (2009).
  * (36) M. Rosvall, C. T. Bergstrom, Maps of random walks on complex networks reveal community structure, _PNAS_ 105, 1118-1123 (2008).
  * (37) R. Aldecoa, I. Marin, Closed benchmarks for network community structure characterization, _Phys. Rev. E_ 85 (026109) (2012).
  * (38) A. Lancichinetti, S. Fortunato, Community detection algorithms: a comparative analysis, _Phys. Rev. E_ 80 (056117) (2009).
  * (39) R. Aldecoa, I. Marin, Surprise maximization reveals the community structure of complex networks, _Sci. Rep._ 3 (1060) (2013).
  * (40) https://github.com/evelinag/StarWars-social-network
  * (41) http://konect.uni-koblenz.de/
  * (42) G. De Masi, G. Iori, G. Caldarelli, Fitness model for the Italian interbank money market, _Phys. Rev. E_ 74 (066112) (2006).
  * (43) K. S. Gleditsch, Expanded trade and GDP data, _J. Conflict Resol._ 46, 712-724 (2002).
  * (44) R. Aldecoa, I. Marin, SurpriseMe: an integrated tool for network community structure characterization using Surprise maximization, _Bioinformatics_ 30 (7), 1041-1042 (2014).
  * (45) G. Casiraghi, V. Nanumyan, Generalised hypergeometric ensembles of random graphs: the configuration model as an urn problem, arXiv:1810.06495 (2018).
  * (46) G. Cimini, A. Carra, L. Didomenicantonio, A. Zaccaria, Meta-validation of bipartite network projections, arXiv:2105.03391 (2021).

##  Appendix A Detecting mesoscale structures in binary networks

As stressed in the main text, community detection on binary networks can be
implemented by carrying out an exact statistical test whose function $f$ must
be identified with a binomial hypergeometric distribution. Hence, we can
compactly write

$f(l_{\bullet})\equiv\text{H}(l_{\bullet}|V,V_{\bullet},L)=\frac{\prod_{i=\bullet,\circ}\binom{V_{i}}{l_{i}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}$
(61)

with a clear meaning of the symbols. The correctedness of such the definition
of can be explicitly checked upon calculating its normalization via the
Vandermonde’s identity
$\sum_{k=0}^{r}\binom{m}{k}\binom{n}{r-k}=\binom{m+n}{r}$. While the binomial
coefficient $\binom{V_{\bullet}}{l_{\bullet}}$ enumerates the number of ways
$l_{\bullet}$ links can be redistributed within communities, i.e. over the
available $V_{\bullet}$ node pairs, the binomial coefficient
$\binom{V-V_{\bullet}}{L-l_{\bullet}}$ enumerates the number of ways the
remaining $L-l_{\bullet}$ links can be redistributed between communities, i.e.
over the remaining $V-V_{\bullet}$ node pairs. The position above, in turn,
induces the definition of the binary surprise, i.e.
$\mathscr{S}\equiv\sum_{l_{\bullet}\geq l_{\bullet}^{*}}f(l_{\bullet})$.

The generalization of the formalism above to detect ‘bimodular’ structure on
binary networks is straightforward. It is enough to consider the multivariate
(or multinomial) version of the hypergeometric distribution, i.e.

$f(l_{\bullet},l_{\circ})\equiv\text{MH}(l_{\bullet},l_{\circ}|V,V_{\bullet},V_{\circ},L)=\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}}{l_{i}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V_{\top}}{l_{\top}}}{\binom{V}{L}}=\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V-(V_{\bullet}+V_{\circ})}{L-(l_{\bullet}+l_{\circ})}}{\binom{V}{L}};$
(62)

the correctedness of such a definition of can be checked upon calculating its
normalization by applying the identity
$\sum_{k=0}^{r}\binom{m}{k}\binom{n}{r-k}=\binom{m+n}{r}$ twice. The position
above, in turn, induces the definition of the binary bimodular surprise, i.e.
$\mathscr{S}_{\sslash}\equiv\sum_{l_{\bullet}\geq
l^{*}_{\bullet}}\sum_{l_{\circ}\geq l^{*}_{\circ}}f(l_{\bullet},l_{\circ})$.

The asymptotic expression of $\mathscr{S}$ and $\mathscr{S}_{\sslash}$ can be
derived by applying the Stirling approximation $n!\simeq\sqrt{2\pi
n}\left(\frac{n}{e}\right)^{n}$ to the binomial coefficients appearing into
the definition of the versions of the surprise described above, leading to
expressions like

$\displaystyle\binom{V}{L}$ $\displaystyle=$
$\displaystyle\frac{V!}{L!(V-L)!}\simeq\frac{\sqrt{2\pi
V}V^{V}e^{-V}}{\sqrt{2\pi L}L^{L}e^{-L}\sqrt{2\pi(V-L)}(V-L)^{V-L}e^{-(V-L)}}$
(63) $\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2\pi
V\frac{L}{V}\left(1-\frac{L}{V}\right)}}\cdot\frac{V^{V}}{L^{L}(V-L)^{V-L}}$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2\pi
Vp(1-p)}}\cdot\frac{1}{\left(\frac{L}{V}\right)^{L}\left(1-\frac{L}{V}\right)^{V-L}}$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2\pi
Vp(1-p)}}\cdot\frac{1}{p^{L}(1-p)^{V-L}}=\frac{1}{\sqrt{2\pi\sigma^{2}}\cdot\text{Ber}(V,L,p)}$

where $\text{Ber}(x,y,z)=z^{y}(1-z)^{x-y}$ is a Bernoulli probability mass
function, $z=\frac{y}{x}$ and $p=\frac{L}{V}$.

The time, in seconds, to find the corresponding optimal partitions is reported
in tables 3 and 4: overall, denser networks require more time.

Network |  # nodes $(N)$ |  # links $(L)$ | Connectance |  Binary (s) |  Weighted (s) |  Enhanced (s)
---|---|---|---|---|---|---
Madrid terrorists | 64 | 486 | 0.24 | 0.05 | 0.05 | 0.07
Star Wars | 110 | 398 | 0.07 | 0.13 | 0.14 | 0.24
Australian college | 217 | 2672 | 0.11 | 0.34 | 0.30 | 0.38
Table 3: Time, in seconds, employed by $\mathscr{S}$, $\mathscr{W}$ and
$\mathscr{E}$ to find the corresponding optimal partition.

##  Appendix B Detecting mesoscale structures in weighted networks

When analyzing weighted networks, one needs to move from the binomial
hypergeometric distribution to its negative counterpart, i.e. to

$f(w_{\bullet})=\text{NH}(w_{\bullet}|V+W,W,V_{\bullet})=\frac{\prod_{i=\bullet,\circ}\binom{V_{i}+w_{i}-1}{w_{i}}}{\binom{V+W-1}{W}}=\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}}{\binom{V+W-1}{W}}=\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V-V_{\bullet}+(W-w_{\bullet})-1}{W-w_{\bullet}}}{\binom{V+W-1}{W}}$
(64)

with a clear meaning of the symbols. The soundness of the definition above can
be cheked upon calculating the normalization of the distribution
$\text{NH}(w_{\bullet}|V+W,W,V_{\bullet})$ via the identity
$\sum_{k=0}^{r}\binom{k+m}{k}\binom{n-m-k}{r-k}=\binom{n+1}{r}$: in fact

$\displaystyle\sum_{w_{\bullet}=0}^{W}\frac{\binom{w_{\bullet}+V_{\bullet}-1}{w_{\bullet}}\binom{w_{\circ}+V_{\circ}-1}{w_{\circ}}}{\binom{V+W-1}{W}}=\sum_{w_{\bullet}=0}^{W}\frac{\binom{w_{\bullet}+V_{\bullet}-1}{w_{\bullet}}\binom{W-w_{\bullet}+V-V_{\bullet}-1}{W-w\bullet}}{\binom{V+W-1}{W}}=\sum_{w_{\bullet}=0}^{W}\frac{\binom{w_{\bullet}+V_{\bullet}-1}{w_{\bullet}}\binom{(W+V-2)-(V_{\bullet}-1)-w_{\bullet}}{W-w_{\bullet}}}{\binom{V+W-1}{W}}=1.$
(65)

The position above induces the definition of the weighted surprise as
$\mathscr{W}=\sum_{w_{\bullet}\geq w_{\bullet}^{*}}f(w_{\bullet})$.

The generalization of the formalism above to detect ‘bimodular’ structure on
weighted networks is straightforward. It is enough to consider the
multivariate (or multinomial) version of the negative hypergeometric
distribution, i.e.

$\displaystyle
f(w_{\bullet},w_{\circ})=\text{MNH}(w_{\bullet},w_{\circ}|V+W,W,V_{\bullet},V_{\circ})$
$\displaystyle=$
$\displaystyle\frac{\prod_{i=\bullet,\circ,\top}\binom{V_{i}+w_{i}-1}{w_{i}}}{\binom{V+W-1}{W}}=\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}\binom{V_{\top}+w_{\top}-1}{w_{\top}}}{\binom{V+W-1}{W}}$
(66) $\displaystyle=$
$\displaystyle\frac{\binom{V_{\bullet}+w_{\bullet}-1}{w_{\bullet}}\binom{V_{\circ}+w_{\circ}-1}{w_{\circ}}\binom{V-(V_{\bullet}+V_{\circ})+W-(w_{\bullet}+w_{\circ})-1}{W-(w_{\bullet}+w_{\circ})}}{\binom{V+W-1}{W}}$

whose correctedness can be checked upon calculating its normalization, by
applying the identity
$\sum_{k=0}^{r}\binom{k+m}{k}\binom{n-m-k}{r-k}=\binom{n+1}{r}$ twice. The
position above, in turn, induces the definition of the weighted bimodular
surprise, as $\mathscr{W}_{\sslash}=\sum_{w_{\bullet}\geq
w^{*}_{\bullet}}\sum_{w_{\circ}\geq w^{*}_{\circ}}f(w_{\bullet},w_{\circ})$.

Network |  # nodes $(N)$ |  # links $(L)$ | Connectance |  Binary (s) |  Weighted (s) |  Enhanced (s)
---|---|---|---|---|---|---
‘Les Miserables’ | 77 | 254 | 0.09 | 0.22 | 0.31 | 0.48
e-MID 2009 | 134 | 4609 | 0.26 | 1.42 | 0.80 | 4.12
WTW 2000 | 187 | 20105 | 0.58 | 12.16 | 12.05 | 21.23
Table 4: Time, in seconds, employed by $\mathscr{S}_{\sslash}$,
$\mathscr{W}_{\sslash}$ and $\mathscr{E}_{\sslash}$ to find the corresponding
optimal partition.

The asymptotic expression of $\mathscr{W}$ and $\mathscr{W}_{\sslash}$ can be
derived by applying the Stirling approximation $n!\simeq\sqrt{2\pi
n}\left(\frac{n}{e}\right)^{n}$ to the binomial coefficients appearing into
the definition of the versions of the surprise described above, leading to
expressions like

$\displaystyle\binom{V+W-1}{W}$ $\displaystyle=$
$\displaystyle\frac{(V+W-1)!}{W!(V-1)!}\simeq\frac{\sqrt{2\pi(V+W-1)}(V+W-1)^{V+W-1}e^{-(V+W-1)}}{\sqrt{2\pi
W}W^{W}e^{-W}\sqrt{2\pi(V-1)}(V-1)^{V-1}e^{-(V-1)}}$ (67) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi(V-1)\left(\frac{W}{V+W-1}\right)}}\cdot\frac{(V+W-1)^{V+W-1}}{W^{W}(V-1)^{V-1}}$
$\displaystyle\simeq$ $\displaystyle\frac{1}{\sqrt{2\pi
Vq}}\cdot\frac{1}{\left(\frac{W}{V+W-1}\right)^{W}\left(\frac{V-1}{V+W-1}\right)^{V-1}}$
$\displaystyle\simeq$ $\displaystyle\frac{1}{\sqrt{2\pi
Vq}}\cdot\frac{1}{q^{W}(1-q)^{V-1}}=\frac{1}{\sqrt{2\pi\mu}\cdot\text{Geo}(V,W,q)}$

where $\text{Geo}(x,y,z)=z^{y}(1-z)^{x}$ is a geometric probability mass
function, $z=\frac{y}{x+y}$ and $q\simeq\frac{W}{V+W}$.

The time, in seconds, to find the corresponding optimal partitions is reported
in tables 3 and 4: overall, denser networks require more time.

##  Appendix C Enhanced detection of mesoscale structures

The enhanced hypergeometric distribution is a negative hypergeometric
distribution ‘conditional’ to the existence of connections within communities,
i.e. $l_{\bullet}>0$. It reads

$\text{EH}(l_{\bullet},w_{\bullet})=\left\\{\begin{array}[]{lr}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}\binom{w_{\circ}-1}{w_{\circ}-l_{\circ}}}{\binom{W-1}{L-1}},&0<l_{\bullet}<L\\\
\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-L}}{\binom{W}{L}},&l_{\bullet}=L\end{array}\right.$
(68)

the first expression being valid when $l_{\bullet}<L$ and the second
expression being valid in the ‘extreme’ case $l_{\bullet}=L$ (i.e. whenever
the total weight has to be redistributed on top of links which are all put
inside the communities). Normalization is ensured by the fact that

$\displaystyle+\sum_{l_{\bullet}=1}^{L-1}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}\cdot\sum_{w_{\bullet}-l_{\bullet}=0}^{W-L}\frac{\binom{(w_{\bullet}-l_{\bullet})+(l_{\bullet}-1)}{w_{\bullet}-l_{\bullet}}\binom{(W-2)-(l_{\bullet}-1)-(w_{\bullet}-l_{\bullet})}{(W-L)-(w_{\bullet}-l_{\bullet})}}{\binom{W-1}{L-1}}$
$\displaystyle+\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}\cdot\sum_{w_{\bullet}=L}^{W}\frac{\binom{w_{\bullet}-1}{w_{\bullet}-L}}{\binom{W}{L}}=\sum_{l_{\bullet}=1}^{L}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}=1-\frac{\binom{V_{\circ}}{L}}{\binom{V}{L}}$
(69)

where we have used the relationships
$\sum_{k=0}^{r}\binom{k+m}{k}\binom{n-m-k}{r-k}=\binom{n+1}{r}$,
$\sum_{k=r}^{n}\binom{k}{r}=\sum_{k=r}^{n}\binom{k}{k-r}=\binom{n+1}{r+1}$ and
$\sum_{k=0}^{r}\binom{m}{k}\binom{n-m}{r-k}=\binom{n}{r}$ and assumed that
$L<V_{\bullet}$. For what concerns the expected value, we have that

$\displaystyle\langle l_{\bullet}\rangle$ $\displaystyle=$
$\displaystyle\sum_{l_{\bullet}=1}^{L-1}l_{\bullet}\cdot\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}\cdot\sum_{w_{\bullet}-l_{\bullet}=0}^{W-L}\frac{\binom{(w_{\bullet}-l_{\bullet})+(l_{\bullet}-1)}{w_{\bullet}-l_{\bullet}}\binom{(W-2)-(l_{\bullet}-1)-(w_{\bullet}-l_{\bullet})}{(W-L)-(w_{\bullet}-l_{\bullet})}}{\binom{W-1}{L-1}}$
(70)
$\displaystyle+L\cdot\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}\cdot\sum_{w_{\bullet}=L}^{W}\frac{\binom{w_{\bullet}-1}{w_{\bullet}-L}}{\binom{W}{L}}=\sum_{l_{\bullet}=1}^{L}l_{\bullet}\cdot\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}=L\frac{V_{\bullet}}{V}$

as it should be: in fact, the expected value of (the number of) connections
should not depend on the weights put on the connections. Moreover,

$\displaystyle\langle w_{\bullet}-l_{\bullet}\rangle$ $\displaystyle=$
$\displaystyle\sum_{l_{\bullet}=1}^{L-1}\sum_{w_{\bullet}-l_{\bullet}=0}^{W-L}(w_{\bullet}-l_{\bullet})\cdot\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}\cdot\frac{\binom{(w_{\bullet}-l_{\bullet})+(l_{\bullet}-1)}{w_{\bullet}-l_{\bullet}}\binom{(W-2)-(l_{\bullet}-1)-(w_{\bullet}-l_{\bullet})}{(W-L)-(w_{\bullet}-l_{\bullet})}}{\binom{W-1}{L-1}}$
(71)
$\displaystyle+(W-L)\cdot\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}\cdot\sum_{w_{\bullet}=L}^{W}\frac{\binom{w_{\bullet}-1}{w_{\bullet}-L}}{\binom{W}{L}}$
$\displaystyle=$
$\displaystyle\sum_{l_{\bullet}=1}^{L-1}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}\cdot\sum_{w_{\bullet}-l_{\bullet}=0}^{W-L}(w_{\bullet}-l_{\bullet})\cdot\frac{\binom{(w_{\bullet}-l_{\bullet})+(l_{\bullet}-1)}{w_{\bullet}-l_{\bullet}}\binom{(W-2)-(l_{\bullet}-1)-(w_{\bullet}-l_{\bullet})}{(W-L)-(w_{\bullet}-l_{\bullet})}}{\binom{W-1}{L-1}}+(W-L)\cdot\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}$
$\displaystyle=$
$\displaystyle\sum_{l_{\bullet}=1}^{L-1}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V-V_{\bullet}}{L-l_{\bullet}}}{\binom{V}{L}}\cdot
l_{\bullet}\cdot\frac{W-L}{L}+(W-L)\cdot\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}$
$\displaystyle=$ $\displaystyle
L\left(\frac{V_{\bullet}}{V}-\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}\right)\left(\frac{W-L}{L}\right)+(W-L)\cdot\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}=(W-L)\frac{V_{\bullet}}{V}.$

The asymptotic expression of $\mathscr{E}$ can be derived by applying the
Stirling approximation $n!\simeq\sqrt{2\pi n}\left(\frac{n}{e}\right)^{n}$ to
the binomial coefficients appearing into the definition of the enhanced
surprise, leading to expressions like

$\displaystyle\binom{V}{L}\binom{W-1}{W-L}$ $\displaystyle=$
$\displaystyle\frac{V!}{L!(V-L)!}\cdot\frac{(W-1)!}{(W-L)!(L-1)!}$ (72)
$\displaystyle\propto$
$\displaystyle\frac{V^{V}}{L^{L}(V-L)^{V-L}}\cdot\frac{(W-1)^{W-1}}{(W-L)^{W-L}(L-1)^{L-1}}$
$\displaystyle\simeq$
$\displaystyle\frac{V^{V}}{L^{L}(V-L)^{V-L}}\cdot\frac{W^{W}}{(W-L)^{W-L}L^{L}}$
$\displaystyle=$
$\displaystyle\frac{1}{\left(\frac{L}{V-L}\right)^{L}\left(1-\frac{L}{V}\right)^{V}}\cdot\frac{1}{\left(1-\frac{L}{W}\right)^{W}\left(\frac{L}{W-L}\right)^{L}}$
$\displaystyle=$
$\displaystyle\frac{1}{\left(\frac{L}{V-L}\right)^{L}\left(1-\frac{L}{V}\right)^{V}}\cdot\frac{1}{\left(1-\frac{L}{W}\right)^{W-L}\left(\frac{L}{W}\right)^{L}}$
$\displaystyle=$
$\displaystyle\frac{1}{p^{L}(1-p)^{V-L}}\cdot\frac{1}{q^{W-L}(1-q)^{L}}$
$\displaystyle=$
$\displaystyle\text{BF}[V,W,L,p,q]=\text{Ber}(V,L,p)\cdot\text{Geo}(L,W-L,q)$

with a clear meaning of the symbols. A graphical representation of the
enhanced hypergeometric distribution is provided in fig. C.1.

Figure C.1: Examples of the enhanced hypergeometric distribution
$\text{EH}(l_{\bullet},w_{\bullet}|V,V_{\bullet},L,W)$. Top-left panel:
instances characterized by the values $V=45$, $V_{\bullet}=20$, $L=30$, $W=50$
and $l_{\bullet}=10,11,12,13,14,15$. Top-right panel: instances characterized
by the values $V=45$, $V_{\bullet}=20$, $L=10$, $W=20$ and
$l_{\bullet}=1,2,3,4,5,6,7$. Bottom panel: full view of the enhanced
hypergeometric distribution with parameters $V=80$, $V_{\bullet}=50$, $L=20$,
$W=50$ and $0<l_{\bullet}<L$.

Let us now briefly consider the distribution inducing the definition of the
enhanced bimodular surprise, i.e. the ‘multivariate’ enhanced hypergeometric
distribution. As before, it is ‘conditional’ to the existence of connections
within either modules one is looking at (i.e. either $l_{\bullet}>0$ or
$l_{\circ}>0$) and reads

$\text{MEH}(l_{\bullet},l_{\circ},w_{\bullet},w_{\circ})=\left\\{\begin{array}[]{lr}\frac{\binom{V_{\bullet}}{l_{\bullet}}\binom{V_{\circ}}{l_{\circ}}\binom{V_{\top}}{l_{\top}}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-l_{\bullet}}\binom{w_{\circ}-1}{w_{\circ}-l_{\circ}}\binom{w_{\top}-1}{w_{\top}-l_{\top}}}{\binom{W-1}{W-L}},&0<l_{\bullet},l_{\circ}<L\\\
\frac{\binom{V_{\bullet}}{L}}{\binom{V}{L}}\cdot\frac{\binom{w_{\bullet}-1}{w_{\bullet}-L}}{\binom{W}{L}},&l_{\bullet}=L,\>l_{\circ}=0\\\
\frac{\binom{V_{\circ}}{l_{\circ}}}{\binom{V}{L}}\cdot\frac{\binom{w_{\circ}-1}{w_{\circ}-L}}{\binom{W}{L}},&l_{\bullet}=0,\>l_{\circ}=L\\\
\end{array}\right.$ (73)

with a clear meaning of the symbols. Notice that the cases $l_{\bullet}=0$,
$0<l_{\circ}<L$ and $l_{\circ}=0$, $0<l_{\bullet}<L$ can be easily recovered
from the first row of the definition above.

##  Appendix D Pseudocodes for surprise optimization

Let us now provide a description of the algorithms implemented in the
‘SurpriseMeMore’ Python package (freely downloadable at the URL
https://github.com/EmilianoMarchese/SurpriseMeMore) for surprise optimization:
to this aim, we will also show the corresponding pseudocodes. Since an
exhaustive exploration of the space of all possible partitions of a network is
not feasible (especially when dealing with large graphs), one often proceeds
heuristically, in either an _agglomerative_ fashion or fixing the number of
clusters to be detected _a priori_.

Community detection on binary and weighted networks. The algorithm for
community detection on binary networks is called PACO (PArtitioning Cost
Optimization) and its pseudocode can be found in Nicolini2016 . The algorithm
implemented there is agglomerative in nature: in the initial configuration of
the algorithm, in fact, nodes are treated as singletons (hence, there are as
many communities as nodes). Then, edges are sorted according to the value of
their Jaccard index, in a decreasing order: as the Jaccard index of an edge is
proportional to the number of common neighbours of the two connected nodes, a
large value of it likely points out that the two nodes must be assigned to the
same module. For further details about PACO, see Nicolini2016 .

Here we adopt a similar strategy to the one implemented by PACO in the binary
case, with two main differences. The first one concerns the edge-sorting step:
edges, in fact, are no longer sorted according to the Jaccard index but
randomly ordered; this strategy avoids getting stuck in local minima more
effectively. The second one concerns the possibility of letting entire
communities merge (instead of moving single nodes from one to another) with a
certain probability $p_{mix}$.

1:function CalculateAndUpdateSurprise($C,C^{\prime}$)

2: $S\leftarrow calculateSurprise(C)$

3: $S^{\prime}\leftarrow calculateSurprise(C^{\prime})$

4: if $S^{\prime}<S$ then

5: $C\leftarrow C^{\prime}$

6: $S\leftarrow S^{\prime}$

7: end if

8:return $C$

9:end function

10:

11:$C\leftarrow$ array of length $N$, randomly initialized with $N$ different
integers, i.e. $1,2\dots N$;

12:$E\leftarrow$ randomly sorted edges;

13:for edge $(u,v)\in E$ do

14: if $C[u]\neq C[v]$ then

15: $C^{\prime}\leftarrow C$

16: if $0\leq uniformDistribution(0,1)\leq p_{mix}$ then

17: for node $t\in C[u]$ do

18: $C^{\prime}[t]\leftarrow C[v]$

19: end for

20: $C\leftarrow$CalculateAndUpdateSurprise($C,C^{\prime}$)

21: end if

22: if $p_{mix}\leq uniformDistribution(0,1)\leq p_{single}$ then

23: $C^{\prime}[u]\leftarrow C[v]$

24: $C\leftarrow$CalculateAndUpdateSurprise($C,C^{\prime}$)

25: end if

26: if $p_{single}\leq uniformDistribution(0,1)\leq 1$ then

27: $C^{\prime}[v]\leftarrow C[u]$

28: $C\leftarrow$CalculateAndUpdateSurprise($C,C^{\prime}$)

29: end if

30: end if

31: $\Rightarrow$ consider each node and swap its membership with the one of
each first neighbour of it; accept the move if $\mathscr{S}$, $\mathscr{W}$ or
$\mathscr{E}$ decreases;

32:end for

33:$\Rightarrow$ repeat the for-loop to improve the chance of finding the
optimal partition

Remarkably, it is also possible to define a version of the surprise-based
community-detection algorithm, where the number of clusters to be detected is
fixed a priori:

1:function CalculateAndUpdateSurprise($C,C^{\prime}$)

2: $S\leftarrow calculateSurprise(C)$

3: $S^{\prime}\leftarrow calculateSurprise(C^{\prime})$

4: if $S^{\prime}<S$ then

5: $C\leftarrow C^{\prime}$

6: $S\leftarrow S^{\prime}$

7: end if

8:return $C$

9:end function

10:

11:$C\leftarrow$ array of length $N$, randomly initialized with $n_{clust}$
different integers, i.e. $1,2\dots n_{clust}$;

12:$E\leftarrow$ randomly sorted edges;

13:for edge $(u,v)\in E$ do

14: $C^{\prime}\leftarrow C$

15: if $C^{\prime}[u]\neq C[v]$ then

16: $C^{\prime}[u]\leftarrow C[v]$

17: $C[v]\leftarrow C^{\prime}[u]$

18: $C\leftarrow$CalculateAndUpdateSurprise($C,C^{\prime}$)

19: else

20: $C^{\prime}[u]\leftarrow randomInteger(1,n_{clust})$

21: $C\leftarrow$CalculateAndUpdateSurprise($C,C^{\prime}$)

22: $C^{\prime\prime}[v]\leftarrow randomInteger(1,n_{clust})$

23: $C\leftarrow$CalculateAndUpdateSurprise($C,C^{\prime\prime}$)

24: end if

25: $\Rightarrow$ consider each node and swap its membership with the one of
each first neighbour of it; accept the move if $\mathscr{S}$, $\mathscr{W}$ or
$\mathscr{E}$ decreases;

26:end for

27:$\Rightarrow$ repeat the for-loop to improve the chance of finding the
optimal partition

‘Bimodular’ structures detection on binary and weighted networks. In this
case, the algorithm we implement fixes the number of clusters to be detected a
priori, since each node of the network can belong to only one _out of two_
possible subsets (e.g. the core and the periphery) - a circumstance that is
reflected into the initialization of the algorithm, i.e. a vector whose
entries are randomly chosen to be either 0 or 1. A second difference with
respect to the community detection case concerns the criterion according to
which edges are sorted (for purely bipartite graphs, in fact, the Jaccard