justified

# Entropy-based random models for hypergraphs

Fabio Saracco<EMAIL_ADDRESS>‘Enrico Fermi’ Research Center (CREF),
Via Panisperna 89A, 00184 Rome (Italy)  Institute for Applied Computing ‘Mauro
Picone’ (IAC), National Research Council, Via dei Taurini 19, 00185 Rome
(Italy)  IMT School for Advanced Studies, Piazza San Francesco 19, 55100 Lucca
(Italy)  Giovanni Petri  Network Science Institute, Northeastern University
London, Devon House 58 St. Katharine’s Way, London E1W 1LP (United Kingdom)
Renaud Lambiotte  Mathematical Institute, University of Oxford, Woodstock
Road, OX2 6GG Oxford (United Kingdom)  Tiziano Squartini  IMT School for
Advanced Studies, Piazza San Francesco 19, 55100 Lucca (Italy)  INdAM-GNAMPA
Istituto Nazionale di Alta Matematica ‘Francesco Severi’, P.le Aldo Moro 5,
00185 Rome (Italy)

###### Abstract

Network theory has primarily focused on pairwise relationships, disregarding
many-body interactions: neglecting them, however, can lead to misleading
representations of complex systems. Hypergraphs represent an increasingly
popular alternative for describing polyadic interactions: our innovation lies
in leveraging the representation of hypergraphs based on the incidence matrix
for extending the entropy-based framework to higher-order structures. In
analogy with the Exponential Random Graphs, we name the members of this novel
class of models Exponential Random Hypergraphs. Here, we focus on two explicit
examples, i.e. the generalisations of the Erdös-Rényi Model and of the
Configuration Model. After discussing their asymptotic properties, we employ
them to analyse real-world configurations: more specifically, _i)_ we extend
the definition of several network quantities to hypergraphs, _ii)_ compute
their expected value under each null model and _iii)_ compare it with the
empirical one, in order to detect deviations from random behaviours.
Differently from currently available techniques, ours is analytically
tractable, scalable and effective in singling out the structural patterns of
real-world hypergraphs differing significantly from those emerging as a
consequence of simpler, structural constraints.

###### pacs:

89.75.Fb; 02.50.Tt

##  I I. Introduction

Networks provide a powerful language to model interacting systems Newman
(2010); Caldarelli (2010). Within a network framework, the basic unit of
interaction, i.e. the _edge_ , involves two nodes and the complexity of the
structure as a whole arises from the combination of these units. Despite its
many successes, network science disregards certain aspects of interacting
systems, notably the possibility that more-than-two constituent units could
interact at a time Lambiotte  _et al._ (2019). Yet, it has been increasingly
shown that, for a variety of systems, interactions cannot be always decomposed
into a pairwise fashion and that neglecting higher-order ones can lead to an
incomplete, if not misleading, representation of them Lambiotte  _et al._
(2019); Battiston  _et al._ (2020); Bick  _et al._ (2023). Examples include
chemical reactions involving several compounds, coordination activities within
small teams of co-working people, brain activities mediated by groups of
neurons. Generally speaking, thus, modelling the joint coordination of
multiple entities calls for a framework that generalises the traditional,
edge-centered one. An increasingly popular alternative to support a science of
many-body interactions is provided by _hypergraphs_ as these versatile,
mathematical objects naturally allow nodes to interact in groups Battiston
_et al._ (2021).

Given the recent interest towards these objects, the definition of analytical
tools to study them is still in its infancy Chodrow (2020); Nakajima  _et al._
(2021); Musciotto  _et al._ (2021); Barthelemy (2022). This paper represents
our contribution to fill the gap: hereby, we extend the class of entropy-based
null models Jaynes (1957); Park and Newman (2004) to hypergraphs. These models
work by preserving a given set of quantities while randomising everything
else, hence destroying all, possible correlations between structural
properties except for those that are genuinely embodied by the constraints
themselves Squartini and Garlaschelli (2011); Cimini  _et al._ (2018);
Squartini  _et al._ (2018). The versatility of such an approach allows it to
be employed either in presence of full information (to quantify the level of
self-organisation of a given configuration by spotting out the patterns that
are incompatible with simpler, structural constraints Squartini  _et al._
(2013); Saracco  _et al._ (2017); Becatti  _et al._ (2019); Caldarelli  _et
al._ (2020); Neal  _et al._ (2021)) or in presence of partial information (to
infer the missing portion of a given configuration Parisi  _et al._ (2020)).

Our strategy for defining null models for hypergraphs is based on the
randomisation of their incidence matrix. The (generally, rectangular) table
contains information about the connectivity of nodes, i.e. the set of
hyperedges they belong to, and the connectivity of hyperedges, i.e. the set of
nodes they cluster. We will explicitly derive two members of the novel class
of models hereby named _Exponential Random Hypergraphs_ , i.e. the Random
Hypergraph Model (RHM, generalising the Erdös-Rényi Model) and the Hypergraph
Configuration Model (HCM, generalising the Configuration Model), and provide
an analytical characterisation of their behaviour. To this aim, we will
exploit the formal equivalence between the incidence matrix of a hypergraph
and the biadjacency matrix of a bipartite graph Nakajima  _et al._ (2021);
Barthelemy (2022). As a next step, we will employ the HCM to assess the
statistical significance of a number of patterns characterising several, real-
word hypergraphs.

##  II II. Formalism and basic quantities

A hypergraph can be defined as a pair $H(\mathcal{V},\mathcal{E}_{H})$ where
$\mathcal{V}$ is the set of vertices and $\mathcal{E}_{H}$ is the set of
hyperedges. Moving from the observation that the edge set $\mathcal{E}_{G}$ of
a traditional, binary, undirected graph $G(\mathcal{V},\mathcal{E}_{G})$ is a
subset of the power set of $\mathcal{V}$, several definitions of the hyperedge
set $\mathcal{E}_{H}$ have been provided: the two most popular ones are those
proposed in Berge (1967, 1973), where hyperedges ‘link’ one or more vertices,
and in Voloshin (2009), where hyperedges are allowed to be empty sets as well.
Hereby, we adopt the definition according to which $\mathcal{E}_{H}$ is a
multiset of the power set of $\mathcal{V}$: in words, we consider non-simple
hypergraphs, admitting loops and parallel edges (i.e. hyperedges involving
exactly the same set nodes) of any size, including 0 (corresponding to empty
hyperedges), and $|\mathcal{V}|$ (corresponding to hyperedges clustering all
vertices together).

$\mathbf{I}\>=\>\begin{tabular}[]{c|c|c|c|c|c|}&$\mathbf{e_{1}}$&$\mathbf{e_{2}}$&$\mathbf{e_{3}}$&$\mathbf{e_{4}}$&$\mathbf{e_{5}}$\\\
\hline\cr$\mathbf{n_{1}}$&0&1&1&0&1\\\ \hline\cr$\mathbf{n_{2}}$&0&0&1&1&1\\\
\hline\cr$\mathbf{n_{3}}$&0&1&1&1&0\\\ \hline\cr$\mathbf{n_{4}}$&1&0&0&1&1\\\
\hline\cr$\mathbf{n_{5}}$&1&1&1&1&0\\\ \hline\cr$\mathbf{n_{6}}$&0&1&0&0&0\\\
\hline\cr\end{tabular}$ Table 1: Example of incidence matrix describing a
binary, undirected hypergraph with $6$ nodes and $5$ hyperedges. The generic
entry $I_{i\alpha}$ is 1 if vertex $i$ belongs to hyperedge $\alpha$ and $0$
otherwise (see also Fig. 1).

As for traditional graphs, an algebraic representation of hypergraphs can be
devised as well. For a formal analogy with the traditional case, let us call
the cardinality of the set of nodes $|\mathcal{V}|=N$ and the cardinality of
the set of hyperedges $|\mathcal{E}_{H}|=L$: then, we consider the $N\times L$
table known as _incidence matrix_ , each row of which corresponds to a node
and each column of which corresponds to a hyperedge. If we indicate the
incidence matrix with $\mathbf{I}$, its generic entry $I_{i\alpha}$ will be
$1$ if vertex $i$ belongs to hyperedge $\alpha$ and $0$ otherwise. Notice that
the number of $1$s along each row can vary between $0$ and $L$, the former
case indicating an isolated node and the latter one indicating a node that
belongs to each hyperedge; similarly, the number of $1$s along each column can
vary between $0$ and $N$, the former case indicating an empty hyperedge and
the latter one indicating a hyperedge that includes all nodes. As explicitly
noticed elsewhere Nakajima  _et al._ (2021); Barthelemy (2022), representing a
hypergraph via its incidence matrix is equivalent at considering the bipartite
graph defined by the sets $\mathcal{V}$ and $\mathcal{E}_{H}$. Table 1
represents the incidence matrix describing the binary, undirected hypergraph
shown in Fig. 1.

Once the incidence matrix has been defined, several quantities needed for the
description of hypergraphs can be defined quite straightforwardly: for
example, the ‘degree of node $i$’ (hereby, _degree_) reads

$k_{i}=\sum_{\alpha=1}^{L}I_{i\alpha}$ (1)

and counts the number of hyperedges that are incident to it; analogously, the
‘degree of hyperedge $\alpha$’ (hereby, _hyperdegree_) reads

$h_{\alpha}=\sum_{i=1}^{N}I_{i\alpha}$ (2)

and counts the number of nodes it clusters. Both the sum of degrees and that
of hyperdegrees equal the total number of $1$s, i.e.
$\sum_{i=1}^{N}k_{i}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}I_{i\alpha}=\sum_{\alpha=1}^{L}\sum_{i=1}^{N}I_{i\alpha}=\sum_{\alpha=1}^{L}h_{\alpha}\equiv
T$. Importantly, a node degree no longer coincides with the number of its
neighbours: instead, it matches the number of hyperedges it belongs to; a
hyperdegree, instead, provides information on the hyperedge size.

Figure 1: Cartoon representation of the hypergraph whose incidence matrix is
shown in Table 1.

##  III III. Binary, undirected hypergraphs randomisation

An early attempt to define randomisation models for hypergraphs can be found
in Ghoshal  _et al._ (2009). Its authors, however, have just considered
hyperedges that are incident to triples of nodes - a framework that has been,
later, applied to the study of the World Trade Network Nakajima  _et al._
(2021). Considering the incidence matrix has two, clear advantages over the
tensor-based representation employed in Ghoshal  _et al._ (2009); Nakajima
_et al._ (2021): _i) generality_ , because the incidence matrix allows
hyperedges of any size to be handled at once; _ii) compactness_ , because the
order of the tensor $\mathbf{I}$ never exceeds two, hence allowing any
hypergraph to be represented as a traditional, bipartite graph.

In order to extend the rich set of null models induced by graph-specific,
global and local constraints to hypergraphs, we, first, need to identify the
quantities that can play this role within the novel setting. In what follows,
we will consider the total number of 1s, i.e. $T$, the degree and the
hyperdegree sequences, i.e. $\\{k_{i}\\}_{i=1}^{N}$ and
$\\{h_{\alpha}\\}_{\alpha=1}^{L}$ \- either separately or in a joint fashion;
moreover, we will distinguish between microcanonical and canonical
randomisation techniques.

###  III.1 A. Homogeneous benchmarks:
the Random Hypergraph Model (RHM)

The model is defined by just one, global constraint, that, in our case, reads
$T=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}I_{i\alpha}$. Its microcanonical version
extends the model by Erdös and Rényi Erdös and Rényi (1959) \- also known as
Random Graph Model (RGM) - to hypergraphs and prescribes to count the number
of incidence matrices that are compatible with a given, total number of 1s,
say $T^{*}$: they are

$\Omega_{\text{RHM}}=\binom{V}{T^{*}}$ (3)

with $V=NL$ being the total number of entries of the incidence matrix
$\mathbf{I}$. Once the total number of configurations composing the
microcanonical ensemble is determined, a procedure to generate them is needed:
in the case of the RHM, it simply boils down to reshuffling the entries of the
incidence matrix, a procedure ensuring that the total number of 1s is kept
fixed while any other correlation is destroyed.

The canonical version of the RHM, instead, extends the model by Gilbert
Gilbert (1959) and rests upon the constrained maximisation of Shannon entropy,
i.e.

$\mathscr{L}\equiv
S[P]-\sum_{i=0}^{M}\theta_{i}[P(\mathbf{I})C_{i}(\mathbf{I})-\langle
C_{i}\rangle]$ (4)

where $S[P]=-\sum_{\mathbf{I}\in\mathcal{I}}P(\mathbf{I})\ln P(\mathbf{I})$,
$C_{0}\equiv\langle C_{0}\rangle\equiv 1$ sums up the normalisation condition
and the remaining $M-1$ constraints represent proper, topological properties.
The sum defining Shannon entropy runs over the set $\mathcal{I}$ of incidence
matrices described in the introductory paragraph and known as canonical
ensemble.

