Hypergraph reconstruction from network data

Networks can describe the structure of a wide variety of complex systems by specifying which pairs of entities in the system are connected. While such pairwise representations are flexible, they are not necessarily appropriate when the fundamental interactions involve more than two entities at the same time. Pairwise representations nonetheless remain ubiquitous, because higher-order interactions are often not recorded explicitly in network data. Here, we introduce a Bayesian approach to reconstruct latent higher-order interactions from ordinary pairwise network data. Our method is based on the principle of parsimony and only includes higher-order structures when there is sufficient statistical evidence for them. We demonstrate its applicability to a wide range of datasets, both synthetic and empirical.


I. INTRODUCTION
Empirical networks are often locally dense and globally sparse [1]. Whether they are social, biological, or technological [2], they comprise large groups of densely interconnected nodes, even when only a small fraction of all possible connections exist. This situation leads to delicate modeling challenges: How can we account for two seemingly contradictory properties of networks-density and sparsity-in our models?
Abundant prior work going back to the early days of social network analysis [3,4] and network science [5,6] suggests that higher-order interactions [7] are a possible explanation for the local density of networks [1,6]. According to this reasoning, entities are connected because they have a shared context-a higher-order interactionwithin which connections can be created [8]. It is clear that a phenomenon along these lines occurs in many social processes: scientists appear as collaborators in the Web of Science because they co-author papers together; colleagues exchange emails because they are part of the same department or the same division of a company. It is also known that similar phenomena explain tie formation in a broader range of networked systems, including biological, technological or informational systems [7].
The ubiquity of higher-order interactions provides a simple and universal explanation for the observed structure of empirical networks. If we assume that most ties are created within contexts of limited scopes, then the resulting networks are locally dense, matching empirical observations [1,9,10].
Despite their tremendous explanatory power, higherorder interactions are seldom used directly to model empirical systems, due to a lack of data [7]. Indeed, while Projected higher-order interactions. We show the extended ego network (circle nodes) of a participant (square node) of the AddHealth study [14]. Friendships are measured between pairs of participants (links and nodes respectively), even when the fundamental units are groups of friends [12]. Multiple combinations of groups and isolated friendships lead to the same network (gray arrows). the context is directly observable for some systems-say, co-authored papers or co-locating species-it is unavailable for several others, including brain data [11], typical social interaction data [12], and ecological competitor data [13] to name only a few.
As a specific motivating example, consider one of the empirical social networks gathered as part of the US National Longitudinal Study of Adolescent to Adult Health [14]. This dataset is constructed using surveys, where participants are asked to nominate their friends. Even though there are good reasons to believe that people often interact because of higher-order groups [12], the survey cannot reveal these groups as it only inquires about pairwise relationships. If we actually need the higher-order interactions to give an appropriate description of the social dynamics at play [12], what should we do with such inadequate survey data? As we show in Fig. 1, there are many kinds of higher-order interactions that are compatible with the same network data. How can one pick among all these possible higher-order descriptions?
Prior work on higher-order interaction discovery in network data often uses cliques-fully connected subgraphs-to identify the interactions [15][16][17]. Cliquebased methods are straightforward to implement because they rely on clique enumeration, a classical problem for which we have exact [18,19] and sampling [20] algorithms that work well in practice. However, clique decompositions do not offer a satisfactory solution to the recovery problem alone. Networks typically admit many possible clique decompositions, which begs the question of which one to pick. For example, a triangle can be decomposed as a single 2-clique, or as three 1-cliques (i.e., as edges); see Fig. 1. In general, the multiplicity of possible solutions implies that higher-order interaction recovery is an ill-posed inverse problem. It becomes well-posed only once we add further constraints on what constitutes a good solution. Thus, existing approaches have sought to address the ill-posed nature of the higher-order interaction recovery problem in various indirect ways. For instance, in graph theory, it is customary to look for a minimal set of cliques covering the network [21,22]. Other methods appeal to notions of randomness and generative modeling to regularize the problem [1,[23][24][25]. These methods describe an explicit process by which one goes from higher-order data to networks, and can therefore assign a likelihood to possible higher-order data representations, allowing the user to single out representations.
In the present work, we develop a Bayesian method for the inference of higher-order of interactions from network. Given a network as input, the method identifies the parts of the network best explained by latent higher-order interactions. Our approach is based on the principle of parsimony and directly addresses the ill-posedness of the reconstruction problem with the methods of information theory. We show that the method can finds compact descriptions of many empirical networked systems by using latent higher-order interactions, thereby demonstrating that such interactions are in complex systems.

A. Generative model
The problem we solve is illustrated in Fig. 1. We have a system we believe is best described with higher-order interactions, but we can only view its structure through the lens of pairwise measurements (an undirected and simple network G); our goal is to reconstruct these higher-order interactions from G only.
For convenience, we encode the higher-order interactions with a hypergraph H [26]. We represent a higherorder interaction between a set of k nodes i 1 , .., i k with a hyperedge of size k. Empirical data often contain repeated interactions between the same group of nodes, so we use hypergraphs with repeated hyperedges and encode the number of hyperedges connecting nodes i 1 , .., i k as A i1,..,i k ≥ 0.
Our method then makes use of a Bayesian generative model to deduce one such hypergraph H from some network dataset G. This generative model gives an explicit description of how the network data G is generated when there are latent higher-order interactions H. With a generative model in place, we can compute the posterior probability that the latent hypergraph is H, given the observed network G. In this equation, P (G|H) and P (H) define our generative model for the data, and its evidence P (G) = H P (G|H)P (H) functions as a normalization constant.
The appeal of such a Bayesian generative formulation is that we can use P (H|G) to make queries about the hypergraph H. What was the most likely set of higherorder interactions? What is the probability that a particular interaction was present in H based on G? How large were the latent higher-order interactions? All of the queries can be answered by computing appropriate averages over P (H|G). As is made evident by Eq. (1), however, we first have to introduce two probability distributions so that we may compute P (H|G) at all. We now define these distributions in detail.

Projection component
The first distribution, P (G|H), is called the projection component of the model. It tells us how likely a particular network G is when the latent hypergraph H is known.
We use a direct projection component and deem two nodes connected in G if and only if these nodes jointly appear in any of the hyperedges of H.
This modeling choice is broadly applicable. For instance, when researchers measure the functional connectivity of two brain regions, they record a connection irrespective of whether the regions peaked as a pair or as jointly with many other regions. Likewise, surveyed social networks contain records of friendships that can be attributed to interactions between pairs of individuals, and to interactions that arise from larger groups.
Certain authors use more nuanced projection components [1,27] and do not assume that the joint participation of two nodes in a hyperedge necessarily leads to measured pairwise interactions. As we have argued in the introduction, we believe that such components blur the lines between community detection and higher-order interaction reconstruction. Hence, we treat measurement as a separate issue to be handled with the methods of Refs. [28,29], for instance.
We formalize the projection component as follows. We set P (G|H) = 1 only when (i) each pair of nodes con-

FIG. 2.
Encoding hypergraphs as factor graphs. The two panels show the factor graph encoding of two different hypergraphs that project to the same network, shown at the bottom of the figure. (a) A hypergraph without higher-order interactions is obtained by associating each edge in G (blue lines) to a node in the factor graph (empty circles); this correspondence is illustrated with a dashed line. Each node in the factor graph is then connected to a factor Aij (squares) where (i, j) is the corresponding edge in G. Higher-order factors A ijk corresponding to possible hyperedges of 3 nodes are also added, and connected to the edges (i, j), (j, k), (i, k) they comprise. Since there is no 4-clique in the graph, the construction stop at this step. Particular combination of hyperedges can then be specified by activating some factors (coloring them blue), here all the factors corresponding to the edges. (b) To encode a hypergraph H with one higher-order interaction, we mark the factor A123 as active.
nected by an edge in G appears jointly in at least one hyperedge of H, and (ii) no two disconnected nodes of G appear together in a hyperedge of H. If either of these conditions is violated, then we set P (G|H) = 0. We can express this definition mathematically as where we use G(H) to denote the projection of H and use G = G(H) to say that H projects to G, or equivalently that (i) and (ii) hold.
Testing G ? = G(H) might appear unwieldy at first but, thankfully, a factor graph encoding of H can help us compute the projection component efficiently by highlighting existing relationships between the edges and cliques of G [30].
To construct this factor graph, we begin by creating two separate sets of nodes: one representing the edges of G, and one representing the cliques of G. Crucially, the second set contains a node for every clique of G, even the included ones like the edges of a triangle, the triangles of a 4-clique, and so on. We call this set the set of factors and refer to nodes in the first set simply as nodes. We obtain a factor graph, by connecting a factor and a node when the corresponding clique contains the corresponding edge.
This construction is illustrated in Fig. 2 for a simple graph of 5 nodes. In Fig. 2, we see that, for example, the edge between nodes 1 and 2 is part of the triangle {1, 2, 3} in G, and it is therefore connected to the factor A 123 . This edge is also part of the 2-clique {1, 2} so it is connected to the factor A 12 , too.
The resulting factor graphs can encode particular hypergraphs H by assigning integers to the factors, corresponding to the number of times every hyperedge appears in H. It is straightforward to check whether G = G(H) holds with this encoding. The first condition-all the connected nodes of G are connected by at least one hyperedge in H-can be verified by checking that every node of the factor graph is connected to at least one active factor, defined as A i1,..,i k > 0. The second conditionno pairs of disconnected nodes in G are connected by a hyperedge of H-is always satisfied by construction, because no factor connects two disconnected nodes of G, so we never represent these forbidden hyperedges with our factor graph.
We note that the factor graph can be stored relatively efficiently, by first enumerating the maximal cliquescliques not included in larger cliques-and then constructing an associative array indexed by cliques, which we expand only when included cliques are needed. Even though enumerating maximal clique is technically an NP-hard problem [31], state-of-the-art enumeration algorithms tend to work well on sparse empirical network data [18,19,32], and indeed we have found that enumeration is not problematic in our experiments.

Hypergraph prior
The second part of Eq. (1), P (H), is the hypergraph prior. Empirical hypergraphs generally have a few properties that a reasonable prior should account for [33]: the size of interactions vary; some of these interactions are repeated, and not all nodes are connected by a hyperedge. It turns out that an existing model [34], known as the Poisson Random Hypergraphs Model (PRHM), reproduces all of these properties. Hence, we adopt it as our hypergraph prior. The PRHM was initially developed to study critical phenomena in hypergraphs [34]; here, we use it to make posterior inferences about networks.
In a nutshell, the PRHM stipulates that the number of hyperedges connecting a set of nodes is a random variable, whose mean λ k only depends on the size k of the set. The variable follows a Poison distribution, such that the number of hyperedges connecting the nodes i 1 , .., i k equals to A i1,..,i k with probability where A i1,..,i k is invariant with respect to permutation of the indexes. The PRHM also models all the hyperedges as independent. Hence, the probability of a particular hypergraph can be calculated as where L is the maximal hyperedge size, C N k denotes all possible subsets of size k of {1, ..., N }, and where λ refers to all the rates collectively. Equation (4) expresses the probability of H in terms of individual hyperedges. To obtain a simpler form, we notice that the number E k of hyperedges of size k, can be calculated as and that there are precisely N k terms in the product over all sets of nodes of size k. We can use these simple observations to rewrite Eq.(4) as where we have defined and where η (k) m is the number of hyperedges of size k that are repeated precisely m times.
In this form, it is clear that the parameters λ control the density of H at all scales. Hence, they more or less determine the kind of hypergraphs we expect to see a priori, and therefore have a major effect on the model output. How can we choose these important parameters carefully?
We propose to a hierarchical empirical Bayes approach, in which we treat λ as unknowns themselves drawn from prior distributions. We use a maximum entropy, or least informative, prior for λ, because we have no information whatsoever about λ a priori. The only thing we know is that these parameters take values in [0, ∞) and are modeled with a finite mean [34]. Hence, the maximal entropy prior of interest is the exponential distribution of mean ν k . We obtain a complete hyperprior for λ by using independent priors for all sizes k, P (λ|ν) = L k=2 P (λ k |ν k ). Integrating over the support of λ, we find that the prior for H is now with ν fixed. It might appear that we have only pushed our problem further ahead-we got rid of λ but we now have a whole new set of parameters on our hands. Notice, however, that the new parameters ν do not have as direct an effect on H. A whole range of densities is now compatible with any choice of ν. And as a result, the model can assign significant probabilities to hypergraphs that project to networks of the correct density, even when the hyperprior is somewhat in error. Hence, we safely fix the new parameters ν with empirical Bayes without risking strongly biased results.
With these precautions in place, we use the observed number of edges E in the network G to choose ν. Our strategy is to equate E to the expected number of edges E(ν) in the network G(H) obtained by projecting H drawn from P (H|ν). This expected density can be approximated as by assuming that hyperedges do not overlap on average. To set the individual values of ν k , we further require that all sizes contribute equally to the final density, with ν k N k = µ for some constant µ. Substituting these equalities in Eq. (10), we obtain such that the prior in Eq. (9) becomes which is the equation we will use henceforth, with µ = E/(L − 1).

B. Properties of the posterior distribution
The model defined in Eqs.
The first noteworthy property is that the model assigns a higher posterior probability to hypergraphs without repeated hyperedges, even though the prior P (H) allows for duplicates. An explicit calculation of how P (H|G) scales with the number of duplicated hyperedges can illustrate this fact. Indeed, consider a hypergraph H 0 with no repeated hyperedges, for which P (G|H 0 ) = 1. Write as α the fraction of k-cliques connected by a hyperedge in H 0 , and consider an experiment in which an average of β ≥ 0 additional hyperedges are placed on top of the hyperedges of size k already present in H 0 . In these hypergraphs, the expected number of hyperedges of size k see Eq. (7). Substituting our various formula in the logarithm of P (H), and using the Stirling approximation log n! ≈ n log n − n, we find that This equation tells us that the log-posterior log P (H|G) decreases with growing β, because the argument of the logarithm is greater or equal to one. Furthermore the likelihood equals one by construction, which implies that the scaling of the prior determines the scaling of the posterior. Hence, the hypergraphs H generated by adding duplicated hyperedges to H 0 -that is by increasing βare less likely than H 0 . A second noteworthy property of the model is that it favors sparser hypergraphs: as long as P (G|H) = 1, the fewer hyperedges, the better. To make this observation precise, suppose we have a hypergraph H m that can be termed minimal for G: every edge of G is covered by exactly one hyperedge of H m and no more. We observe that we cannot improve on the posterior probability of H m by adding a hyperedge, even when this new hyperedge does not fully repeat an existing one. Indeed, consider the hypergraph H m created by adding a hyperedge of size k to H m . (For example, we could add a hyperedge of size 3 on a triangle whose sides were already covered by edges, but did not yet participate in any larger hyperedge together.) By direct calculation, the ratio of posterior probability for H m and H m equals where Z k is the quantity in Eq. (7) for the modified minimal hypergraph, and Z k is the same quantity for the minimal hypergraph. One can show that this ratio is always smaller than one and that, as a result, adding a spurious hyperedge to a minimal hypergraph decreases the posterior probability. The proof is straightforward and relies on the observation that for a minimal hypergraph, we have E k ≤ N k , Z k = 1, and Z k = 1 or Z k = 2. The result follows by direct computation when E k < N k and uses the fact that that Z k = 2 when E k = N k (because adding a single hyperedge to a completely connected minimal hypergraph means one has to double-up one hyperedge).
As a corollary of the two above observations, we conclude that the minimal hypergraphs are high-quality local maxima of P (H|G). We cannot simply pick one of these optima as our reconstruction, however, because there may exist multiple ones of comparable quality. Further, non-optimal hypergraphs may account for a significant fraction of the posterior probability in principle. Instead, we handle these possibly conflicting descriptions by combining them.

C. Posterior estimation
In the Bayesian formulation of hypergraph inference, estimating a given quantity of interests always amount to computing an expectations over the posterior distribution P (H|G). For example, the expected number of hyperedges of size k can be computed as E k = H E k (H)P (H|G). More generally, we are interested in averages of the form for arbitrary functions f that map hypergraphs to vectors or scalars. The summation in Eq. (13) is unfortunately intractable: the set of possible hypergraphs grows exponentially in size with both the number of nodes and the maximal size of the hyperedges. Hence, we propose a Markov Chain Monte-Carlo (MCMC) algorithm to evaluate Eq. (13). This kind of approach generates a random walk over the space of all hypergraphs, with a limiting distribution identical to P (H|G). We use the Metropolis-Hastings (MH) construction to implement the random walk. As is usual, the algorithm consists of proposing a move from H to H with probability Q(H ← H ) and accepting it with probability [35] a = min 1, We use the factor graph representation of H to define these Monte Carlo moves this encoding facilitates checking the value of P (G|H ). Hence, we can state the moves as modifications to the value of the factors A i1,....i k , i.e., the number of hyperedges connecting particular sets of nodes.
The specific set of moves we use goes as follow. For every move, we begin by choosing a maximal factor node uniformly at random from the set of all such factors. We select a size uniformly at random from {2, 3, ..., k} where k is the size of the clique corresponding to the current maximal factor. Then we select one of the subfactors A i1,..,i of size uniformly at random, among the k factors of that size, and we update the selected factors as either A i1,..,i = A i1,..,i +1 or A i1,..,i = A i1,..,i −1 (with probability 1/2). If A i1,..,i was already equal to zero, we force A i1,..,i = A i1,..,i + 1. Therefore we have that Finally, we check whether P (G|H) = 1 using the factor representation, and compute the ratio P (H )/P (H) to obtain the acceptance probability a. We test for acceptance and, if the move is accepted, we record the update. Otherwise, we do nothing. The posterior distribution is rugged so the initialization of the MCMC algorithm matters a great deal in practice. Building on our observations about the properties of P (H|G), we select as our initialization the hypergraph with one hyperedge for every maximal clique of G. This starting point is not a known optimum of P (H|G) but it is close to many of them. Hence, chains initialized at this point have a fairly good chance of converging to a good optimum. And indeed, in our experiments, we find that the maximal clique initialization works much better than a random initialization, an edge initialization or an empty one.

D. Recovery of planted higher-order interactions in synthetic data
To develop an intuition for the workings of our method, we first use our algorithm to uncover higher-order interactions in synthetic data generated by the model appearing in Eq. (2)- (12), altered slightly to facilitate the interpretation of the results. In this experiment, we create a hypergraph that comprises a few large disconnected hyperedges, and we add several random edges (chosen uniformly from the set of all edges) to create a noisy hy-pergraphH. We then project this noise hypergraph to obtain a network G(H), which we feed to our recovery algorithm as input. Our goal in this experiment is to find the hypergraph H * that maximizes the posterior probability P (H|G(H)) (we do not use the full samples given by our MCMC algorithm just yet). We can consider the experiment successful if H * contains all the higher-order interactions planted inH.
The results of this experiment are reported in Fig. 3. At the bottom of Fig. 3a, we show a typical example of what the projected networks G(H) look like when there are very few added random edges. In this regime, the recovered higher-order interactions (in blue) correspond perfectly to those planted inH, for the entire range of parameters we investigated. For the sake of comparison, we also generate an equivalent random network, obtained by completely rewiring the edges of G(H), see the top of Fig. 3a. (Equivalently, we generate an Erdős-Rényi graph with an equal number of edges [36].) This network has the same number of edges as G(H) but is otherwise unstructured. As expected, we find no higher-order interactions beyond the random triangles that occur at this density [37].
If we add many more random edges, we obtain the results shown in Fig. 3c. Again, we can recover the planted higher-order interactions, but we also start to find additional ones, due to the appearance of random triangles formed by triplets of edges added at random [38]. To understand this behavior we turn to the minimum description length (MDL) interpretation of Bayesian inference [39,40].
In a nutshell, the description length is the number of bits that a receiver and a sender with shared knowledge of the model P (G, H) would need to communicate the network G to one another. This communication costs can be minimized by finding a hypergraph H * that is both cheap to communicate and projects to G; receivers who know P (G, H) also know that they can project H * to find G = G(H) * . From this communication perspective, hypergraphs with as few hyperedges as possible are good candidates because they are cheaper to send [24]. The connection with Bayesian inference is that H * happens to coincide with the hypergraph which maximizes the posterior probability P (H|G) (see Supplementary Notes 1 and 2 for a detailed discussion). Hence, maximum a posteriori inference is equivalent to compression.
Reviewing our experiment with the MDL interpretation in mind illuminates the results. In Fig. 3b, we plot the description length provided by our model, for levels of randomness that interpolate between the easy regime shown in Fig. 3a, and the much more random regime appearing in Fig. 3c. We find that the model compresses those networks that have planted interactions much better than their randomized equivalents. These results make intuitive sense: networks with planted interactions contain large cliques, and these cliques can be harnessed to communicate regularities in G. As can be expected, these savings disappear once the large cliques are destroyed by rewiring.

E. Recovery of planted higher-order interactions in empirical data
Having verified that the method works when the higher-order network possess little structure beyond disjoint planted cliques, we turn to more complicated problems. We ask: can our method identify relevant higherorder interactions when the data is (plausibly) more structured? To answer this question, we use empirical bipartite networks and create higher-order networks, by representing the bipartite networks as hypergraphs H [2]. We then project the hypergraphs with Eq. (2), and attempt to recover the planted higher-order interactions in G(H) with our method.
In Fig. 4, we report the results of this experiment for 11 hypergraph constructed with empirical networks datasets [41][42][43][44][45][46][47][48][49][50]. (See also Supplementary Table 1 for a detailed numerical account of the results.) The figure depicts the accuracy of the reconstruction, as quantified a Jaccard similarity J defined as the number of hyperedges found in both the original and the reconstructed hypergraph, divided by the number of hyperedges found in either of them. A similarity of J = 1 denotes perfect agreement, while J = 0 would mean that the hypergraphs share no hyperedges. (When computing J, we ignore duplicate hyperedges since they are impossible to distinguish from the projection. For instance, if the board of many companies comprises the exact same directors, than we encode their association with a single hyperedge.) To obtain a baseline, we also attempt a reconstruc- tion by identifying the maximal cliques of the projected graph to hyperedges-a maximal clique reconstruction. We find that the reconstruction given by our method is good but imperfect, which is expected as the problem is under-determined. However, we also find that our method systematically outperforms the maximal clique decomposition, often by a sizable margin. In many cases, the maximal clique decomposition recovers nearly none of the common interactions, whereas our method reconstructs the interactions to a great extent.

F. Detailed case study of higher-order interactions in an empirical network
To understand why our method works well in practice, it is useful to analyze a small empirical dataset in details. For this example, we will consider the well-known football network [51]. The nodes of this network represent teams playing in Division I-A of the NCAA (now the NCAA Division I Football Bowl Subdivision), and two teams are connected if they played at least one game during the regular season of Fall 2000. The relationships between teams are viewed through the lens of pairwise interactions, but higher-order phenomena shape the system. For example, the teams of a conference all play each other during a season. Other higher-order phenomena such as geography also intervene: teams in different conferences are likely to meet during the regular season if they are in close-by states. There might also be more subtle phenomena like historical rivalries that survived conference changes. Which of these higher-order organizing principles best determine the structure is not that clear, so there is no single natural bipartite representation of the system-it is best to work with the projected network and let the data guide us.

Best model fit
In Fig. 5 we show the interactions that our method uncovers when we look for the single best higher-order description H * . We find a large number of interactions that are not pairwise: 30 of the hyperedges of H * involve more than 2 nodes.
The higher-order interactions uncovered by our method are not merely the maximal cliques of G (see Fig. 5a). We argue that interlocked maximal cliquescliques that share edges-are the reason why these descriptions differ. When two maximal cliques interlock, the hypergraph constructed directly from maximal cliques contains two overlapping hyperedges. This choice is wasteful from a compression perspective: the edges in the intersection of the two cliques are part of two hyperedges, and therefore contribute twice to the description length Σ = − log P (H). Our method instead looks for a more parsimonious description of the data. In doing so, it can identify trade-offs and, for example, represent one of the two cliques as a higher-order interaction and break down the other as a series of smaller interactions, thereby avoiding redundancies. These trade-offs culminate into much better compression: we find a hypergraph H * with a description length of 2411.3 bits, which represents a 43.2% saving over the description length of the maximal cliques hypergraph (4246.5 bits). The interactions in the optimal hypergraphs do not necessarily map to obvious suspect like sub-divisions or geographical clusters; instead, they interact in non-obvious ways and reveal, for example, that one of the sub-divisions (top left of Fig. 5b) is best described as two interlocking large hyperedges with a few interactions.

Probabilistic descriptions
Being Bayesian, our method provides complete estimation procedures, beyond maximum a posteriori estimation. For example, a quantity of particular interest is the posterior probability that a set of nodes is connected by at least one hyperedge [28], once we account for the full distribution over hypergraphs P (H|G). By computing this probability for all sets of nodes with a non-negligible connection probability, we can encode the probabilistic structure of H in a compact way [52], with a few probabilities only.
In practice, we evaluate the connection probabilities by generating samples from P (H|G) and counting the samples in which a set of interest is connected by at least one hyperedge (recall that the model defines a distribution over hypergraph with repeated hyperedges). Mathematically this is computed as where H 1 , ..., H n are n hypergraph sampled from P (H|G), and where X i1,..,i k (H) = 1 A i 1 ,..,i k ≥1 is a presence/absence variable, equal to 1 if and only if there is at least one hyperedge connecting nodes i 1 , .., i k in hypergraph H.
Applying this technique to the Football data, we find that many of the hyperedges of H * have a presence probability close to 1, even once we account for the full dis- tribution over hypergraphs. The hypergraph is not reconstructed with absolute certainty, however. Observing that probabilities P (X i1,...i k = 1|G) close to 1 or 0 both indicate confidence in the presence/absence of edge, we define a certainty threshold α and classify all hyperedges with existence probabilities in [α, 1 − α] as uncertain. With a threshold of α = 0.05, we find 16 uncertain triangles (hyperedges on 3 nodes), 70 uncertain edges, and 9 additional uncertain interactions of higher orders.
To go beyond a simple threshold analysis, we compute the entropy of the probabilitiesp := P (X i1,...i k = 1|G), defined as The entropy provides a useful transformation ofp because it grows asp moves away from the extremesp = 0, 1, with a maximum of S = 1 atp = 1/2, the point of maximal uncertainty. The distribution of entropy is shown in Fig. 6a for the Football data. The figure shows that while the majority of hyperedges are certain (i.e., their entropy is greater than S * ≈ 0.286 corresponding to α = 0.05), making the certainty criterion slightly more stringent would add many more uncertain hyperedges to the ones we already have. In Fig. 6b, we show the location of the uncertain hyperedges in H. We observe that these uncertain hyper-edges are often concentrated together. The minimality properties of P (H|G) discussed above can explain these results. Hypergraphs that have a sizable posterior probability are typically sparse and include as few hyperedges as possible. But they also need to cover the whole graph, meaning than at every edge of G needs to appear as a subset of at least one hyperedge H (due to the constraint G = G(H)). Co-located uncertain triangles and hyperedges are hence the result of competing solutions of roughly equal qualities, that cover a specific part of the hypergraph with hyperedges of different sizes.
For each empirical network in our list, we first search for the hypergraph H * that maximizes P (H|G), as we have done in our two previous examples. This search gives us a minimum description length Σ. For the sake of comparison, we also compute the description length Σ that we would obtain if we were to use the maximum clique decomposition to construct H naively. We note that Σ cannot be smaller than Σ because it is the description length of the starting point of the MCMC algorithm-at best, the algorithm cannot improve on Σ , and we then have Σ = Σ . The difference Σ − Σ gives the compression factor or, in other words, the number of bits we save by using the best hypergraph instead of a hypergraph of maximal cliques.
In Fig. 7a, we show the description lengths of the networks in our collection of datasets. We observe a broad range of outcomes. Compression of multiple orders of magnitude is possible in some cases, like with the political blogs data [54] highlighted in blue, while the best description is directly the maximal cliques in others, like with the Southern women interaction data [53] highlighted in yellow. We find that the average degree correlates with compression (Kendall's τ = 0.52); see Fig. 7b. This result is expected: the denser a network, the more likely it is that interlocking cliques are present, and therefore that a parsimonious description can be obtained by optimizing over P (H|D). The average local clustering coefficient C [2] is not correlated with compression, however (τ = 0.03); see Fig. 7c. Local clustering quantifies the density of closed triangles in the neighborhood of a node and is, as such, a proxy for the density of cliques. However, as our results show, C fails to capture the correct type of redundancy necessary for good compression with our model.
We note that clustering nonetheless predicts the absence of compression well: If C = 0, then there are no closed triangles in G, and it is impossible to compress the network with our method-there are no cliques, and Higher-order interaction in empirical networks. A few datasets are highlighted with colors: the Southern women interaction data [53] (yellow), the Football data [51] (orange), and the political blogs [54] (blue). (a) Description length of the hypergraph H whose hyperedges are the maximal cliques of the input network G, compared with the description length found with our method. (b,c) Compression, defined as the difference in description length, as a function of the average degree and clustering coefficient [2]. (d,e) Interaction size, averaged over all interactions in H * , as a function of the average degree and clustering coefficient. Detailed numerical results are reported in Supplementary Table 2. therefore no higher-order interactions in the data. The Southern Women [53] falls in this category because it is a bipartite network.
In Figs. 7d and 7e, we show the size of the higherorder interactions founds by our method, averaged over the hyperedges of H * . We again observe a wide range of outcomes. As a sanity check, we can confirm that the incompressible network has an average interaction size of 2. All hypergraphs are just networks in this case and therefore have no higher-order interactions. Other datasets yield hypergraphs with large interactions on average, involving as many as 5 nodes the airport network. The correlation between local properties and interaction size is also weak, though there are some dependencies (τ = 0.12 and τ = 0.09 for the degree and local clustering, respectively). These might be partly explained by constrained on the possible values that the average interaction size s can adopt. For instance, to have an average size s , a network must have an average degree of at least s − 1. Likewise, some level of clustering is required to obtain large interactions.
Summarizing these results, we find that some level of compression is always possible, except when the network has no clustering whatsoever. Furthermore, we find that a high average degree is related to more compression and larger higher-order interactions. Finally, we find that some minimal level of clustering is necessary for compression, but that results vary otherwise.

III. CONCLUSION
Higher-order interactions shape most relational data sets [7,26], even when they are not explicitly encoded. In this work, we have shown that it is possible to recover these interactions from data. We have argued that while the problem is ill-defined, one can introduce regularization in the form of a Bayesian generative model, and obtain a principled recovery method.
The framework we have presented is close in spirit to precursors who have used generative model to find small patterns in networks, so it is worth pointing out where it differs, both in its methodological details and philosophical underpinning. Closely related work include that of Wegner [23], who used a notion of probabilistic subgraphs covers to induce distributions over possible decomposition in motifs, and more recent works in graph machine learning that solve graph compression by decomposing the network in small building blocks [24,25]. Unlike these authors, however, we focused on higherorder interactions, so we considered decompositions in hyperedges rather than in general motif grammars. We also differ on a methodological ground: we embraced the complexity of the problem and proposed a fully Bayesian method that can account for the multiplicity of descriptions, in contrast with the greedy optimization favored in Ref. [23][24][25]. As a result of these methodological choices, our work is perhaps closest to that of Williamson and Tec [1], who also solved a similar problem from using Bayesian nonparametric techniques [1], and view a net-work as collections of overlapping cliques. Unlike these authors, however, we have formalized network data as uncorrupted; in our framework, latent higher-order interactions always show up in network data as fully connected cliques. In contrast, they think of this process as noisy, so latent higher-order interactions can translate into relatively sparsely connected sets of nodes. Their proposed methods thus bear a resemblance to community detection techniques that formalize communities as noisily measured cliques [27,[65][66][67][68][69].
The method we have proposed here is undoubtedly one of the simplest instantiations of the broader idea of uncovering higher-order interactions in empirical relational data. There are many ways in which one could expand on the method. On the modeling front, for example, it would be worthwhile to study the interplay of the projection component P (G|H) of Eq. (2) and inference: can it be defined in a way that does not turn higher-order interaction discovery into overlapping community detection? The hypergraph prior, too, will have to be expanded as the Poisson Random Hypergraphs Model (PRHM) we have used is pretty simple. Interesting models could include degree heterogeneity as part of the reconstruction [70][71][72], or community structure [73]. One could also envision a simplicial analog to these models, leading to probabilistic simplicial complex recovery [74,75]. Finally, it would be interesting to explore the connection between different forms of regularizations that make the problem well-defined.
On the technical front, it will be interesting to see whether more refined MCMC methods can lead to more robust convergence and faster mixing. Our proposed move-set is among the simplest one can propose for the problem and could be improved. Another interesting avenue of research will be to harness the known proper-ties of P (H|G) to construct efficient inference algorithms, and perhaps connect the method to algorithms in the graph theory of clique covers.
The higher-order interaction data we need to inform the development of higher-order network science [7] are often inaccessible. Our methods provide the tools needed to extract higher-order structures from much more accessible and abundant relation data. With this work, we hope to have shown that moving to principled techniques is possible, and that ad hoc reconstruction methods should be avoided, in favor of those based on informationtheoretic parsimony and statistical evidence.

DATA AVAILABILITY
The network data that support the findings of this study are available online in the Netzschleuder network repository [76].

CODE AVAILABILITY
A fast implementation of the Markov Chain Monte Carlo algorithm described in this study is freely available as part of the graph-tool Python library [77].