Such an optimisation procedure defines the Exponential Random Hypergraph (ERH)
framework, described by the expression

$P(\mathbf{I})=\frac{e^{-H(\mathbf{I})}}{Z}=\frac{e^{-H(\mathbf{I})}}{\sum_{\mathbf{I}\in\mathcal{I}}e^{-H(\mathbf{I})}}=\frac{e^{-\sum_{i=1}^{M}\theta_{i}C_{i}(\mathbf{I})}}{\sum_{\mathbf{I}\in\mathcal{I}}e^{-\sum_{i=1}^{M}\theta_{i}C_{i}(\mathbf{I})}}.$
(5)

In the simplest case, the only, global constraint is represented by $T$ and
leads to the expression

$\displaystyle P(\mathbf{I})$ $\displaystyle=\frac{e^{-\theta
T(\mathbf{I})}}{\sum_{\mathbf{I}\in\mathcal{I}}e^{-\theta T(\mathbf{I})}}$
$\displaystyle=\frac{e^{-\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\theta
I_{i\alpha}}}{\sum_{\mathbf{I}\in\mathcal{I}}e^{-\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\theta
I_{i\alpha}}}$
$\displaystyle=\prod_{i=1}^{N}\prod_{\alpha=1}^{L}x^{I_{i\alpha}}\prod_{i=1}^{N}\prod_{\alpha=1}^{L}(1+x)^{-1}$
(6)

that can be re-written as

$\displaystyle
P(\mathbf{I})=\prod_{i=1}^{N}\prod_{\alpha=1}^{L}p^{I_{i\alpha}}(1-p)^{1-I_{i\alpha}}=p^{T(\mathbf{I})}(1-p)^{NL-T(\mathbf{I})}$
(7)

with $e^{-\theta}\equiv x$ and $p\equiv x/(1+x)$. The canonical ensemble, now,
includes all $N\times L$ rectangular matrices whose number of entries equating
1 ranges from 0 to $NL$. According to such a model, the entries of the
incidence matrix are i.i.d. Bernoulli random variables, i.e.
$I_{i\alpha}\sim\text{Ber}(p)$, $\forall\>i,\alpha$: as a consequence, the
total number of 1s, the degrees and the hyperdegrees obey Binomial
distributions - being all defined as sums of i.i.d. Bernoulli random
variables. Specifically, $T\sim\text{Bin}(NL,p)$, $k_{i}\sim\text{Bin}(L,p)$,
$\forall\>i$ and $h_{\alpha}\sim\text{Bin}(N,p)$, $\forall\>\alpha$, in turn,
implying that $\langle T\rangle_{\text{RHM}}=NLp$, $\langle
k_{i}\rangle_{\text{RHM}}=Lp$, $\forall\>i$ and $\langle
h_{\alpha}\rangle_{\text{RHM}}=Np$, $\forall\>\alpha$.

In order to ensure that $\langle T\rangle_{\text{RHM}}=T^{*}$, parameters have
to be tuned opportunely. To this aim, the likelihood maximisation principle
can be invoked Garlaschelli and Loffredo (2008): it prescribes to maximise the
function $\mathcal{L}(\theta)\equiv\ln P(\mathbf{I}^{*}|\theta)$ with respect
to the unknown parameter(s) that define it. Such a recipe leads us to find
$p=T^{*}/NL$, with $T^{*}=T(\mathbf{I}^{*})$ indicating the empirical value of
the constraint defining the RHM.

The RHM (also considered in Barthelemy (2022), although a derivation from
first principles is missing, there) is formally equivalent to the Bipartite
Random Graph Model (BiRGM) Saracco  _et al._ (2015). Such an identification is
guaranteed by our focus on non-simple hypergraphs.

####  III.1.1 1\. Estimation of the number of empty hyperedges

As non-simple hypergraphs admit the presence of empty as well as parallel
hyperedges, let us evaluate how frequently they appear in the ensembles
induced by our benchmarks.

To this aim, let us start by considering the probability for the generic
hyperedge $\alpha$ to be empty or, in other terms, that its hyperdegree
$h_{\alpha}$ is zero. Upon remembering that $h_{\alpha}\sim\text{Bin}(N,p)$,
one finds

$p_{0}\equiv(1-p)^{N}.$ (8)

Such a quantity leads us to find the probability of observing at least one,
empty hyperedge (or, in other terms, the complementary of the probability that
no hyperedge is empty), i.e.

$p_{1}\equiv 1-(1-p_{0})^{L}=1-[1-(1-p)^{N}]^{L},$ (9)

and the expected number of empty hyperedges: upon remembering that they are
i.i.d. Binomial random variables, one finds

$\langle N_{0}\rangle\equiv Lp_{0}=L(1-p)^{N}.$ (10)

Let us, now, inspect the behaviour of the aforementioned quantities on the
ensemble induced by the RHM, as the density of 1s in the incidence matrix,
i.e. $p=T/NL$, varies. The two regimes of interest are the _dense_ one,
defined by $T\rightarrow NL$, and the _sparse_ one, defined by $T\rightarrow
0$. In the dense case, one finds

$\displaystyle\lim_{T\to NL}p_{0}$ $\displaystyle=\lim_{T\to
NL}\left(1-\frac{T}{NL}\right)^{N}=0,$ (11)

a relationship inducing $p_{1}\xrightarrow[]{T\to NL}0$ and $\langle
N_{0}\rangle\xrightarrow[]{T\to NL}0$: in words, the probability of observing
empty hyperedges progressively vanishes as the density of 1s increases.
Consistently, in the sparse case one finds

$\displaystyle\lim_{T\to 0}p_{0}=\lim_{T\to
0}\left(1-\frac{T}{NL}\right)^{N}=1,$ (12)

a relationship inducing $p_{1}\xrightarrow[]{T\to 0}1$ and $\langle
N_{0}\rangle\xrightarrow[]{T\to 0}L$: in words, the probability of observing
empty hyperedges progressively rises as the density of 1s decreases.

To evaluate the density of 1s in the incidence matrix in correspondence of
which the transition from the sparse to the dense regime happens, let us
consider the case $N\gg 1$: more formally, this amounts to consider the
asymptotic framework defined by letting $N\to+\infty$ while posing $T=O(L)$,
i.e. $p=O(1/N)$. Since $p=T/NL$ and $h\equiv T/L$ remains finite, the
probability for a generic hyperedge to be empty remains finite as well:
consistently, the expected number of empty hyperedges becomes

$\lim_{N\to+\infty}\langle
N_{0}\rangle=\lim_{N\to+\infty}L\left(1-\frac{h}{N}\right)^{N}=Le^{-h};$ (13)

upon imposing $Le^{-h}\leq 1$, i.e. that the expected number of empty
hyperedges is _at most_ 1, one derives what may be called _filling threshold_
, corresponding to $h_{f}^{\text{RHM}}=\ln L$. In words, a value
$p>p_{f}^{\text{RHM}}\equiv h_{f}^{\text{RHM}}/N=\ln L/N$ ensures that the
expected number of empty hyperedges in our random hypergraph is strictly less
than one. The calculation above complements the result derived in Barthelemy
(2022) about the percolation threshold for hypergraphs (see also Appendix A).

####  III.1.2 2\. Estimation of the number of parallel hyperedges

Let us, now, move to considering the issue of parallel hyperedges. By
definition, two, parallel hyperedges $\alpha$ and $\beta$ are characterised by
identical columns: hence, their Hamming distance, defined as the number of
positions at which the corresponding symbols are different, is zero. More
formally,

$d_{\alpha\beta}\equiv\sum_{i=1}^{N}[I_{i\alpha}(1-I_{i\beta})+I_{i\beta}(1-I_{i\alpha})],$
(14)

a sum whose generic addendum is 1 in just two cases: either $I_{i\alpha}=1$
and $I_{i\beta}=0$ or $I_{i\alpha}=0$ and $I_{i\beta}=1$. Since
$d_{\alpha\beta}\sim\text{Bin}(N,2p(1-p))$, one finds that
$P(d_{\alpha\beta}=0)=[1-2p(1-p)]^{N}$ and $\langle
d_{\alpha\beta}\rangle=2p(1-p)N$, $\forall\>\alpha\neq\beta$. Notice that

$\displaystyle\lim_{T\to NL}P(d_{\alpha\beta}=0)$ $\displaystyle=\lim_{T\to
0}P(d_{\alpha\beta}=0)=1$ (15)

and

$\displaystyle\lim_{T\to NL}\langle d_{\alpha\beta}\rangle$
$\displaystyle=\lim_{T\to 0}\langle d_{\alpha\beta}\rangle=0$ (16)

as $(2T/NL)(1-T/NL)\to 0$ in both regimes: in words, the probability of
observing parallel hyperedges progressively rises both as a consequence of
having many 1s and as a consequence of having few 1s. Analogously for the
expected Hamming distance.

Let us, now, evaluate the expected Hamming distance between any, two
hyperedges $\alpha$ and $\beta$, within the asymptotic framework defined by
letting $N\to+\infty$ while posing $T=O(L)$, i.e. $p=O(1/N)$. Since $p=T/NL$
and $h\equiv T/L$ remains finite, one finds that

$\lim_{N\to+\infty}\langle
d_{\alpha\beta}\rangle=\lim_{N\to+\infty}\frac{2h}{N}\left(1-\frac{h}{N}\right)N=2h;$
(17)

upon imposing $2h\geq 1$, i.e. that the expected Hamming distance between any,
two hyperedges $\alpha$ and $\beta$ is _at least_ 1, one derives what may be
called _resolution threshold_ , corresponding to $h_{r}^{\text{RHM}}=1/2$. In
words, a value $p>p_{r}^{\text{RHM}}\equiv h_{r}^{\text{RHM}}/N=1/2N$ ensures
that, on average, any, two hyperedges $\alpha$ and $\beta$ differ by at least
one element (see also Appendix A).

####  III.1.3 3\. Estimation of the overlap between hyperedges

Finally, let us move to considering the issue of overlapping hyperedges. To
this aim, let us require that the probability for any, two hyperedges to
overlap amounts to 1. In formulas,

$1-(1-p^{2})^{N}=1,$ (18)

a relationship leading to the equation

$(1-p^{2})^{N}\simeq 1-Np^{2}=0$ (19)

that holds true in the sparse regime and further implies

$h^{2}=N;$ (20)

such a value, firstly derived in Barthelemy (2022), individuates the so-called
_percolation threshold_ , corresponding to $h_{p}^{\text{RHM}}=\sqrt{N}$. In
words, a value $p>p_{p}^{\text{RHM}}\equiv h_{p}^{\text{RHM}}/N=1/\sqrt{N}$
ensures that any, two hyperedges in our random hypergraph share at least one
node.

Since $h_{r}^{\text{RHM}}\leq h_{f}^{\text{RHM}}\leq h_{p}^{\text{RHM}}$, such
a result contributes to describe the following situation: by progressively
rising the parameter $p$, three, different thresholds are met, i.e. _i)_ the
resolution threshold $p_{r}^{\text{RHM}}=1/2N$, ensuring that, on average,
any, two hyperedges differ by at least one element; _ii)_ the filling
threshold $p_{f}^{\text{RHM}}=\ln L/N$, ensuring that the expected number of
empty hyperedges is strictly less than one; _iii)_ the percolation threshold
$p_{p}^{\text{RHM}}=1/\sqrt{N}$, ensuring that any, two hyperedges overlap.

###  III.2 B. Heterogeneous benchmarks:
the Hypergraph Configuration Model (HCM)

The number of constraints can be enlarged to include the degrees, i.e. the
sequence $\\{k_{i}\\}_{i=1}^{N}$, and the hyperdegrees, i.e. the sequence
$\\{h_{\alpha}\\}_{\alpha=1}^{L}$. Counting the number of configurations on
which both sequences match their empirical values is a hard task - although
numerical recipes to shuffle the entries of a rectangular matrix, while
preserving its marginals, exist Strona  _et al._ (2014); Carstens (2015);
Chodrow (2020). Solving the corresponding problem in the canonical framework
is, instead, straightforward. Indeed, Shannon entropy maximisation leads to

$\displaystyle P(\mathbf{I})$
$\displaystyle=\frac{e^{-\sum_{i=1}^{N}\alpha_{i}k_{i}(\mathbf{I})-\sum_{\alpha=1}^{L}\beta_{\alpha}h_{\alpha}(\mathbf{I})}}{\sum_{\mathbf{I}\in\mathcal{I}}e^{-\sum_{i=1}^{N}\alpha_{i}k_{i}(\mathbf{I})-\sum_{\alpha=1}^{L}\beta_{\alpha}h_{\alpha}(\mathbf{I})}}$
$\displaystyle=\frac{e^{-\sum_{i=1}^{N}\sum_{\alpha=1}^{L}(\alpha_{i}+\beta_{\alpha})I_{i\alpha}}}{\sum_{\mathbf{I}\in\mathcal{I}}e^{-\sum_{i=1}^{N}\sum_{\alpha=1}^{L}(\alpha_{i}+\beta_{\alpha})I_{i\alpha}}}$
$\displaystyle=\prod_{i=1}^{N}x_{i}^{k_{i}(\mathbf{I})}\prod_{\alpha=1}^{L}y_{\alpha}^{h_{\alpha}(\mathbf{I})}\prod_{i=1}^{N}\prod_{\alpha=1}^{L}(1+x_{i}y_{\alpha})^{-1},$
(21)

an expression that can be re-written as

$\displaystyle
P(\mathbf{I})=\prod_{i=1}^{N}\prod_{\alpha=1}^{L}p_{i\alpha}^{I_{i\alpha}}(1-p_{i\alpha})^{1-I_{i\alpha}}$
(22)

with $e^{-\alpha_{i}}\equiv x_{i}$, $\forall\>i$, $e^{-\beta_{\alpha}}\equiv
y_{\alpha}$, $\forall\>\alpha$ and $p_{i\alpha}\equiv
x_{i}y_{\alpha}/(1+x_{i}y_{\alpha})$, $\forall\>i,\alpha$. According to such a
model, the entries of the incidence matrix of a hypergraph are independent
random variables that obey different Bernoulli distributions, i.e.
$I_{i\alpha}\sim\text{Ber}(p_{i\alpha})$, $\forall\>i,\alpha$. As a
consequence, both degrees and hyperdegrees obey Poisson-Binomial
distributions, i.e.
$k_{i}\sim\text{PoissBin}(L,\\{p_{i\alpha}\\}_{\alpha=1}^{L})$, $\forall\>i$
and $h_{\alpha}\sim\text{PoissBin}(N,\\{p_{i\alpha}\\}_{i=1}^{N})$,
$\forall\>\alpha$ Saracco  _et al._ (2017).

In this case, solving the likelihood maximisation problem amounts to solve the
system of coupled equations

$\displaystyle k_{i}^{*}$
$\displaystyle=\sum_{\alpha=1}^{L}\frac{x_{i}y_{\alpha}}{1+x_{i}y_{\alpha}}=\sum_{\alpha=1}^{L}p_{i\alpha}=\langle
k_{i}\rangle,\>\forall\>i$ (23) $\displaystyle h_{\alpha}^{*}$
$\displaystyle=\sum_{i=1}^{N}\frac{x_{i}y_{\alpha}}{1+x_{i}y_{\alpha}}=\sum_{i=1}^{N}p_{i\alpha}=\langle
h_{\alpha}\rangle,\>\forall\>\alpha$ (24)

ensuring that $\langle k_{i}\rangle=k_{i}^{*}$, $\forall\>i$, $\langle
h_{\alpha}\rangle=h_{\alpha}^{*}$, $\forall\>\alpha$ (and, as a consequence,
$\langle T\rangle=T^{*}$). In case hypergraphs are sparse and in absence of
hubs

$p_{i\alpha}\simeq
x_{i}y_{\alpha}=\frac{k_{i}^{*}h_{\alpha}^{*}}{T^{*}},\>\forall\>i,\alpha.$
(25)

The HCM reduces to a ‘partial’ Configuration Model Saracco  _et al._ (2017)
when either the degree or the hyperdegree sequence is left unconstrained (see
also Appendix B). The canonical ensemble of each randomisation model (Table 2
in Appendix B sums up the set of constraints defining them) can be explicitly
sampled by considering each entry of $\mathbf{I}$, drawing a real number
$u_{i\alpha}\in U[0,1]$ and posing $I_{i\alpha}=1$ if $u_{i\alpha}\leq
p_{i\alpha}$, $\forall\>i,\alpha$.

The HCM is formally equivalent to the Bipartite Configuration Model (BiCM)
Saracco  _et al._ (2015). Such an identification is guaranteed by our focus on
non-simple hypergraphs.

####  III.2.1 1\. Estimation of the number of empty hyperedges

Let us, now, consider the probability for the generic hyperedge $\alpha$ to be
empty or, in other terms, that its hyperdegree $h_{\alpha}$ is zero. Upon
remembering that
$h_{\alpha}\sim\text{PoissBin}(N,\\{p_{i\alpha}\\}_{i=1}^{N})$, one finds

$p_{0}^{\alpha}\equiv\prod_{i=1}^{N}(1-p_{i\alpha}).$ (26)

Such a quantity leads to the probability of observing at least one, empty
hyperedge, now reading

$p_{1}\equiv
1-\prod_{\alpha=1}^{L}\left[1-p_{0}^{\alpha}\right]=1-\prod_{\alpha=1}^{L}\left[1-\prod_{i=1}^{N}(1-p_{i\alpha})\right],$
(27)

and to the expected number of empty hyperedges

$\langle
N_{0}\rangle\equiv\sum_{\alpha=1}^{L}p_{0}^{\alpha}=\sum_{\alpha=1}^{L}\prod_{i=1}^{N}(1-p_{i\alpha}).$
(28)

As in the previous section, let us, now, inspect the behaviour of the
aforementioned quantities on the ensemble induced by the HCM, as the density
of 1s in the incidence matrix varies. Although it depends on (the
heterogeneity of) the sets of coefficients $\\{x_{i}\\}_{i=1}^{N}$ and
$\\{y_{\alpha}\\}_{\alpha=1}^{L}$, general conclusions can still be drawn
within a simpler framework. To this aim, let us consider the functional form
reading

$p_{i\alpha}=\frac{zf_{i}g_{\alpha}}{1+zf_{i}g_{\alpha}},\>\forall\>i,\alpha$
(29)

where the vector of fitnesses $\\{f_{i}\\}_{i=1}^{N}$ accounts for the
heterogeneity of nodes, the vector of fitnesses
$\\{g_{\alpha}\\}_{\alpha=1}^{L}$ accounts for the heterogeneity of hyperedges
and $z$ tunes the density of 1s in the incidence matrix111Partial
Configuration Models are recovered upon posing either $f_{i}=1$, $\forall\>i$
or $g_{\alpha}=1$, $\forall\>\alpha$.. Within such a framework, the fitnesses
of the nodes and the fitnesses of the hyperedges can be drawn from any
distribution. The dense and sparse regimes are, now, defined by the positions
$z\to+\infty$ and $z\to 0$, respectively. In the dense case, one finds

$\displaystyle\lim_{z\to+\infty}p_{0}^{\alpha}$
$\displaystyle=\lim_{z\to+\infty}\prod_{i=1}^{N}(1-p_{i\alpha})$
$\displaystyle=\lim_{z\to+\infty}\prod_{i=1}^{N}\left(1-\frac{zf_{i}g_{\alpha}}{1+zf_{i}g_{\alpha}}\right)=0,$
(30)

a relationship inducing $p_{1}\xrightarrow[]{z\to+\infty}0$ and $\langle
N_{0}\rangle\xrightarrow[]{z\to+\infty}0$: in words, the probability of
observing empty hyperedges progressively vanishes as the density of 1s
increases. Consistently, in the sparse case one finds

$\displaystyle\lim_{z\to 0}p_{0}^{\alpha}$ $\displaystyle=\lim_{z\to
0}\prod_{i=1}^{N}(1-p_{i\alpha})$ $\displaystyle=\lim_{z\to
0}\prod_{i=1}^{N}\left(1-\frac{zf_{i}g_{\alpha}}{1+zf_{i}g_{\alpha}}\right)=1,$
(31)

a relationship inducing $p_{1}\xrightarrow[]{z\to 0}1$ and $\langle
N_{0}\rangle\xrightarrow[]{z\to 0}L$. In words, the probability of observing
empty hyperedges progressively rises as the density of 1s decreases.

####  III.2.2 2\. Estimation of the number of parallel hyperedges

Let us, now, focus on the issue of parallel hyperedges. As in the case of the
RHM, we consider the Hamming distance between the columns representing the two
hyperedges $\alpha$ and $\beta$. Since, now,
$d_{\alpha\beta}\sim\text{PoissBin}(N,\\{q_{i}^{\alpha\beta}\\}_{i=1}^{N})$,
where $q_{i}^{\alpha\beta}\equiv
p_{i\alpha}(1-p_{i\beta})+p_{i\beta}(1-p_{i\alpha})$ with
$p_{i\alpha}=zf_{i}g_{\alpha}/(1+zf_{i}g_{\alpha})$ and
$p_{i\beta}=zf_{i}g_{\beta}/(1+zf_{i}g_{\beta})$, one finds that
$P(d_{\alpha\beta}=0)=\prod_{i=1}^{N}(1-q_{i}^{\alpha\beta})$ and $\langle
d_{\alpha\beta}\rangle=\sum_{i=1}^{N}q_{i}^{\alpha\beta}$,
$\forall\>\alpha\neq\beta$. Notice that

$\displaystyle\lim_{z\to+\infty}P(d_{\alpha\beta}=0)$
$\displaystyle=\lim_{z\to 0}P(d_{\alpha\beta}=0)=1$ (32)

and

$\displaystyle\lim_{z\to 0}\langle d_{\alpha\beta}\rangle$
$\displaystyle=\lim_{z\to 0}\langle d_{\alpha\beta}\rangle=0$ (33)

as $q_{i}^{\alpha\beta}\to 0$ in both regimes. As in the case of the RHM, the
probability of observing parallel hyperedges progressively rises both as a
consequence of having many 1s and as a consequence of having few 1s.
Analogously for the expected Hamming distance.

When considering heterogeneous benchmarks, closed-form expressions for the
asymptotic estimates can be provided only in specific cases (see also Appendix
B).

Figure 2: Left panel: trend of
$p_{1}\xrightarrow[]{N\to+\infty}1-(1-e^{-h})^{L}$, i.e. the probability of
observing at least one empty hyperedge as a function of $T/NL$; evaluating it
in correspondence of $h_{f}^{\text{RHM}}=\ln L\simeq 6.907$ (whose related,
critical density value $p_{f}^{\text{RHM}}=h_{f}^{\text{RHM}}/N=\ln L/N\simeq
0.023$ individuates a steep transition separating the sparse regime from the
dense one) returns the value $1-(1-1/L)^{L}=0.6323$. Right panel: trend of
$p_{1}$ as defined by Eq. (27). Notice that the trend induced by the HCM is
shifted on the right with respect to the trend induced by the RHM: this is
probably due to the presence of small fitnesses increasing the probability of
observing at least one, empty hyperedge, hence requiring a larger value of $z$
to let $p_{1}$ vanish. The dense (sparse) regime is recovered for large
(small) values of $T/NL$. Each dot represents an average taken over an
ensemble of $1000$ configurations (explicitly sampled from either the RHM or
the HCM) and is accompanied by the corresponding $95\%$ confidence interval.

##  IV IV. Results

###  IV.1 A. Hypergraphs in the dense and sparse regime

Let us start by verifying the correctness of the estimations of the number of
empty and parallel hyperedges provided by our benchmarks: to this aim, we have
considered the values $N=300$ and $L=1000$ (for other values, see Appendix A
and Appendix B).

####  IV.1.1 1\. The Random Hypergraph Model

Let us start focusing on our homogeneous benchmark. The left panel of Fig. 2
depicts the trend of
$p_{1}=1-\left[1-\left(1-\frac{T}{NL}\right)^{N}\right]^{L}$. The left panel
of Fig. 3 depicts the trend of $\frac{\langle
N_{0}\rangle}{L}=\left(1-\frac{T}{NL}\right)^{N}=p_{0}$ and
$P(d_{\alpha\beta}=0)=\left[1-\frac{2T}{NL}\left(1-\frac{T}{NL}\right)\right]^{N}$.
Each quantity has been plotted as a function of $T/NL\in[10^{-6},1]$. The
dense (sparse) regime is recovered for large (small) values of $T/NL$. Each
dot of Figs. 2 and 3 represents an average taken over an ensemble of $1000$
configurations explicitly sampled from the RHM and is accompanied by the
corresponding $95\%$ confidence interval. A first, general comment concerns
the agreement between the analytical trends (solid lines) and the numerical
estimations (dots), confirming the correctness of our formulas.

##### The filling threshold.

Although the value of the filling threshold has been determined by inspecting
the asymptotic behaviour of $\langle N_{0}\rangle$, the quantity showing the
neatest transition from the sparse to the dense regime is
$p_{1}\xrightarrow[]{N\to+\infty}1-(1-e^{-h})^{L}$; evaluating it in
correspondence of $h_{f}^{\text{RHM}}=\ln L\simeq 6.907$ returns the value
$1-(1-1/L)^{L}\simeq 0.6323$ \- interestingly enough, letting $L$ grow
asymptotically as well, the previous value tends to $1-e^{-1}\simeq 0.6321$.

Upon considering that $p_{0}\xrightarrow[]{N\to+\infty}e^{-h}$, evaluating it
in correspondence of $h_{f}^{\text{RHM}}=\ln L\simeq 6.907$ returns the value
$1/L=0.001$: in words, although the filling threshold ensures that each,
single hyperedge is empty with an overall small probability, the likelihood of
observing at least one, empty hyperedge is still large (i.e. $\simeq 2/3)$.
The steepness of the trend of $p_{1}$, however, suggests it to quickly vanish
as the density of 1s in the incidence matrix crosses the value
$p_{f}^{\text{RHM}}=h_{f}^{\text{RHM}}/N=\ln L/N\simeq 0.023$.

Figure 3: Left panel: trends of $p_{0}\xrightarrow[]{N\to+\infty}e^{-h}$
(aquamarine line) and $P(d_{\alpha\beta}=0)\xrightarrow[]{N\to+\infty}e^{-2h}$
(dark-blue line) as a function of $T/NL$; evaluating them in correspondence of
$h_{r}^{\text{RHM}}=1/2$ (the related, critical density value is indicated in
light-blu) returns the values $e^{-{h_{r}^{\text{RHM}}}}=1/\sqrt{e}\simeq
0.606$ and $e^{-{2h_{r}^{\text{RHM}}}}=e^{-1}\simeq 0.367$; evaluating them in
correspondence of $h_{f}^{\text{RHM}}=\ln L\simeq 6.907$ (the related,
critical density value is indicated in magenta) returns the values
$e^{-{h_{f}^{\text{RHM}}}}=1/L=0.001$ and
$e^{-{2h_{f}^{\text{RHM}}}}=1/L^{2}=10^{-6}$. The red line indicates the
density value corresponding to $h=1$, ensuring that each, expected hyperdegree
amounts to 1. Right panel: trends of $p_{0}$, as defined by Eq. (26), and
$P_{0}$, i.e. the percentage of pairs of parallel hyperedges. As already
noticed, the shift on the right of the trend induced by the HCM is probably
due to the presence of small fitnesses, decreasing the magnitude of the set of
probability coefficients. The dense (sparse) regime is recovered for large
(small) values of $T/NL$. Each dot represents an average taken over an
ensemble of $1000$ configurations (explicitly sampled from either the RHM or
the HCM) and is accompanied by the corresponding $95\%$ confidence interval.
The yellow line indicates the density value in correspondence of which the
minimum, expected hyperdegree amounts to 1.

##### The resolution threshold.

The value of the resolution threshold has been determined by inspecting the
asymptotic behaviour of $\langle d_{\alpha\beta}\rangle$: evaluating
$P(d_{\alpha\beta}=0)\xrightarrow[]{N\to+\infty}e^{-2h}$ in correspondence of
$h_{r}^{\text{RHM}}=1/2$ returns the value $e^{-1}\simeq 0.367$.

Upon considering that $p_{0}\xrightarrow[]{N\to+\infty}e^{-h}$, evaluating it
in correspondence of $h_{r}^{\text{RHM}}=1/2$ returns the value
$1/\sqrt{e}\simeq 0.606$: in words, the resolution threshold ensures the
likelihood of observing any, two parallel hyperedges to be almost half of the
probability of each, single hyperedge to be empty, hence suggesting that
getting rid of parallel hyperedges is ‘easier’ than getting rid of empty ones.

##### Crossing the thresholds.

As Fig. 3 confirms, the probability for any, two hyperedges to be parallel
decreases faster than the probability of each, single hyperedge to be empty:
in words, while filling a large portion of the incidence matrix is not
required to observe a limited amount of parallel hyperedges, when empty
hyperedges are no longer observed, parallel hyperedges are no longer observed,
too. More formally, $h_{r}^{\text{RHM}}=1/2<h_{f}^{\text{RHM}}=\ln L$: hence,
progressively rising the density of 1s in the incidence matrix leads us
crossing the resolution threshold _before_ crossing the filling threshold. In
more quantitative terms, evaluating $\langle
N_{0}\rangle\xrightarrow[]{N\to+\infty}Le^{-h}$ in correspondence of
$h_{r}^{\text{RHM}}=1/2$ returns the value $L/\sqrt{e}\simeq 606$ which is
larger than $\langle d_{\alpha\beta}\rangle\xrightarrow[]{N\to+\infty}2h$
evaluated in correspondence of $h_{r}^{\text{RHM}}=1/2$ \- returning the value
1 by definition.

####  IV.1.2 2\. The Hypergraph Configuration Model

Let us, now, focus on the trends induced by our heterogeneous benchmarks.
Following the procedure described in the previous sections, both the fitnesses
of nodes and those of hyperedges were drawn from a Pareto distribution222Other
fat-tailed distributions were considered: qualitatively, results do not
change. with $\alpha=2$. The right panel of Fig. 2 depicts the trend of
$p_{1}$, as defined by Eq. (27). The right panel of Fig. 3 depicts the trend
of $p_{0}$, as defined by Eq. (26), and the arithmetic mean $P_{0}\equiv
2\sum_{\alpha<\beta}P(d_{\alpha\beta}=0)/L(L-1)$ with
$P(d_{\alpha\beta}=0)=\prod_{i=1}^{N}(1-q_{i}^{\alpha\beta})$, i.e. the
average of the value $P(d_{\alpha\beta}=0)$ over all pairs of hyperedges - a
choice dictated by the will of accounting for the heterogeneity of different
pairs of hyperedges. Both quantities have been plotted as a function of
$\langle T\rangle/NL\in[10^{-6},1]$ where $\langle
T\rangle=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}zf_{i}g_{\alpha}/(1+zf_{i}g_{\alpha})$.
The dense (sparse) regime is recovered for large (small) values of $z$. Each
dot of Figs. 2 and 3 represents an average taken over an ensemble of $1000$
configurations explicitly sampled from the HCM and is accompanied by the
corresponding $95\%$ confidence interval. As in the case of the RHM, the
agreement between the analytical trends (solid lines) and the numerical
estimations (dots) confirms the correctness of our formulas.

Comparing the analytical estimation of (the trend of) the percentage of pairs
of parallel hyperedges with its sample counterpart reveals the former to
provide a quite accurate estimation of the latter. Besides, we explicitly
notice that the trends characterising the HCM are sistematically shifted on
the right: this is probably due to the presence of small fitnesses that, at
least in case of Fig. 2, increase the probability of observing at least one,
empty hyperedge, hence requiring a larger value of $z$ to let $p_{1}$ vanish.
Such a result seems to be confirmed by the trends depicted in Fig. 3.

Deriving an explicit expression for the filling threshold in the case of the
HCM is a rather difficult task; still, we can proceed in a purely numerical
fashion and individuate the value of the density of 1s in the incidence matrix
guaranteeing that the minimum, expected hyperdegree amounts to 1: as evident
from Fig. 3, such a threshold still lies in the right tail of the trends
induced by the HCM - although its, precise, numerical value depends on those
of the fitnesses.

###  IV.2 B. Solving the HCM on real-world hypergraphs

In order to test our benchmarks on real-world configurations, we have focused
on a number of data sets taken from Austin R. Benson’s
website333https://www.cs.cornell.edu/~arb/data/, i.e. the contact-primary-
school, the email-Enron and the NDC-classes ones.

To have an explicit expression for the expected value of any quantity of
interest, let us comment on the goodness of the approximation

$p_{i\alpha}\simeq
x_{i}y_{\alpha}=\frac{k_{i}^{*}h_{\alpha}^{*}}{T^{*}},\>\forall\>i,\alpha.$
(34)

As Fig. 7 in Appendix A shows, it is is quite accurate for each data set
considered here - in fact, one can safely assume that $x_{i}\simeq
k_{i}^{*}/\sqrt{T^{*}}$, $\forall\>i$ and $h_{\alpha}^{*}/\sqrt{T^{*}}$,
$\forall\>\alpha$.

###  IV.3 C. ‘Hypergraph to graph’ projection

The canonical formalism that we have adopted leads to factorisable
distributions, i.e. distributions that can be written as a product of pair-
wise probability distributions; this allows the expectation of several
quantities of interest to be evaluated analytically. To this aim, let us
consider the matrix introduced in Battiston  _et al._ (2020), reading

$\mathbf{W}=\mathbf{I}\cdot\mathbf{I}^{T}-\mathbf{K}$ (35)

with $\mathbf{K}$ being the diagonal matrix whose $i$-th entry reads $k_{i}$;
according to the definition above

$w_{ij}=\sum_{\alpha=1}^{L}I_{i\alpha}I_{j\alpha}-\delta_{ij}k_{i}$ (36)

counts the number of hyperedges nodes $i$ and $j$ are connected by - more
explicitly, $w_{ij}=\sum_{\alpha=1}^{L}I_{i\alpha}I_{j\alpha}$, $i\neq j$ and
$w_{ii}=\sum_{\alpha=1}^{L}I_{i\alpha}I_{i\alpha}-k_{i}=\sum_{\alpha=1}^{L}I_{i\alpha}-k_{i}=k_{i}-k_{i}=0$.
The null models discussed so far can be employed to calculate $\langle
w_{ij}\rangle$, $i\neq j$ that, in a perfectly general fashion, reads

$\langle w_{ij}\rangle=\sum_{\alpha=1}^{L}\langle
I_{i\alpha}I_{j\alpha}\rangle=\sum_{\alpha=1}^{L}\langle
I_{i\alpha}\rangle\langle
I_{j\alpha}\rangle=\sum_{\alpha=1}^{L}p_{i\alpha}p_{j\alpha};$ (37)

the total number of hyperedges shared by node $i$ with any other node in the
hypergraph (in a sense, its ‘strength’ - see also Fig. 4) can be computed as

$\sigma_{i}=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}w_{ij}=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\sum_{\alpha=1}^{L}I_{i\alpha}I_{j\alpha}$ (38)

whose expected value reads

$\displaystyle\langle\sigma_{i}\rangle$
$\displaystyle=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\langle
w_{ij}\rangle=\sum_{\alpha=1}^{L}p_{i\alpha}[\langle h_{\alpha}\rangle-
p_{i\alpha}].$ (39)

As further confirmed by Fig. 8 in Appendix A, the approximation provided by
Eq. (34) allows us to pose $\langle\sigma_{i}\rangle_{\text{HCM}}\simeq
k_{i}^{*}\sum_{\alpha=1}^{L}(h_{\alpha}^{*})^{2}/T^{*}$, $\forall\>i$.
Interestingly, as Fig. 5a shows, the HCM overestimates the extent to which any
two nodes of the email-Enron data set overlap: in words, such a real-world
hypergraph is more compartmentalised than expected.

Let us, now, extend the concept of _assortativity_ to hypergraphs. To this
aim, we consider the quantity named _average incident hyperedges degree_ ,
defined as

$k_{i}^{nn}=\sum_{\alpha=1}^{L}\frac{I_{i\alpha}h_{\alpha}}{k_{i}}=\frac{\sigma_{i}+k_{i}}{k_{i}}=\frac{\sigma_{i}}{k_{i}}+1\simeq\frac{\sigma_{i}}{k_{i}}$
(40)

and representing the arithmetic mean of the degrees of the hyperedges
including node $i$. An analytical approximation of its expected value can be
provided as well:

$\displaystyle\langle k_{i}^{nn}\rangle$
$\displaystyle\simeq\sum_{\alpha=1}^{L}\frac{p_{i\alpha}[\langle
h_{\alpha}\rangle+1-p_{i\alpha}]}{\langle k_{i}\rangle}$
$\displaystyle=\frac{\langle\sigma_{i}\rangle+\langle k_{i}\rangle}{\langle
k_{i}\rangle}=\frac{\langle\sigma_{i}\rangle}{\langle
k_{i}\rangle}+1\simeq\frac{\langle\sigma_{i}\rangle}{\langle k_{i}\rangle}.$
(41)

Figure 4: Graphical representation of the ‘hypergraph to graph’ projection.
Notice that while the degree of each node (i.e. the number of hyperedges that
are incident to it) reads $\\{k_{i}\\}=(3,3,3,3,4,1)$, the sigma of each node
(i.e. the total number of hyperedges shared by it with any other node or,
equivalently, the strength on the projection) reads
$\\{\sigma_{i}\\}=(8,8,9,6,10,3)$ and the kappa of each node (i.e. the number
of nodes it shares an hyperedge with or, equivalently, the degree on the
projection) reads $\\{\kappa_{i}\\}=(5,4,5,4,5,3)$.

###  IV.4 D. Disparity ratio and degree in the projection

More information about the patterns shaping real-world hypergraphs can be
obtained upon defining the ratio $f_{ij}=w_{ij}/\sigma_{i}$, $i\neq j$ that
induces the quantity

$Y_{i}=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}f_{ij}^{2}=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\frac{w_{ij}^{2}}{\sigma_{i}^{2}},$ (42)

known as _disparity ratio_ and quantifying the (un)evenness of the
distribution of the weights constituting the strength of node $i$ over the
$\kappa_{i}=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\Theta[w_{ij}]\equiv\sum_{\begin{subarray}{c}j=1\\\
(j\neq i)\end{subarray}}^{N}a_{ij}$ links characterising its connectivity -
since $a_{ij}=1$ if nodes $i$ and $j$ share, at least, one hyperedge,
$\kappa_{i}$ is the degree of node $i$ in the projection of the hypergraph
(see also Fig. 4 and Fig. 5b). Since, under the RHM,
$w_{ij}\sim\text{Bin}(L,p^{2})$, we find that $\langle
a_{ij}\rangle=1-(1-p^{2})^{L}$, i.e. the expected value of $a_{ij}$ coincides
with the probability of observing a non-zero overlap. Under the HCM, instead,
$w_{ij}\sim\text{PoissBin}(L,\\{p_{i\alpha}p_{j\alpha}\\}_{\alpha=1}^{L})$,
hence

$\langle a_{ij}\rangle=1-\prod_{\alpha=1}^{L}(1-p_{i\alpha}p_{j\alpha}).$ (43)

Let us also notice that

$Y_{i}=\frac{1}{\kappa_{i}}$ (44)

in case weights are equally distributed among the connections established by
node $i$, i.e. $w_{ij}=a_{ij}\sigma_{i}/\kappa_{i}$, $i\neq j$. Any larger
value signals an excess concentration of weight in one or more links. An
analytical approximation of the expected value of the disparity ratio of node
$i$ can be provided as well:

$\langle Y_{i}\rangle\simeq\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\frac{\langle
w_{ij}^{2}\rangle}{\langle\sigma_{i}^{2}\rangle}.$ (45)

Figure 5: Scatter plots between the empirical values $\\{\sigma_{i}\\}$,
$\\{\kappa_{i}\\}$, $\\{Y_{i}\\}$, $\\{\text{CEC}_{i}\\}$ and the expected
values $\\{\langle\sigma_{i}\rangle_{\text{HCM}}\\}$,
$\\{\langle\kappa_{i}\rangle_{\text{HCM}}\\}$,
$\langle\\{Y_{i}\rangle_{\text{HCM}}\\}$,
$\\{\langle\text{CEC}_{i}\rangle_{\text{HCM}}\\}$. The HCM overestimates the
extent to which any, two nodes overlap, as well as the CEC; the disparity
ratio, instead, is underestimated by it. These results can be understood by
considering that the HCM just constrains the degree sequences, hence inducing
an ensemble where connections are ‘distributed’ more evenly than observed.

Contrarily to what has been previously observed, the expected value of the
disparity ratio cannot be always safely decomposed as a ratio of expected
values - not even if the ‘full’ HCM is employed. In fact, while this
approximation works relatively well for the contact-primary-school data set,
it does not for the email-Enron and the NDC-classes ones (see also Fig. 10 in
Appendix C). For this reason, the expected value of the disparity ratio has
been evaluated by explicitly sampling the ensemble of incidence matrices
induced by the ‘full’ HCM. In any case, as Fig. 5c shows, such a null model
underestimates the disparity ratio characterising each node of the email-Enron
data set: in words, the empirical overlap between any two nodes is (much) less
evenly ‘distributed’ than expected. A similar conclusion can be drawn by
considering
$\langle\kappa_{i}\rangle=\sum_{i=1}^{N}[1-\prod_{\alpha=1}^{L}(1-p_{i\alpha}p_{j\alpha})]$:
as Fig. 5b shows, the degree of the nodes in the projection tend to be
significantly smaller than expected, meaning that hyperedges concentrate on
fewer edges than expected. This observation is in line with recent works
showing the encapsulation and ‘simpliciality’ of real-world hypergraphs LaRock
and Lambiotte (2023); Landry  _et al._ (2024).

###  IV.5 F. Eigenvector centrality

Centrality measures for hypergraphs have been defined as well. An example is
provided by the _clique motif eigenvector centrality_ (CEC), defined in Benson
(2019) (see also Appendix D): $\text{CEC}_{i}$ corresponds to the $i$-th entry
of the Perron-Frobenius eigenvector of $\mathbf{W}$. As Fig. 5d shows, the HCM
underestimates the CEC as well: such a result can be understood by considering
that the HCM constrains only the degree sequences, hence inducing an ensemble
where connections are ‘distributed’ more evenly than observed, an evidence
letting the nodes overlap more, thus causing the entries of
$\langle\mathbf{W}\rangle$ to be overall larger and less dissimilar, as well
as those of its Perron-Frobenius eigenvector.

Figure 6: Top panels: projections of the contact-primary-school, email-Enron
and NDC-classes data sets onto the layer of nodes. Bottom panels: validated
counterparts of the aforementioned projections: any two nodes are linked if
they share a significantly large number of hyperedges. Communities have been
detected by running the Louvain algorithm.

###  IV.6 G. Confusion matrix

Let us, now, consider the set of properties composing the so-called _confusion
matrix_ (see also Appendix E). They are known as the _true positive rate_
(TPR), i.e. the percentage of 1s correctly recovered by a given method and
whose expected value reads

$\langle\text{TPR}\rangle=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{I_{i\alpha}p_{i\alpha}}{T};$
(46)

the _specificity_ (SPC), i.e. the percentage of 0s correctly recovered by a
given method and whose expected value reads

$\langle\text{SPC}\rangle=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{(1-I_{i\alpha})(1-p_{i\alpha})}{NL-T};$
(47)

the _positive predictive value_ (PPV), i.e. the percentage of 1s correctly
recovered by a given method with respect to the total number of 1s predicted
by it and whose expected value reads

$\langle\text{PPV}\rangle=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{I_{i\alpha}p_{i\alpha}}{\langle
T\rangle};$ (48)

and the _accuracy_ (ACC), measuring the overall performance of a given method
in correctly placing both 1s and 0s and whose expected value reads

$\langle\text{ACC}\rangle=\frac{\langle\text{TP}\rangle+\langle\text{TN}\rangle}{NL}.$
(49)

Results on the confusion matrix of a number of real-world hypergraphs reveal
that the large sparsity of the latter ones makes it difficult to reproduce the
TPR and the PPV (see also Table 3 in Appendix E); on the other hand, the
capability of the HCM (both in its ‘full’ and approximated version) to
reproduce the density of 1s - and, as a consequence, the density of 0s -
ensures the SPC to be recovered quite precisely, in turn ensuring the overall
accuracy of the model to be large (for an overall evaluation of the
performance of the HCM in reproducing real-world hypergraphs, see also Table 4
in Appendix F).

###  IV.7 H. Community detection

Communities are commonly understood as densely connected groups of nodes.
Representing an hypergraph via its incidence matrix allows this statement to
be made more precise from a statistical perspective: in fact, the null models
discussed so far can be employed to test if any two nodes share a
significantly large number of hyperedges - hence can be clustered together,
should this be the case.

To this aim, we can adapt the recipe proposed in Saracco  _et al._ (2017) to
project bipartite networks and summarised in the following. One, first,
computes

$\text{p-value}(w_{ij}^{*})=\sum_{x\geq w_{ij}^{*}}f(x)$ (50)

for each pair of nodes; $f(x)$ depends on the chosen null model: in case the
RHM is employed, it coincides with the Binomial distribution
$\text{Bin}(x|L,p)$; in case the HCM is employed, it coincides with the
Poisson-Binomial distribution
$\text{PoissBin}(x|L,\\{p_{i\alpha}p_{j\alpha}\\}_{\alpha=1}^{L})$. Second,
one implements the FDR procedure, designed to handle multiple tests of
hypothesis Benjamini and Hochberg (1995): in practice, after ranking the
p-values in increasing order, i.e.
$\text{p-value}_{1}\leq\text{p-value}_{2}\leq\dots\leq\text{p-value}_{n}$, one
individuates the largest integer $\hat{i}$ satisfying the condition

$\text{p-value}_{\hat{i}}\leq\frac{\hat{i}t}{n}$ (51)

where $n=N(N-1)/2$ and $t$ is the single-test significance level, set to
$0.01$ in the present analysis. Third, one links the (pairs of) nodes whose
related p-value is smaller than the aforementioned threshold.

Figure 6 shows the result returned by the Louvain algorithm, run on the
validated projections: as noticed elsewhere Saracco  _et al._ (2017); Pratelli
_et al._ (2024); Guarino  _et al._ (2024), the detection of mesoscale
structures is enhanced if carried out on filtered topologies.

##  V V. Conclusions

Our paper contributes to current research on hypergraphs by extending the
constrained entropy-maximisation framework to incidence matrices, i.e. their
simplest, tabular representation. Differently from the currently-available
techniques Chodrow (2020), our methodology has the advantage of being
analytically tractable, scalable and versatile enough to be straightforwardly
extensible to directed and/or weighted hypergraphs.

Beside leading to results whose relevance is mostly theoretical (i.e. the
individuation of different regimes for higher-order structures and the
estimation of the actual impact of empty and parallel hyperedges on the
analysis of empirical systems), our models prove to be particularly useful
when employed as benchmarks for real-world systems, i.e. for detecting
patterns that are not imputable to purely random effects. Specifically, our
results suggest that real-world hypergraphs are characterised by a degree of
self-organisation that is absolutely non-trivial (see also Appendix F).

This is even more surprising when considering that our results are obtained
under a benchmark such as the HCM, i.e. a null model constraining both the
degree and the hyperdegree sequences: since it overestimates the extent to
which any two nodes overlap - a result whose relevance becomes evident as soon
as one considers the effects that higher-order structures have on spreading
and cooperation processes de Arruda  _et al._ (2021); Iacopini  _et al._
(2022); St-Onge  _et al._ (2022) \- our future efforts will be directed
towards the analysis of benchmarks constraining non-linear quantities such as
the co-occurrences between nodes and/or hyperedges.

##  VI VI. Acknowledgements

RL acknowledges support from the EPSRC grants EP/V013068/1, EP/V03474X/1 and
EP/Y028872/1. TS acknowledges support from SoBigData.it that receives funding
from European Union – NextGenerationEU – National Recovery and Resilience Plan
(Piano Nazionale di Ripresa e Resilienza, PNRR) – Project: ‘SoBigData.it –
Strengthening the Italian RI for Social Mining and Big Data Analytics’ - Prot.
IR0000013 – Avviso n. 3264 del 28/12/2021. FS acknowledges support from the
project ‘CODE - Coupling Opinion Dynamics with Epidemics’, funded under PNRR
Mission 4 ‘Education and Research’ - Component C2 - Investment 1.1 - Next
Generation EU ‘Fund for National Research Program and Projects of Significant
National Interest’ PRIN 2022 PNRR, grant code P2022AKRZ9.

## References

  * Newman (2010) M. Newman,  _Networks: an introduction_ (Oxford University Press, New York, USA, 2010).
  * Caldarelli (2010) G. Caldarelli,  _Scale-free networks_ (Oxford University Press, Oxford, UK, 2010).
  * Lambiotte  _et al._ (2019) R. Lambiotte, M. Rosvall, and I. Scholtes, Nature Physics 15, 313 (2019).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J. G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Bick  _et al._ (2023) C. Bick, E. Gross, H. A. Harrington, and M. T. Schaub, SIAM Review 65, 686 (2023).
  * Battiston  _et al._ (2021) F. Battiston, E. Amico, A. Barrat, G. Bianconi, G. F. de Arruda, B. Franceschiello, I. Iacopini, S. Kéfi, V. Latora, Y. Moreno, M. M. Murray, T. P. Peixoto, F. Vaccarino, and G. Petri, Nature Physics 17, 1093–1098 (2021).
  * Chodrow (2020) P. S. Chodrow, Journal of Complex Networks 8, 3 (2020).
  * Nakajima  _et al._ (2021) K. Nakajima, K. Shudo, and N. Masuda, arXiv preprint arXiv:2106.12162 (2021).
  * Musciotto  _et al._ (2021) F. Musciotto, F. Battiston, and R. N. Mantegna, Communications Physics 4, 218 (2021).
  * Barthelemy (2022) M. Barthelemy, Phys. Rev. E 106, 064310 (2022).
  * Jaynes (1957) E. Jaynes, The Physical Review 106, 181 (1957).
  * Park and Newman (2004) J. Park and M. E. J. Newman, Physical Review E 70, 66117 (2004).
  * Squartini and Garlaschelli (2011) T. Squartini and D. Garlaschelli, New Journal of Physics 13, 083001 (2011).
  * Cimini  _et al._ (2018) G. Cimini, T. Squartini, F. Saracco, D. Garlaschelli, A. Gabrielli, and G. Caldarelli, Nature Reviews Physics 1, 58 (2018).
  * Squartini  _et al._ (2018) T. Squartini, G. Caldarelli, G. Cimini, A. Gabrielli, and D. Garlaschelli, Physics Reports 757, 1 (2018).
  * Squartini  _et al._ (2013) T. Squartini, I. van Lelyveld, and D. Garlaschelli, Scientific Reports 3, 3357 (2013).
  * Saracco  _et al._ (2017) F. Saracco, M. J. Straka, R. Di Clemente, A. Gabrielli, G. Caldarelli, and T. Squartini, New Journal of Physics 19, 16 (2017).
  * Becatti  _et al._ (2019) C. Becatti, G. Caldarelli, R. Lambiotte, and F. Saracco, Palgrave Communications 5, 1 (2019).
  * Caldarelli  _et al._ (2020) G. Caldarelli, R. D. Nicola, F. D. Vigna, M. Petrocchi, F. Saracco, F. D. Vigna, M. Petrocchi, and F. Saracco, Communications Physics 3, 1 (2020).
  * Neal  _et al._ (2021) Z. P. Neal, R. Domagalski, and B. Sagan, Scientific Reports 11, 23929 (2021).
  * Parisi  _et al._ (2020) F. Parisi, T. Squartini, and D. Garlaschelli, New Journal of Physics 22, 053053 (2020).
  * Berge (1967) C. Berge,  _Graphes et hypergraphes_ (Dunod, Paris, FR, 1967).
  * Berge (1973) C. Berge,  _Graphes et hypergraphes_ (North-Holland, Amsterdam, NL, 1973).
  * Voloshin (2009) V. I. Voloshin,  _Introduction to graph and hypergraph theory_ (Nova Science Publishers, 2009).
  * Ghoshal  _et al._ (2009) G. Ghoshal, V. Zlatic, G. Caldarelli, and M. E. J. Newman, Physical Review E 79, 066118 (2009).
  * Erdös and Rényi (1959) P. Erdös and A. Rényi, Publicationes Mathematicae 6, 290 (1959).
  * Gilbert (1959) E. N. Gilbert, The Annals of Mathematical Statistics 30, 1141 (1959).
  * Garlaschelli and Loffredo (2008) D. Garlaschelli and M. I. Loffredo, Physical Review E 78, 1 (2008).
  * Saracco  _et al._ (2015) F. Saracco, R. D. Clemente, A. Gabrielli, and T. Squartini, Scientific Reports 5, 10595 (2015).
  * Strona  _et al._ (2014) G. Strona, D. Nappo, F. Boccacci, S. Fattorini, and J. San-Miguel-Ayanz, Nature Communications 5, 4114 (2014).
  * Carstens (2015) C. J. Carstens, Physical Review E 91, 1 (2015).
  * LaRock and Lambiotte (2023) T. LaRock and R. Lambiotte, Journal of Physics: Complexity 4, 045007 (2023).
  * Landry  _et al._ (2024) N. W. Landry, J.-G. Young, and N. Eikmeier, EPJ Data Science 13, 17 (2024).
  * Benson (2019) A. R. Benson, SIAM Journal on Mathematics of Data Science 1, 293 (2019).
  * Benjamini and Hochberg (1995) Y. Benjamini and Y. Hochberg, Journal of the Royal Statistical Society. Series B (Methodological) 57, 289 (1995).
  * Pratelli  _et al._ (2024) M. Pratelli, F. Saracco, and M. Petrocchi, PNAS Nexus 3, pgae177 (2024), https://academic.oup.com/pnasnexus/article-pdf/3/5/pgae177/58004712/pgae177.pdf .
  * Guarino  _et al._ (2024) S. Guarino, A. Mounim, G. Caldarelli, and F. Saracco, “Verified authors shape x/twitter discursive communities,” (2024), arXiv:2405.04896 [cs.SI] .
  * de Arruda  _et al._ (2021) G. F. de Arruda, G. Petri, P. M. Rodriguez, and Y. Moreno, arXiv preprint arXiv:2112.04273 (2021).
  * Iacopini  _et al._ (2022) I. Iacopini, G. Petri, A. Baronchelli, and A. Barrat, Communications Physics 5, 64 (2022).
  * St-Onge  _et al._ (2022) G. St-Onge, I. Iacopini, V. Latora, A. Barrat, G. Petri, A. Allard, and L. Hébert-Dufresne, Communications Physics 5, 25 (2022).

##  VII Appendix A.
The Random Hypergraph Model

###  VII.1 1\. Asymptotic results: estimation of the number of empty
hyperedges

In the main text, we derived the expression

$p_{0}=(1-p)^{N}=\left(1-\frac{T}{NL}\right)^{N}=\left(1-\frac{h}{N}\right)^{N}$
(52)

describing the probability that a generic hyperedge is empty. It can be
manipulated more easily upon considering the case $N\gg 1$. In fact,

$p_{0}=\left(1-\frac{h}{N}\right)^{N}\to e^{-h}$ (53)

where the notation is meant to indicate $\lim_{N\to+\infty}p_{0}$. Such a
result affects the asymptotic behaviour of quantities like $p_{1}$ and
$\langle N_{0}\rangle$ since

$p_{1}=1-[1-(1-p)^{N}]^{L}=1-[1-p_{0}]^{L}\to 1-[1-e^{-h}]^{L}$ (54)

and

$\langle N_{0}\rangle=L(1-p)^{N}=Lp_{0}\to Le^{-h}.$ (55)

Since $p=T/NL$ and $h\equiv T/L$ remains finite, we are implicitly posing
$T=O(L)$, i.e. $p=O(1/N)$. In words, all quantities of interest remains
finite, allowing us to derive an explicit expression for the filling
threshold. Requiring $\langle N_{0}\rangle=1$ leads to

$Le^{-h}=1$ (56)

which is solved by $h=\ln L$, further inducing $p=h/N=\ln L/N$. Numerical
simulations confirm the correctness of the aforementioned estimations for a
wide range of different combinations of $N$ and $L$ (see Fig. 6).

Figure 7: Trends of $\frac{\langle N_{0}\rangle}{L}=p_{0}\to e^{-h}$ and
$P_{0}\equiv 2\sum_{\alpha<\beta}P(d_{\alpha\beta}=0)/L(L-1)$ as a function of
$T/NL$ under the RHM (left panels) and the HCM (right panels). The dense
(sparse) regime is recovered for large (small) values of $T/NL$. Each dot
represents an average taken over an ensemble of $1000$ configurations
(explicitly sampled from either the RHM or the HCM) and is accompanied by the
corresponding $95\%$ confidence interval.

###  VII.2 2\. Asymptotic results: estimation of the number of parallel
hyperedges

Let us, now, focus on the issue of parallel hyperedges. The position $h\equiv
T/L$ leads to the results

$P(d_{\alpha\beta}=0)=[1-2p(1-p)]^{N}=\left[1-\frac{2T}{NL}\left(1-\frac{T}{NL}\right)\right]^{N}=\left[1-\frac{2h}{N}\left(1-\frac{h}{N}\right)\right]^{N}\to
e^{-2h}$ (57)

and

$\langle
d_{\alpha\beta}\rangle=2p(1-p)N=\frac{2T}{NL}\left(1-\frac{T}{NL}\right)N=\frac{2h}{N}\left(1-\frac{h}{N}\right)N\to
2h;$ (58)

all quantities of interest remain finite, allowing us to derive an explicit
expression for the resolution threshold. Requiring $\langle
d_{\alpha\beta}\rangle=1$ leads to

$2h=1$ (59)

which is solved by $h=1/2$, further inducing $p=h/N=1/2N$. Numerical
simulations confirm the correctness of the aforementioned estimations for a
wide range of different combinations of $N$ and $L$ (see Fig. 7).

##  VIII Appendix B.
The Hypergraph Configuration Model

###  VIII.1 1\. Approximating the Hypergraph Configuration Model

Let us, first, consider that the generic probability coefficient induced by
the HCM can be Taylor-expanded as

$p_{i\alpha}^{\text{HCM}}=\frac{x_{i}y_{\alpha}}{1+x_{i}y_{\alpha}}=x_{i}y_{\alpha}-(x_{i}y_{\alpha})^{2}+(x_{i}y_{\alpha})^{3}\dots$
(60)

one may, thus, wonder at which order the expansion can be safely truncated. In
case hypergraphs are sparse, i.e. $T/NL\ll 1$, and in absence of hubs it turns
out that just considering the first order is enough to obtain predictions as
accurate as those achievable under the ‘full’ HCM: in other words, one can
simply put $p_{i\alpha}^{\text{FOA}}\simeq x_{i}y_{\alpha}$. Upon doing so,
the system of equations defining the HCM simplifies to

$\displaystyle k_{i}^{*}$
$\displaystyle=\sum_{\alpha=1}^{L}x_{i}y_{\alpha},\>\forall\>i$ (61)
$\displaystyle h_{\alpha}^{*}$
$\displaystyle=\sum_{i=1}^{N}x_{i}y_{\alpha},\>\forall\>\alpha$ (62)

expressions leading us to find $x_{i}=k_{i}^{*}/\sqrt{T^{*}}$, $\forall\>i$
and $y_{\alpha}=h_{\alpha}^{*}/\sqrt{T^{*}}$, $\forall\>\alpha$. As a
consequence,
$p_{i\alpha}^{\text{FOA}}=x_{i}y_{\alpha}=k_{i}^{*}h_{\alpha}^{*}/T^{*}$ \- a
position that is commonly named Chung-Lu approximation (CLA). Figure 8 shows
the goodness of the latter one in approximating
$\langle\sigma_{i}\rangle_{\text{HCM}}$ for the contact-primary-school, the
email-Enron and the congress-bills data sets, i.e. the accuracy of the
following chain of equalities

$\displaystyle\langle\sigma_{i}\rangle_{\text{HCM}}$
$\displaystyle=\sum_{\alpha=1}^{L}p_{i\alpha}^{\text{HCM}}[h_{\alpha}-p_{i\alpha}^{\text{HCM}}]\overset{\text{FOA}}{\simeq}\sum_{\alpha=1}^{L}x_{i}y_{\alpha}[h_{\alpha}-x_{i}y_{\alpha}]\overset{\text{CLA}}{=}\frac{k_{i}^{*}}{T^{*}}\left(1-\frac{k_{i}^{*}}{T^{*}}\right)\sum_{\alpha=1}^{L}(h_{\alpha}^{*})^{2}\simeq
k_{i}^{*}\sum_{\alpha=1}^{L}\frac{(h_{\alpha}^{*})^{2}}{T^{*}},\>\forall\>i.$
(63)

Figure 8: Left panels: numerical values of the node-specific Lagrange
multipliers obtained by solving the ‘full’ HCM scattered versus the ones
derived by adopting the CLA. Middle panels: numerical values of the hyperedge-
specific Lagrange multipliers obtained by solving the ‘full’ HCM scattered
versus the ones derived by adopting the CLA. Right panels: expected values
$\\{\langle\sigma_{i}\rangle_{\text{HCM}}\\}$ scattered versus the expected
values $\\{\langle\sigma_{i}\rangle_{\text{CLA}}\\}$.

###  VIII.2 2\. Partial Configuration Models for binary, undirected
hypergraphs

The HCM reduces to a ‘partial’ Configuration Model in case one of the two
degree sequences is left unconstrained.

Hypergraph Partial Configuration Model - layer of nodes. For instance, let us
constrain only the degree sequence $\\{k_{i}\\}_{i=1}^{N}$, by posing
$\beta_{\alpha}=0$, $\forall\>\alpha$ (or, equivalently, $y_{\alpha}=1$,
$\forall\>\alpha$). Our canonical probability distribution becomes

$P(\mathbf{I})=\prod_{i=1}^{N}\prod_{\alpha=1}^{L}\frac{x_{i}^{I_{i\alpha}}}{1+x_{i}}=\prod_{i=1}^{N}p_{i}^{k_{i}}(1-p_{i})^{L-k_{i}}$
(64)

where $e^{-\alpha_{i}}\equiv x_{i}$ and $p_{i}\equiv x_{i}/(1+x_{i})$. The
entries of the incidence matrix are, now, independent random variables obeying
the Bernoulli distributions reading

$I_{i\alpha}\sim\text{Ber}(p_{i}),\>\forall\>i,\alpha;$ (65)

in words, the entries along the same row obey the same distribution while the
ones along the same column obey different distributions. As a consequence, the
$i$-th node degree (a sum of i.i.d. Bernoulli random variables), obeys the
Binomial distribution

$k_{i}\sim\text{Bin}(L,p_{i})$ (66)

while the $\alpha$-th hyperedge degree (a sum of independent random variables
that obey different Bernoulli distributions) obeys the Poisson-Binomial
distribution

$h_{\alpha}\sim\text{PoissBin}(N,p_{1}\dots p_{N}).$ (67)

Interestingly, while the degrees of the nodes obey different Binomial
distributions, the degrees of the hyperedges obey the same Poisson-Binomial
distribution.

The resolution of the likelihood maximisation problem leads us to find the
values $p_{i}=k_{i}^{*}/L$, $\forall\>i$ which, in turn, ensure that $\langle
k_{i}\rangle_{\text{HPCM}_{\text{nodes}}}=\sum_{\alpha=1}^{L}p_{i}=k_{i}^{*}$,
$\forall\>i$ and that $\langle T\rangle_{\text{HPCM}_{\text{nodes}}}=T^{*}$;
instead, $\langle h_{\alpha}\rangle_{\text{HPCM}_{\text{nodes}}}$ will, in
general, differ from $h_{\alpha}^{*}$ \- in fact, $\langle
h_{\alpha}\rangle_{\text{HPCM}_{\text{nodes}}}=T^{*}/L$, $\forall\>\alpha$.

From a microcanonical perspective, the number of configurations satisfying the
requirement that the degrees of the nodes match their empirical values amounts
to

$\Omega_{\text{HPCM}_{\text{nodes}}}=\prod_{i=1}^{N}\binom{L}{k_{i}^{*}};$
(68)

reshuffling the 1s along each row of the incidence matrix separately ensures
the degrees of the nodes to be preserved, while destroying any other
correlation.

Hypergraph Partial Configuration Model - layer of hyperedges. Analogously, the
canonical probability distribution describing the case in which only the
hyperdegree sequence $\\{h_{\alpha}\\}_{\alpha=1}^{L}$ is constrained reads

$P(\mathbf{I})=\prod_{i=1}^{N}\prod_{\alpha=1}^{L}\frac{y_{\alpha}^{I_{i\alpha}}}{1+y_{\alpha}}=\prod_{\alpha=1}^{L}p_{\alpha}^{h_{\alpha}}(1-p_{\alpha})^{N-h_{\alpha}}$
(69)

where $e^{-\beta_{\alpha}}\equiv y_{\alpha}$ and $p_{\alpha}\equiv
y_{\alpha}/(1+y_{\alpha})$. As for the previous null model, the entries of the
incidence matrix are independent random variables obeying different Bernoulli
distributions, i.e.

$I_{i\alpha}\sim\text{Ber}(p_{\alpha}),\>\forall\>i,\alpha;$ (70)

in words, the entries along the same column obey the same distribution while
the ones along the same row obey different distributions. As a consequence,
the $i$-th node degree (a sum of independent random variables that obey
different Bernoulli distributions) obeys the Poisson-Binomial distribution

$k_{i}\sim\text{PoissBin}(L,p_{1}\dots p_{L})$ (71)

while the $\alpha$-th hyperedge degree (a sum of i.i.d. Bernoulli random
variables) obeys the Binomial distribution

$h_{\alpha}\sim\text{Bin}(N,p_{\alpha}).$ (72)

Interestingly, while the degrees of the nodes obey the same Poisson-Binomial
distribution, the degrees of the hyperedges obey different Binomial
distributions.

The resolution of the likelihood maximisation problem leads us to find the
values $p_{\alpha}=h_{\alpha}^{*}/N$, $\forall\>\alpha$ which, in turn, ensure
that $\langle
h_{\alpha}\rangle_{\text{HPCM}_{\text{hyperedges}}}=\sum_{i=1}^{N}p_{\alpha}=h_{\alpha}^{*}$,
$\forall\>\alpha$ and that $\langle
T\rangle_{\text{HPCM}_{\text{hyperedges}}}=T^{*}$; instead, $\langle
k_{i}\rangle_{\text{HPCM}_{\text{hyperedges}}}$ will, in general, differ from
$k_{i}^{*}$ \- in fact, $\langle
k_{i}\rangle_{\text{HPCM}_{\text{hyperedges}}}=\sum_{\alpha=1}^{L}p_{\alpha}=T^{*}/N$,
$\forall\>i$.

From a microcanonical perspective, the number of configurations satisfying the
requirement that the degrees of the hyperedges match their empirical values
amounts to

$\Omega_{\text{HPCM}_{\text{hyperedges}}}=\prod_{\alpha=1}^{L}\binom{N}{h_{\alpha}^{*}};$
(73)

reshuffling the 1s along each column of the incidence matrix separately
ensures the degrees of the hyperedges to be preserved, while destroying any
other correlation.

|  Total number of 1s, $\langle T\rangle$ |  Degrees, $\langle k_{i}\rangle$ |  Hyperdegrees, $\langle h_{\alpha}\rangle$
---|---|---|---
Random Hypergraph Model | $T^{*}$ | $T^{*}/N$ | $T^{*}/L$
Hypergraph Partial Configuration Model (node layer) | $T^{*}$ | $k_{i}^{*}$ | $T^{*}/L$
Hypergraph Partial Configuration Model (hyperedge layer) | $T^{*}$ | $T^{*}/N$ | $h_{\alpha}^{*}$
Hypergraph Configuration Model | $T^{*}$ | $k_{i}^{*}$ | $h_{\alpha}^{*}$
Table 2: Summary of the constraints defining the canonical version of each
benchmark considered here: while the total number of 1s is preserved by each
model, the degrees are separately preserved by the two, partial HCMs and
jointly preserved only by the HCM. Stars indicate the empirical values of the
constraints.

###  VIII.3 3\. Asymptotic results: estimation of the number of empty
hyperedges

As noticed in the main text, the asymptotic behaviour of the HCM depends on
the whole set of coefficients $\\{p_{i\alpha}\\}_{i,\alpha}$: closed form
expressions can be, thus, obtained only after having specified the functional
dependence of each of them on $N$. General conclusions can be drawn more
easily for the two, ‘partial’ Configuration Models. Let us start by
considering the behaviour of HPCM${}_{\text{hyperedges}}$, defined by the
coefficients $p_{i\alpha}=h_{\alpha}/N$, $\forall\>i,\alpha$. Such a position
leads to

$p_{0}^{\alpha}=\prod_{i=1}^{N}(1-p_{i\alpha})=\left(1-\frac{h_{\alpha}}{N}\right)^{N}\to
e^{-h_{\alpha}},$ (74)

that, in turn, induces the expression

$p_{1}=1-\prod_{\alpha=1}^{L}\left[1-\prod_{i=1}^{N}(1-p_{i\alpha})\right]=1-\prod_{\alpha=1}^{L}\left[1-\left(1-\frac{h_{\alpha}}{N}\right)^{N}\right]\to
1-\prod_{\alpha=1}^{L}(1-e^{-h_{\alpha}}),$ (75)

for the probability of observing at least one, empty hyperedge, and the
expression

$\langle
N_{0}\rangle=\sum_{\alpha=1}^{L}\prod_{i=1}^{N}(1-p_{i\alpha})=\sum_{\alpha=1}^{L}\left(1-\frac{h_{\alpha}}{N}\right)^{N}\to\sum_{\alpha=1}^{L}e^{-h_{\alpha}},$
(76)

for the expected number of empty hyperedges. It is interesting to consider the
‘full’ HCM and calculate $p_{0}^{\alpha}$ in the sparse regime. We get

$p_{0}^{\alpha}=\prod_{i=1}^{N}(1-p_{i\alpha})\simeq
1-\sum_{i=1}^{N}p_{i\alpha}=1-h_{\alpha}\simeq e^{-h_{\alpha}},$ (77)

an expression proving that the trends returned by HPCM${}_{\text{hyperedges}}$
are compatible with those returned by the ‘full’ HCM (see also Fig. 9).

Let us, now, inspect the behaviour of HPCM${}_{\text{nodes}}$, defined by the
coefficients $p_{i\alpha}=k_{i}/L$, $\forall\>i,\alpha$. Such a position, in
the sparse regime, leads to

$p_{0}=\prod_{i=1}^{N}(1-p_{i\alpha})=\prod_{i=1}^{N}\left(1-\frac{k_{i}}{L}\right)\simeq
1-\sum_{i=1}^{N}\frac{k_{i}}{L}=1-\frac{T}{L}=1-h\simeq e^{-h}$ (78)

(where the third passage comes from the choice of considering only the leading
order of the product, in the sparse regime, and the last one comes from the
choice of truncating the exponential function to the first order of its Taylor
expansion). The equation above allows us to recover the threshold value $h=\ln
L$ and $p_{f}^{\text{RHM}}=h/N=\ln L/N$. Consistently, both the probability of
observing at least one, isolated hyperedge and the expected number of isolated
hyperedges coincide with the ones returned by the RHM (see also Fig. 9).

Figure 9: Trends of $\frac{\langle N_{0}\rangle}{L}=p_{0}$ and $P_{0}\equiv
2\sum_{\alpha<\beta}P(d_{\alpha\beta}=0)/L(L-1)$ as a function of $T/NL$ under
the ‘partial’ Configuration Models: while the trends of
HPCM${}_{\text{nodes}}$ overlap with those of the RHM, the trends of
HPCM${}_{\text{hyperedges}}$ overlap with those of the HCM. The dense (sparse)
regime is recovered for large (small) values of $T/NL$. Each dot represents an
average taken over an ensemble of $1000$ configurations (explicitly sampled
from either the RHM or the HCM) and is accompanied by the corresponding $95\%$
confidence interval. The light-blue line indicates the density value
corresponding to the resolution threshold. The magenta line indicates the
density value corresponding to the filling threshold. The red line indicates
the density value corresponding to $h=1$. The yellow line indicates the
density value in correspondence of which the minimum, expected hyperdegree
amounts to 1.

###  VIII.4 4\. Asymptotic estimation of the number of parallel hyperedges

Let us, now, focus on the issue of parallel hyperedges. Let us start by
considering the behaviour of HPCM${}_{\text{hyperedges}}$, defined by the
coefficients $p_{i\alpha}=h_{\alpha}/N$, $\forall\>i,\alpha$. Such a position
leads to

$\displaystyle
P(d_{\alpha\beta}=0)=\prod_{i=1}^{N}(1-q_{i}^{\alpha\beta})=\prod_{i=1}^{N}[1-p_{i\alpha}(1-p_{i\beta})-p_{i\beta}(1-p_{i\alpha})]$
$\displaystyle=\prod_{i=1}^{N}\left(1-\frac{h_{\alpha}+h_{\beta}}{N}+\frac{2h_{\alpha}h_{\beta}}{N^{2}}\right)$
$\displaystyle=\left(1-\frac{h_{\alpha}+h_{\beta}}{N}+\frac{2h_{\alpha}h_{\beta}}{N^{2}}\right)^{N}\to
e^{-(h_{\alpha}+h_{\beta})}$ (79)

and to

$\langle
d_{\alpha\beta}\rangle=\sum_{i=1}^{N}q_{i}^{\alpha\beta}=\sum_{i=1}^{N}[p_{i\alpha}(1-p_{i\beta})+p_{i\beta}(1-p_{i\alpha})]=\left(\frac{h_{\alpha}+h_{\beta}}{N}-\frac{2h_{\alpha}h_{\beta}}{N^{2}}\right)N\to
h_{\alpha}+h_{\beta}.$ (80)

It is interesting to consider the ‘full’ HCM and calculate
$q_{i}^{\alpha\beta}=p_{i\alpha}(1-p_{i\beta})+p_{i\beta}(1-p_{i\alpha})$ in
the sparse regime, i.e. by truncating it to the first addendum: in words,
considering $q_{i}^{\alpha\beta}\simeq p_{i\alpha}+p_{i\beta}$ allows us to
recover the trends returned by HPCM${}_{\text{hyperedges}}$ (see also Fig. 9).

Let us, now, inspect the behaviour of HPCM${}_{\text{nodes}}$, defined by the
coefficients $p_{i\alpha}=k_{i}/L$, $\forall\>i,\alpha$. Such a position, in
the sparse regime, leads to

$\displaystyle
P(d_{\alpha\beta}=0)=\prod_{i=1}^{N}(1-q_{i}^{\alpha\beta})=\prod_{i=1}^{N}[1-p_{i\alpha}(1-p_{i\beta})+p_{i\beta}(1-p_{i\alpha})]$
$\displaystyle=\prod_{i=1}^{N}\left[1-\frac{2k_{i}}{L}\left(1-\frac{k_{i}}{L}\right)\right]$
$\displaystyle\simeq
1-\sum_{i=1}^{N}\frac{2k_{i}}{L}=1-\frac{2T}{L}=1-2h\simeq e^{-2h},$ (81)

(where the fourth passage comes from the choice of considering only the
leading order of the product, in the sparse regime, and the last one comes
from the choice of truncating the exponential function to the first order of
its Taylor expansion) and to

$\langle
d_{\alpha\beta}\rangle=\sum_{i=1}^{N}q_{i}^{\alpha\beta}=\sum_{i=1}^{N}[p_{i\alpha}(1-p_{i\beta})+p_{i\beta}(1-p_{i\alpha})]=\sum_{i=1}^{N}\frac{2k_{i}}{L}\left(1-\frac{k_{i}}{L}\right)\simeq\sum_{i=1}^{N}\frac{2k_{i}}{L}=\frac{2T}{L}=2h$
(82)

(where the fourth passage comes from the choice of considering only the
leading order of the sum, in the sparse regime, and the last one comes from
the choice of truncating the exponential function to the first order of its
Taylor expansion). The equation above allows us to recover the threshold value
$h=1/2$ and $p_{r}^{\text{RHM}}=h/N=1/2N$ (see also Fig. 9).

##  IX Appendix C.
Expected value of topological quantities

Let us, now, provide the explicit expression of the expected value of some
topological properties of interest. Let us start with the entries of the
matrix $\mathbf{W}$, for which the following results hold true

$\displaystyle\langle w_{ij}\rangle$ $\displaystyle=\sum_{\alpha=1}^{L}\langle
I_{i\alpha}I_{j\alpha}\rangle=\sum_{\alpha=1}^{L}\langle
I_{i\alpha}\rangle\langle
I_{j\alpha}\rangle=\sum_{\alpha=1}^{L}p_{i\alpha}p_{j\alpha},$ (83)
$\displaystyle\text{Var}[w_{ij}]$
$\displaystyle=\sum_{\alpha=1}^{L}\text{Var}[I_{i\alpha}I_{j\alpha}]+2\sum_{\beta<\gamma}\text{Cov}[I_{i\beta}I_{j\beta},I_{i\gamma}I_{j\gamma}]=\sum_{\alpha=1}^{L}p_{i\alpha}p_{j\alpha}(1-p_{i\alpha}p_{j\alpha}).$
(84)

Therefore,

$\displaystyle\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\langle w_{ij}\rangle^{2}$
$\displaystyle=\sum_{\alpha=1}^{L}p_{i\alpha}^{2}\left[\sum_{j=1}^{N}p_{j\alpha}^{2}-p_{i\alpha}^{2}\right]+2\sum_{\beta<\gamma}p_{i\beta}p_{i\gamma}\left[\sum_{j=1}^{N}p_{j\beta}p_{j\gamma}-p_{i\beta}p_{i\gamma}\right],$
(85) $\displaystyle\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\text{Var}[w_{ij}]$
$\displaystyle=\sum_{\alpha=1}^{L}p_{i\alpha}(\langle h_{\alpha}\rangle-
p_{i\alpha})-\sum_{\alpha=1}^{L}p_{i\alpha}^{2}\left[\sum_{j=1}^{N}p_{j\alpha}^{2}-p_{i\alpha}^{2}\right].$
(86)

Figure 10: Left panels: scatter plots between the analytical and the sample
estimates of the numerator of the disparity ratio, under the ‘full’ HCM.
Middle panels: scatter plots between the analytical and the sample estimates
of the denominator of the disparity ratio, under the ‘full’ HCM. Right panels:
scatter plots between the expected values $\\{\langle
Y_{i}\rangle_{\text{HCM}}\\}$, computed as in Eq. (89), and the ones obtained
by explicitly sampling the ensemble induced by the ‘full’ HCM. As analytical
and sampled values of numerators and denominators nicely agree, the result
concerning $\\{\langle Y_{i}\rangle\\}$ must be imputed to the correlation
between the two, disregarded by Eq. (89).

Let us, now, consider the strength sequence of matrix $\mathbf{W}$, for which
the following results hold true

$\displaystyle\langle\sigma_{i}\rangle$
$\displaystyle=\sum_{\alpha=1}^{L}p_{i\alpha}(\langle h_{\alpha}\rangle-
p_{i\alpha}),$ (87) $\displaystyle\text{Var}[\sigma_{i}]$
$\displaystyle=\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\text{Var}[w_{ij}]+2\sum_{\begin{subarray}{c}j<k\\\
(j,k\neq
i)\end{subarray}}\text{Cov}[w_{ij},w_{ik}]=\sum_{\begin{subarray}{c}j=1\\\
(j\neq
i)\end{subarray}}^{N}\text{Var}[w_{ij}]+2\sum_{\begin{subarray}{c}j<k\\\
(j,k\neq i)\end{subarray}}[\langle w_{ij}w_{ik}\rangle-\langle
w_{ij}\rangle\langle w_{ik}\rangle].$ (88)

As a consequence, the expected value of the disparity ratio reads

$\displaystyle\langle Y_{i}\rangle$
$\displaystyle\simeq\sum_{\begin{subarray}{c}j=1\\\ (j\neq
i)\end{subarray}}^{N}\frac{\langle
w_{ij}^{2}\rangle}{\langle\sigma_{i}^{2}\rangle}=\sum_{\begin{subarray}{c}j=1\\\
(j\neq i)\end{subarray}}^{N}\frac{\langle
w_{ij}\rangle^{2}+\text{Var}[w_{ij}]}{\langle\sigma_{i}\rangle^{2}+\text{Var}[\sigma_{i}]}$
(89)

where the $\simeq$ symbol is understood to approximate the expected value of a
ratio as the ratio of expected values. This approximation just represents the
first order of the Taylor expansion

$\text{E}\left[\frac{X}{Y}\right]\simeq\frac{\text{E}[X]}{\text{E}[Y]}-\frac{\text{Cov}[X,Y]}{\text{E}[Y]^{2}}+\frac{\text{E}[X]}{\text{E}[Y]^{3}}\text{Var}[X];$
(90)

due to the correlation between numerator and denominator, the second and the
third terms of Eq. (90) may not be negligible: as Fig. 10 shows, while Eq.
(89) holds true for the contact-primary-school data set, it does not for the
others.

Figure 11: Scatter plot between the ensemble average (of each component) of
the Perron-Frobenius eigenvector of a given, empirical matrix $\mathbf{W}^{*}$
and the numerical value (of each component) of the Perron-Frobenius
eigenvector of the corresponding, expected matrix
$\langle\mathbf{W}\rangle_{\text{HCM}}$. The extremely good agreement suggests
a faster way to calculate the CEC that avoids to explicitly sample the
ensemble induced by a null model.

##  X Appendix D.
Expected value of the clique motif eigenvector centrality

The clique motif eigenvector centrality (CEC) has been proposed in Benson
(2019): its $i$-th entry, $\text{CEC}_{i}$, is nothing else that the
corresponding entry of the Perron-Frobenius eigenvector of $\mathbf{W}$.

In order to evaluate the expected value of the CEC, we have explicitly sampled
the ensemble of incidence matrices induced by the HCM; then, we have
calculated the matrix $\tilde{\mathbf{W}}$ induced by each
$\tilde{\mathbf{I}}$, according to Eq. (35). Afterwards, we have calculated
the Perron-Frobenius eigenvector of each ‘projected’ matrix and taken their
average in an entry-wise fashion.

Remarkably, we found that the aforementioned ensemble average basically
coincides with the Perron-Frobenius eigenvector of the ensemble average of
$\mathbf{W}$ itself, i.e. $\langle\mathbf{W}\rangle_{\text{HCM}}$, as Fig. 11
shows.

##  XI Appendix E.
Expected value of the entries of the confusion matrix

Before defining the entries of the confusion matrix, let us call
$\tilde{\mathbf{I}}$ the generic member of the ensemble of incidence matrices
induced by one of the benchmarks considered in the present paper: the expected
value of any quantity which is a function of $\tilde{\mathbf{I}}$ is readily
calculated by averaging it over the ensemble itself.

The true positive rate (TPR) is defined as the percentage of 1s correctly
recovered by a given model and reads

$\text{TPR}=\frac{\text{TP}}{T}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{I_{i\alpha}\tilde{I}_{i\alpha}}{T}\Longrightarrow\langle\text{TPR}\rangle=\frac{\langle\text{TP}\rangle}{T}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{I_{i\alpha}p_{i\alpha}}{T};$
(91)

the specificity (SPC) is defined as the percentage of 0s correctly recovered
by a given model and reads

$\text{SPC}=\frac{\text{TN}}{NL-T}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{(1-I_{i\alpha})(1-\tilde{I}_{i\alpha})}{NL-T}\Longrightarrow\langle\text{SPC}\rangle=\frac{\langle\text{TN}\rangle}{NL-T}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{(1-I_{i\alpha})(1-p_{i\alpha})}{NL-T}$
(92)

the positive predictive value (PPV) is defined as the percentage of 1s
correctly recovered by a given model with respect to the total number of
predicted of 1s and reads

$\text{PPV}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{I_{i\alpha}\tilde{I}_{i\alpha}}{\tilde{T}}\Longrightarrow\langle\text{PPV}\rangle=\frac{\langle\text{TP}\rangle}{\langle
T\rangle}=\sum_{i=1}^{N}\sum_{\alpha=1}^{L}\frac{I_{i\alpha}p_{i\alpha}}{\langle
T\rangle};$ (93)

lastly, the accuracy (ACC) measures the overall performance of a given
reconstruction method in correctly placing both 0s and 1s and reads

$\text{ACC}=\frac{\text{TP}+\text{TN}}{NL}\Longrightarrow\langle\text{ACC}\rangle=\frac{\langle\text{TP}\rangle+\langle\text{TN}\rangle}{NL}.$
(94)

As stressed in the main text, results on the confusion matrix of a number of
real-world hypergraphs reveal that the large sparsity of the latter ones makes
it difficult to reproduce the TPR and the PPV (see Table 3); on the other
hand, the capability of the HCM to reproduce the density of 1s makes it
capable of recovering the density of 0s as well, thus ensuring the overall
accuracy of the model to be quite large.

Data set | $\langle\text{TPR}\rangle$ | $\langle\text{SPC}\rangle$ | $\langle\text{PPV}\rangle$ | $\langle\text{ACC}\rangle$ | $T/NL$
---|---|---|---|---|---
contact-primary-school | $3.57\times 10^{-8}$ | 0.991 | $3.57\times 10^{-8}$ | 0.983 | $8.66\times 10^{-3}$
email-Enron | $8.15\times 10^{-8}$ | 0.982 | $8.15\times 10^{-8}$ | 0.966 | $1.72\times 10^{-2}$
NDC-classes | $2.40\times 10^{-10}$ | 0.997 | $2.40\times 10^{-10}$ | 0.995 | $2.71\times 10^{-3}$
Table 3: Confusion matrix and density of 1s for the three real-world
hypergraphs considered in the present paper.

##  XII Appendix F.
Comparing the RHM and the HCM performance

Let us, now, inspect the performance of the RHM and that of the HCM in
reproducing the patterns of our real-world hypergraphs. To this aim, let us
consider the coefficient of determination that, in our case, can be re-defined
as

$R^{2}_{\text{null}}(\mathbf{X})=1-\frac{\sum_{i}\left(X_{i}-\langle
X_{i}\rangle_{\text{null}}\right)^{2}}{\sum_{i}\left(X_{i}-\bar{X}\right)^{2}};$
(95)

it measures the goodness of a null model in reproducing empirical data when
compared to a baseline model whose only prediction reads $\bar{X}$, i.e. the
arithmetic mean of the empirical data themselves. While such a model is
characterised by $R^{2}=0$, a model that is capable of matching the observed
values exactly is characterised by $R^{2}=1$; models returning predictions
whose discrepancies are larger than the discrepancies accompanying the
predictions of the baseline model are characterised by $R^{2}<0$.

As Table 4 shows, the HCM is often more accurate than the RHM in reproducing
the quantities we have considered. More precisely, the RHM is found to
steadily perform worse than it when the CEC is considered and to perform quite
similar to it when the $Y$ is considered; the opposite tendency is observed
when the $\sigma$ is considered. Overall, these results confirm that just
constraining the total number of 1s of a given incidence matrix is not enough
to achieve an accurate description of the corresponding hypergraph but also
suggest that constraining the degree and the hyperdegree sequences may be not
enough as well - as we have already noticed when commenting the highly non-
trivial degree of self-organisation of the real-world hypergraphs (whose
trends are) depicted in Fig. 4.

Data set | $R^{2}_{\text{HCM}}(\sigma)$ | $R^{2}_{\text{HCM}}(Y)$ | $R^{2}_{\text{HCM}}(\text{CEC})$ | $R^{2}_{\text{RHM}}(\sigma)$ | $R^{2}_{\text{RHM}}(Y)$ | $R^{2}_{\text{RHM}}(\text{CEC})$
---|---|---|---|---|---|---
contact-primary-school | -2.9697 | -1.7248 | -0.1489 | -2.7349 | -1.8189 | -0.8478
email-Enron | -0.1344 | -0.2591 | -0.0519 | -0.0131 | -0.3462 | -0.1705
NDC-classes | -0.0044 | -0.1107 | -0.0482 | -0.0008 | -0.1222 | -0.0183
Table 4: The $R^{2}$ index allows us to compare the performance of the RHM
with that of the HCM in reproducing the patterns of our real-world
hypergraphs. While both models perform worse than the baseline model, the HCM
is, overall, more accurate than the RHM in reproducing all quantities we have
considered, the only exception being $\sigma$: therefore, just constraining
the total number of 1s of a given incidence matrix is, in general, not enough
to achieve an accurate description of the hypergraph it represents - although
constraining the degree and the hyperdegree sequences may be not enough as
well.