Growing Hypergraphs with Homophily
July 17, 2026
Abstract
There are many extant models of hypergraphs with interactions governed by attribute-based homophily between nodes, but most assume independence between edges conditional on node parameters. Relaxing this assumption, we study a mechanistic model of growing hypergraphs in which edge formation is influenced by both previous edges and binary node labels. Edges form in this model as noisy copies of previous edges, where the transmission of nodes from one edge to the next depends on multiple homophilic mechanisms between labels. These homophilic mechanisms give rise to tunable assortative structure in the hypergraph. We derive a power law for the degree distribution in this model and describe the long-term dynamics of the joint distribution of labels contained in edges. Our model defines a likelihood over a labeled hypergraph, allowing us to use standard maximum-likelihood techniques to structure algorithms. We estimate the model parameters on synthetic and real data via (stochastic) expectation maximization. These estimates give statistically-principled descriptions of the operation of homophily in empirical polyadic systems. We also demonstrate an approach to community detection via simulated annealing which, though computationally expensive, achieves competitive results on both synthetic data and certain empirical data sets known to be challenging to community detection techniques based on the edge-independence assumption. Our findings highlight the benefits of incorporating edge- and label-dependence in higher-order modeling and data analysis, and point to several directions for future work.
1 Introduction↩︎
Hypergraphs provide a natural, flexible data representation for complex systems with polyadic (multi-way) interactions; such systems form one class of higher-order networks [1], [2]. Systems naturally modeled as hypergraphs include social interaction networks (both in-person [3] and online [4]), academic co-authorship neteworks, [5], and networks of ecological relations [6]. Much recent work has emphasized the importance of faithfully representing polyadic structure rather than reducing it to dyadic interactions [1], [2] via graph projections.
Our focus in this work is the mechanism of homophily—the tendency of nodes to interact with other nodes with similar attributes—in hypergraph data. One way to model homophily in hypergraphs is via stochastic generative models, which define probability distributions over hypergraphs. In one standard framework, a set \(\mathcal{N}\) of \(n\) nodes is initialized along with a vector \(\mathbf{z}\in \{1,\ldots,k\}^n\) that assigns to each node one of \(k\) community labels. One then forms the edge set \(\mathcal{E}\) iteratively by visiting each subset \(R \subseteq \mathcal{N}\) and independently assigning \(a_R \in \mathbb{N}\) (possibly repeated) edges to \(R\) according to a distribution \(p(a_R | \mathbf{z}_R ;\boldsymbol{\theta})\), where \(\mathbf{z}_R\) is the vector of labels of the nodes in \(R\) and \(\boldsymbol{\theta}\) is a set of model parameters. Independence of sampling results in a probability distribution over hypergraphs which can be factorized into a product of edge probabilities conditioned on the node labels and model parameters: \[\begin{align} p(\mathcal{H}| \mathbf{z}; \boldsymbol{\theta}) = \prod_{R \subseteq \mathcal{N}} p(a_R | \mathbf{z}_R; \boldsymbol{\theta}) \;. \label{eq:edge-independence} \end{align}\tag{1}\] Many generative models fall into this broad scheme, including hypergraph stochastic block models [7]–[9] and some of their generalizations [10], [11].
Models of this type often underly algorithms for hypergraph community detection (also called clustering or partitioning), in which one seeks to estimate the latent node labels \(\mathbf{z}\) given an observed hypergraph \(\mathcal{H}\). The conditional independence property 1 is extremely useful for computation in this context. Factorizability of the likelihood enables many efficiencies in various inference algorithms, including approximate maximum-likelihood estimation via direct optimization [7] or expectation-maximization [8]; spectral methods [12]–[14]; and message passing [15]. Another recent optimization-based approach [16] is equivalent to a microcanonical version of the stochastic block model, which satisfies an asymptotic form of edge-independence in the large sparse limit. Latent-space models (e.g. [17]) also satisfy a form of edge independence, albeit in the context of a very different stochastic generating process.
Edge independence, however, is a strong and often unrealistic assumption for modeling hypergraph data. A parallel line of work has found that polyadic data sets often contain significant edge correlations [4], [18], [19]. These correlations are often observed as similarity between interactions. When an interaction occurs on a set \(R \subseteq \mathcal{N}\), this increases the likelihood of observing a new interaction on a set \(S \subseteq \mathcal{N}\) correlated with \(R\) in the sense that \(R\cap S\) is larger than would be expected by chance. That is, new interactions are likely to have large overlap with the interactions we have already observed. These overlaps can be shown to affect the kinds of dynamics which can unfold on the hypergraph [20], [21]. Several generative models for the growth of hypergraphs with edge-dependence have been proposed [22]–[26], including several with explicit nodel labels and corresponding homophilic mechanisms [27]–[29]; see [30] for a recent review. We note that the above are all models of hypergraph formation rather than models of dynamics on hypergraphs; recent work argues that the latter may be better expressed via dyadic graph models [31] but does not treat the former.
Of the generative models of hypergraph growth with both edge-dependence and node labels, we are unaware of any which allow for statistically-motivated inference of either the dynamical parameters or the latent node labels. In this paper, we aim to connect stochastic generative modeling of edge-dependent, node-labeled hypergraphs to inference algorithms. To this end, we introduce a model of hypergraphs with binary node labels in which edges form via a noisy edge copy process tuned by homophilic interactions between the node labels. We refer to this model as Copying Hyperedges Influenced by Label Interactions (CHILI). Our presentation of CHILI proceeds as follows. In 2, we describe CHILI and derive limiting descriptions of the joint distribution of node labels within a given edge, as well as the degree distribution of the hypergraph. We turn to optimization-based inference methods in 3, considering both parameter inference and community detection. Evaluating the likelihood of our model requires marginalizing over a set of latent variables that govern the noisy copy mechanism, which is computationally expensive. Therefore, to estimate our model’s parameters when node labels are known, we describe a stochastic expectation-maximization algorithm with favorable scaling. Community detection in this model is much more computationally challenging. As an early step in this direction,we describe a community detection algorithm for estimating node labels via simulated annealing on the CHILI likelihood. In 4, we demonstrate our algorithms on synthetic and empirical data. We show that, on synthetic data sampled from CHILI, our stochastic expectation maximization algorithm is both consistent and scalable, and that our simulated annealing algorithm for community detection, though computationally expensive, outperforms other standard methods in certain parameter regimes. On empirical data, our parameter estimates yield insight into the mechanisms governing edge formation. CHILI community detection achieves competitive recovery of observed metadata labels on some data sets known to be challenging to other recent community detection methods. We conclude with some discussion of future directions for models of node-labeled hypergraphs with edge-dependence. Our findings highlight the benefits of using explicit stochastic generative models to guide data analysis algorithms for studying homophily and assortativity in polyadic systems.
2 The CHILI Model↩︎
We propose Copying Hyperedges Influenced by Label Interactions (CHILI), a generative model of growing hypergraphs with binary node labels. CHILI is an extension of the Hyperedge Copy Model (HCM) of He et al. [24]. In the HCM, a new edge \(e\) is formed through a combination of three distinct steps: a copy step, in which a previous edge \(f\) is sampled and imperfectly copied into \(e\); an extant node addition step, in which nodes in the remainder of the node set, \(\mathcal{N}\setminus f\), are added to \(e\); and a novel node addition step, in which nodes are added to both \(e\) and \(\mathcal{N}\) for the first time. CHILI replicates these mechanisms, but adds label-aware parameters to regulate each.
Formally, a CHILI model state at time \(t\) is a hypergraph \(\mathcal{H}_{t} = \left( \mathcal{N}_{t}, \mathcal{E}_{t} \right)\), where \(\mathcal{N}_{t}\) is the set of nodes and \(\mathcal{E}_{t}\) is the set of hyperedges, as well as a binary label vector \(\mathbf{z}_{t} \in \{0, 1\}^{|\mathcal{N}_{t}|}\) which assigns to each node a label of either 0 or 1; \(z_u\) denotes the label of node \(u\). We let \(\bar{z} = 1 - z\) denote the opposite label of \(z\). At timestep \(t+1\), we add edge \(e\) to \(\mathcal{E}_{t}\) to obtain \(\mathcal{E}_{t+1}\), where \(e\) is formed according to the following three-step process:
Edge-Copying: An edge \(f\), called the seed edge, is sampled uniformly at random from \(\mathcal{E}_{t}\), and a node \(u\), called the focal node, is sampled uniformly at random from \(f\).
We refer to \(f\) as the seed edge and \(u\) as the focal node. The node \(u\) is included in \(e\). Then, each node \(v\) in \(f\setminus \{u\}\) is included in \(e\) with probability \(\rho{}_+\) if \(z_{v} = z_{u}\) and with probability \(\rho{}_-\) if \(z_{v} \neq z_{u}\).
Extant Node Addition: We sample \(X_{z_u} \sim \mathrm{Poisson}\left( \gamma{}_+ \right)\) nodes from \(\mathcal{N}_{t} \setminus f\) with label \(z_u\) and add them to \(e\). We similarly sample \(X_{\bar{z}_u} \sim \mathrm{Poisson}\left( \gamma{}_- \right)\) nodes from \(\mathcal{N}_{t} \setminus f\) with label \(\bar{z}_u\) and add them to \(e\). In the event that there are not enough nodes with the required label in \(\mathcal{N}_{t} \setminus f\), we add all available nodes with that label to \(e\).
Novel Node Addition: We sample \(Y_{z_u} \sim \mathrm{Poisson}\left( \eta{}_+ \right)\) and add \(Y_{z_u}\) nodes with label \(z_u\) to both \(e\) and \(\mathcal{N}_{t+1}\). We similarly sample \(Y_{\bar{z}_u} \sim \mathrm{Poisson}\left( \eta{}_- \right)\) and add \(Y_{\bar{z}_u}\) nodes with label \(\bar{z}_u\) to both \(e\) and \(\mathcal{N}_{t+1}\). The associated labels are also added to \(\mathbf{z}_{t}\) to form \(\mathbf{z}_{t+1}\).
We let \(\boldsymbol{\theta}= (\rho{}_+, \rho{}_-, \gamma{}_+, \gamma{}_-, \eta{}_+, \eta{}_-)\) denote the complete vector of model parameters. Formally, the probability of realizing a given hypergraph \(\mathcal{H}\) with label vector \(\mathbf{z}\) under CHILI at time \(t\) can be described as a time-dependent probability distribution \(p_{t}(\mathcal{H}, \mathbf{z}; \boldsymbol{\theta})\) over the space of hypergraphs with labeled nodes; we omit the time index when otherwise clear from context. Our implementation of CHILI the hypergraph data structure supplied by the XGI package for polyadic data analysis in Python [32].
1 illustrates the role of the model parameters in tuning the community structure of realized hypergraphs. We measure community structure by computing Newman-Girvan modularity [33] on the weighted clique-projection of the hypergraph. We also show several sample hypergraphs at different parameter combinations for visual illustration. For fixed \(\gamma_-\), \(\eta{}_+\), and \(\eta{}_-\), increasing the difference \(\rho{}_+ - \rho{}_-\) leads to higher modularity, as does increasing \(\gamma{}_+\). Hypergraphs sampled from our model using high values of \(\rho{}_+\) and \(\gamma{}_+\) tend to have label-homogeneous cores of densely overlapping hyperedges, which are weakly connected to one another by edges with more heterogeneous label compositions. Hypergraphs with lower values of \(\rho{}_+\) and \(\gamma{}_+\) tend to have less-structured distributions of labels over hyperedges. The CHILI model can also be used to generate hypergraphs under heterophilic conditions.
2.1 Asymptotic Properties↩︎
2.1.1 Joint Distribution of Labels in Edges↩︎
Many asymptotic properties of our CHILI model can be described in terms of the joint distribution \(q(k_0, k_1)\) of node labels in an edge; here \(q(k_0, k_1)\) gives the probability that an edge selected uniformly at random from \(\mathcal{E}\) contains exactly \(k_0\) nodes of label 0 and \(k_1\) nodes of label 1. Under the \(t\rightarrow \infty\) limit this joint distribution can be obtained as the leading eigenvector of a linear map. Let \(q_{t}(k_0, k_1)\) be the joint distribution of label counts in a uniformly random edge after \(t\) timesteps, and let \(q(k_0, k_1) = q_{\infty} \left( k_0, k_1 \right)\) be the stationary distribution. Then, the probability \(r(k_0', k_1')\) of realizing an edge with \(k_0'\) nodes of label 0 and \(k_1'\) nodes of label 1 in timestep \(t+1\) is given by \[\begin{align} r(k_0', k_1') = \sum_{k_0, k_1} T(k_0', k_1' | k_0, k_1)q_{t}(k_0, k_1) \;, \end{align}\] where \(T(k_0', k_1' | k_0, k_1)\) gives the probability that the edge \(e\) realized in step \(t+1\) has \(k_0'\) nodes of label 0 and \(k_1'\) nodes of label 1 given that the edge \(f\) sampled to form \(e\) has \(k_0\) nodes of label 0 and \(k_1\) nodes of label 1. At stationarity, we must have \(r(k_0', k_1') = q_{t}(k_0', k_1') = q\left( k_0', k_1' \right)\), so the stationarity condition can be obtained by solving the linear system \[\begin{align} q\left( k_0', k_1' \right) = \sum_{k_0, k_1} T(k_0', k_1' | k_0, k_1)q\left( k_0, k_1 \right) \;. \end{align}\] The stationary joint distribution \(q\left( k_0, k_1 \right)\) is therefore given by the leading eigenvector of the doubly-indexed matrix \(T\) with entries \(T(k_0', k_1' | k_0, k_1)\). The entries of this matrix follow from the parametric probability descriptions used to define the CHILI model. Conditional on choosing an edge \(f\) with \(k_0\) nodes of label 0 and \(k_1\) nodes of label 1, a node of label \(z\) is selected as the focal node \(u\) with probability \(\frac{k_z}{k_0 + k_1}\). Then, the number of nodes of labels \(z\) and \(\bar{z}\) in \(e\) have distributions \[\begin{align} k_z' &\sim 1 + \text{Binom}\left( k_z - 1, \rho{}_+ \right) + \text{Poisson}(\gamma{}_+ + \eta{}_+) \;, \\ k_{\bar{z}}' &\sim \text{Binom}\left( k_{\bar{z}}, \rho{}_- \right) + \text{Poisson}(\gamma{}_- + \eta{}_-) \;. \end{align}\] The probability of obtaining \(k_0'\) nodes of label 0, conditioning on \(k_0\) and \(k_1\), is
\[\begin{align} t(k_0' | k_0, k_1) =& \frac{k_0}{k_0 + k_1} \sum_{i=0}^{k_0' - 1} \binom{k_0 - 1}{i} \rho{}_+^i (1 - \rho{}_+)^{k_0 - 1 - i} e^{-(\gamma{}_+ + \eta{}_+)} \frac{(\gamma{}_+ + \eta{}_+)^{k_0' - 1 - i}}{(k_0' - 1 - i)!} + \\ & \frac{k_1}{k_0 + k_1} \sum_{i=0}^{k_0'} \binom{k_0}{i} \rho{}_-^i (1 - \rho{}_-)^{k_0 - i} e^{-(\gamma{}_- + \eta{}_-)} \frac{(\gamma{}_- + \eta{}_-)^{k_0' - i}}{(k_0' - i)!}\;. \end{align}\]
A parallel expression describes \(t(k_1' | k_0, k_1)\), and the transition matrix entries are given by the product \[\begin{align} T(k_0', k_1' | k_0, k_1) = t(k_0' | k_0, k_1) t(k_1' | k_0, k_1) \;. \end{align}\] In 2, we study the dynamics of our CHILI model as described by the matrix \(T\). The spectral gap \(1 - \lambda_2\) of \(T\) is a measure of the rate of convergence to the stationary distribution \(q\left( k_0, k_1 \right)\), where \(\lambda_2\) is the second-largest eigenvalue of \(T\). For fixed values of the parameters \(\gamma{}_+\), \(\gamma{}_-\), \(\eta{}_+\) and \(\eta{}_-\), the spectral gap \(1 - \lambda_2\) of \(T\) is highest when \(\rho_-\) is large and \(\rho_+\) is small. High values of the spectral gap correspond to fast convergence to the stationary distribution \(q\left( k_0, k_1 \right)\) given by the Perron eigenvector \(\mathbf{v}_1\) of \(T\) with eigenvalue \(1\) (panel (b)). We instead illustrate an alternative regime of relatively slow convergence in which \(\rho_+\) is large and \(\rho_-\) is small. In this regime, the simulation may maintain transient dynamics in the direction of the second eigenvector \(\mathbf{v}_2\) of \(T\) with eigenvalue \(\lambda_2 < 1\) (panel (c)) for relatively long periods of time before converging to the stationary distribution. We illustrate this behavior in panel (d), where the system converges relatively slowly to the leading eigenvector \(\mathbf{v}_1\) and away from the transient component \(\mathbf{v}_2\). The resulting final simulation state matches the stationary distribution \(q\left( k_0, k_1 \right)\) relatively closely (panel (e)).
2 shows that, when \(\rho{}_+\) is sufficiently large, the stationary distribution \(q\left( k_0, k_1 \right)\) is concentrated on hypergraphs in which typical edges have large majorities of nodes with the same label (panel b), the signature of higher-order structural assortativity. This distribution is symmetric with respect to the two node labels, however, indicating that although local edges may have dominant majority labels, the long-run dynamics converge to equal global label populations. The second eigenvector \(v_2\) (panel c), however, is antisymmetric, describing a transient in which a single label group tends to dominate most edges. When \(\rho{}_+\) is sufficiently large, the spectral gap \(1 - \lambda_2\) is small (panel a), indicating that this transient decays slowly and that the system may be dominated by edges with a single majority label for relatively long periods of time before equilibrating to the leading eigenvector (panels (d-e)). This behavior suggests that, although the CHILI model ultimately converges to a stationary state with symmetry between the two node labels, it can also serve through its transients as a reasonable generative model of systems with label asymmetry.
2.1.2 Degree Distribution↩︎
One model feature we can obtain from the stationary joint distribution \(q\left( k_0, k_1 \right)\) is the asymptotic degree distribution generated from CHILI. This asymptotic distribution follows a power law with an exponent which can be computed in closed form from the model parameters and the stationary joint distribution \(q\left( k_0, k_1 \right)\). Our derivation is similar to that of He et al. [24], which in turn is based on a derivation by Mitzenmacher [34]. The argument requires a mean-field assumption: when a node \(u\) is selected from edge \(f\) in the edge-sampling step, the degree \(d\) of \(u\) is independent of the ratio \(\frac{k_0}{k}\) of label-0 nodes in the edge \(f\) sampled to form the new edge \(e\). Under this assumption, it is possible to show that the degree distribution follows a power law with exponent \[\begin{align} \zeta &= 1 + \frac{\left\langle k_0 \right\rangle + \left\langle k_1 \right\rangle}{1 + \rho{}_+(\mu_{00} + \mu_{11} - 1) + 2\rho{}_- \mu_{01}} \label{eq:power-law-exponent}\\ \mu_{ij} &= \left\langle \frac{k_i k_j}{k} \right\rangle \;. \nonumber \end{align}\tag{2}\] In the expression for \(\mu_{ij}\) (and for \(\left\langle k_i \right\rangle\)) the expectation is taken with respect to the stationary joint distribution \(q\left( k_0, k_1 \right)\). A derivation of this result and a computational illustration are provided in Appendix 6 and 8.
3 Model Inference↩︎
In this section, we consider the problems of parameter inference (learning \(\boldsymbol{\theta}\) given \(\mathcal{H}\) and \(\mathbf{z}\)) and community detection (learning \(\mathbf{z}\) given \(\mathcal{H}\) and \(\boldsymbol{\theta}\)).
3.1 Parameter Estimation via Stochastic EM↩︎
We first consider the problem of learning the parameter vector \(\boldsymbol{\theta}= (\rho{}_+, \rho{}_-, \gamma{}_+, \gamma{}_-, \eta{}_+, \eta{}_-)\) given an observed hypergraph \(\mathcal{H}\) with sequential edge indices and known node labels \(\mathbf{z}\). We aim to do this via maximum-likelihood estimation. Computing the complete data likelihood requires marginalizing over the model latent variables, which in this case are the edge \(f\) and node \(u\) selected at each timestep of the growth model to form the new edge \(e\). For \(e\in \mathcal{E}\), let \(t_e\) denote the time at which \(e\) was added to the hypergraph. We say that \(f\prec e\) if \(t_{f} < t_{e}\). Then, the complete data likelihood can be written
\[\begin{align} p(\mathcal{H}, \mathbf{z}; \boldsymbol{\theta}) = \prod_{e\in \mathcal{E}} \sum_{f\prec e} \sum_{u\in f} p(e, \mathbf{z}'_e | f, u, \mathbf{z}_{t_e-1}; \boldsymbol{\theta}) p(f,u| \mathcal{E}_{t_e-1}, \mathbf{z}_{t_e-1}; \boldsymbol{\theta}) \;, \label{eq:likelihood} \end{align}\tag{3}\]
where \(\mathbf{z}_{t_e-1}\) is the vector of node labels immediately before \(e\) was added, and \(\mathbf{z}'_e\) is the vector of node labels added with edge \(e\) at time \(t_e\). The term \(p(f,u| \mathcal{E}_{t_e-1}, \mathbf{z}_{t_e-1}; \boldsymbol{\theta})\) gives the likelihood of selecting edge \(f\) and node \(u\) at time \(t - 1\); its value is \(\frac{1}{\left| f \right|m_{t-1}}\) if \(f\in \mathcal{E}_{t_e-1}\) and \(0\) otherwise, where \(m_{t-1} = \left| \mathcal{E}_{t_e-1} \right|\) is the number of edges at time \(t - 1\). The term \(p(e, \mathbf{z}'_e | f, u, \mathbf{z}_{t_e-1}; \boldsymbol{\theta})\) gives the likelihood of realizing edge \(e\) and the associated new node labels \(\mathbf{z}'_e\) given that edge \(f\) and node \(u\) were selected at time \(t - 1\). A helpful note for computation is that \(p(e, \mathbf{z}'_e | f, u, \mathbf{z}_{t_e-1}; \hat{\boldsymbol{\theta}}^{[\ell]}) = 0\) if either \(u\notin e\cap f\) or \(f\not\prec e\).
The requirement in 3 to sum over the latent variables at each timestep implies a computational cost to evaluate \(p(\mathcal{H}, \mathbf{z}; \boldsymbol{\theta})\) (or its derivatives) which is quadratic in the number of edges \(m\), making direct optimization of the likelihood expensive for data sets of even modest size. Batch expectation-maximization [35] is an alternative approach which maximizes the likelihood without explicitly computing 3 , but even this approach requires forming a belief distribution over the values of all latent variables at each timestep, which again incurs quadratic cost per iteration.
We therefore use a stochastic expectation-maximization (SEM) algorithm [36] to estimate the parameters \(\boldsymbol{\theta}\). Each iteration of SEM consists of two steps: the stochatic E-step, in which noisy exponentially weighted moving averages of the sufficient statistics are calculated using the current parameter estimates, and the stochastic M-step, in which the parameter estimates are updated using the expected sufficient statistics.
To start the SEM algorithm, we initialize a vector \(\mathbf{s}^{[0]}\) of sufficient statistics to arbitrary values; in the CHILI model there are two sufficient statistics per copy parameter and one sufficient statistic per each of the four addition parameters, so \(\mathbf{s}\) has length \(8\). From \(\mathbf{s}^{[0]}\), we compute an initial estimate of the parameter vector \(\hat{\boldsymbol{\theta}}^{[0]} = g(\mathbf{s}^{[0]})\) via a function \(g\) which maps sufficient statistics to parameter estimates; details on both the sufficient statistics \(\mathbf{s}\) and the mapping \(g\) are provided in 7. We then update a moving average of \(\mathbf{s}\) as follows. In each algorithmic step \(\ell\), we pick an edge \(e\) uniformly at random from \(\mathcal{H}\). Then, given current estimates \(\hat{\boldsymbol{\theta}}^{[\ell]}\), we form a belief distribution over the latent identities of \(f\in \mathcal{E}\) and \(u\in f\) which could have been selected in the edge- step to form \(e\). Using Bayes’ rule, this probability is
\[\begin{align} p(f,u| e, \mathbf{z}'_e; \hat{\boldsymbol{\theta}}^{[\ell]}) = \frac{p(e, \mathbf{z}'_e | f, u, \mathbf{z}_{t_e-1}; \hat{\boldsymbol{\theta}}^{[\ell]}) \left| f \right|^{-1}}{\sum_{f' \prec e} \sum_{u' \in f'} p(e, \mathbf{z}'_e | f', u', \mathbf{z}_{t_e-1}; \hat{\boldsymbol{\theta}}^{[\ell]}) \left| f' \right|^{-1} } \;. \label{eq:bayes} \end{align}\tag{4}\]
We then form the vector of expected sufficient statistics \[\begin{align} s' = \sum_{f\prec e} \sum_{u\in f\cap e} p(f, u| e, \mathbf{z}'_e; \hat{\boldsymbol{\theta}}^{[\ell]}) \boldsymbol{\psi}(e, f, u, \mathbf{z}) \;, \label{eq:sufficient-stats-update} \end{align}\tag{5}\] where \(\boldsymbol{\psi}(e, f, u, \mathbf{z})\) is the vector of sufficient statistics for the parameters \(\boldsymbol{\theta}\) given \(e\), \(f\), \(u\), and \(\mathbf{z}\). The entries of \(\boldsymbol{\psi}(e, f, u, \mathbf{z})\) can be computed as functions of the node labels in \(e\) and \(f\); we supply explicit formulae in 7. We then perform the update \[\begin{align} \mathbf{s}^{[\ell+1]} = (1 - \alpha^{[\ell]} )\mathbf{s}^{[\ell]} + \alpha^{[\ell]} \mathbf{s}' \;, \end{align}\] where \(\alpha^{[\ell]}\) is the learning rate at iteration \(\ell\); we use an exponentially decaying learning rate in our experiments. Our new estimate of the parameter vector \(\boldsymbol{\theta}\) is then given by a function \(\hat{\boldsymbol{\theta}}^{[\ell+1]} = g(\mathbf{s}^{[\ell+1]})\) of the sufficient statistics, which we also supply in 7. 3 summarizes this process. Our version of the algorithm uses an online computation of \(\mathbf{s}'\), which allows us to perform the double-sum appearing in 5 4 only once, and without explicitly storing the resulting belief distribution at any point.
We perform SEM by repeating 3 until a stopping criterion is met. For our criterion, we define the function \[\begin{align} h(\hat{\boldsymbol{\theta}}^{[\ell]}, \hat{\boldsymbol{\theta}}^{[\ell']}) = \min_{i} \left\{\left| \hat{\theta}_i^{[\ell]} - \hat{\theta}_i^{[\ell']} \right| - \epsilon_\mathrm{rel} - \epsilon_\mathrm{abs} \hat{\theta}_i^{[\ell']}\right\} \;. \end{align}\] The function \(h\) here can be viewed as combining absolute and relative error tolerances, enabling the algorithm to declare convergence even in cases when the parameter estimates are near zero. In our implementation, we declare convergence at timestep \(\ell\) when \(h(\hat{\boldsymbol{\theta}}^{[\ell]}, \hat{\boldsymbol{\theta}}^{[\ell']}) \geq 0\) for all \(\ell'\) in the most recent 800 iterations of the algorithm, with tolerances \(\epsilon_\mathrm{rel}\) and \(\epsilon_\mathrm{abs}\).
3.2 Community Detection via Simulated Annealing↩︎
We also consider the community detection problem of learning the node labels \(\mathbf{z}\) given an observed hypergraph \(\mathcal{H}\) and a fixed parameter vector \(\boldsymbol{\theta}\).
Since 3 defines a likelihood function over both the parameter vector \(\boldsymbol{\theta}\) and the node labels \(\mathbf{z}\), we can attempt to also learn \(\mathbf{z}\) via maximum-likelihood estimation. Unlike in stochastic block model variants in which there are several heuristics with favorable computational properties [7], [8], [11], [12], [14], [16], [37], the edge dependence structure encoded by CHILI makes it difficult to apply standard efficient heuristics. The development of alternative efficient heuristics is an important direction of future work. In our present approach, we therefore attempt to learn \(\mathbf{z}\) via simulated annealing [38], a general-purpose optimization algorithm which is typically used when other heuristics are not available. Our primary aim is to demonstrate that likelihood maximization in this model can detect label structure in at least some circumstances where other community detection methods perform poorly, thereby motivating the development of more efficient optimization heuristics in future work.
A direct evaluation of 3 requires \(O(m^2)\) terms, making it expensive to evaluate even for data sets of modest size. An important computational limitation of simulated annealing in this context is that the label \(z_i\) of node \(i\) appears in all terms of the likelihood 3 corresponding to edges after which node \(i\) arrives. This includes edges in which node \(i\) does not participate, implying that we must perform \(O(m^2)\) computations even to evaluate the change in likelihood under a proposed change to \(z_i\). This contrasts to, for example, modularity [33], where the computational cost of moving a node between communities is constant under appropriate use of data structures.
To mitigate the computational burden, we use an intersection-based approximation to the likelihood of forming edge \(e\) in which we only sum over candidate copied edges \(f\) such that \(\left| f\cap e \right|\) is large. The intuition of this approximation, which we justify in detail in 9, is that \(p(e, \mathbf{z}_e'|f, u, \mathbf{z}; \boldsymbol{\theta})\in O(n^{-x})\), where \(n\) is the number of nodes in the hypergraph and \(x = \left| e \right| - \left| e\cap f \right| - \left| e\setminus \mathcal{N} \right|\) is the number of extant nodes required to form \(e\) given \(f\). The essence of the argument is that forming \(x\) extant nodes requires \(x\) independent samples from the sets \(\mathcal{N}_0\setminus f\) and \(\mathcal{N}_1\setminus f\), each of which has probability \(O(n^{-1})\) of being selected. We therefore design our algorithm to maximize the objective function \(\tilde{\mathcal{L}}(\mathbf{z}) = \log \tilde{p}(\mathcal{H}, \mathbf{z}; \boldsymbol{\theta})\) where \(\tilde{p}(\mathcal{H}, \mathbf{z}; \boldsymbol{\theta})\) is the likelihood of \(\mathcal{H}\) and \(\mathbf{z}\) under the intersection-based approximation. A single epoch consists of \(n\) individual update steps. In an update step in epoch \(\ell\), we propose a change to the labels \(\mathbf{z}\) by generating a node index \(j\) uniformly at random and evaluating the change in likelihood of “flipping” the binary node label (\(z_j \gets \bar{z}_j\)). If the proposal increases \(\tilde{\mathcal{L}}(\mathbf{z})\), it is accepted. If the proposal decreases \(\tilde{\mathcal{L}}(\mathbf{z})\), however, we may also accept it with an acceptance probability \(a^{[\ell]}\) that decreases with the epoch \(\ell\) and considers the magnitude of the decrease in \(\tilde{\mathcal{L}}(\mathbf{z})\). It is possible to store computational components required to evaluate \(\tilde{\mathcal{L}}(\mathbf{z})\) in a data structure that allows us to compute the change in likelihood from the prior state without recomputing the entire likelihood. Although this does not result in an asymptotic improvement, experimentally we observe a roughly 30-fold reduction in compute time on our empirical data sets at the expense of increasing memory requirements.
We use a custom acceptance schedule \(a^{[\ell]}\) for our simulated annealing implementation. Our custom approach is motivated by the observation that, in our application, the change in the objective \(\tilde{\mathcal{L}}\) of a proposal can vary significantly depending on the size of the hypergraph under study. Thus, in our early experiments, traditional annealing schedules such as harmonic or exponential decay struggled to still occasionally accept proposals that decrease \(\tilde{\mathcal{L}}\). Our custom implementation instead takes an adaptive approach based on the distribution of proposals previously observed. In the first epoch, we accept all proposals, while storing the changes to the objective \(\tilde{\mathcal{L}}\) produced by those proposals. At the end of the epoch, we compute the standard deviation \(\sigma(\Delta \tilde{\mathcal{L}})\) of the changes in \(\tilde{\mathcal{L}}\) observed during the epoch. For the remaining epochs, we accept proposals that decrease \(\tilde{\mathcal{L}}\) by an amount \(\Delta \tilde{\mathcal{L}}\) stochastically by evaluating a Gaussian density with mean 0 and a standard deviation proportional to \(\sigma(\Delta \tilde{\mathcal{L}})\) which decays with epoch \(\ell\). A full specification of the algorithm is provided in 10.
4 Results↩︎
We now test the ability of our proposed CHILI model to recover parameters (via SEM) and node labels (via simulated annealing) on both synthetic and empirical data sets.
4.1 Parameter Estimation on Synthetic Data↩︎
4 illustrates the behavior of SEM in recovering parameters from synthetic hypergraphs generated by the CHILI model. Each synthetic hypergraph is generated by simulating CHILI with a known parameter vector \(\boldsymbol{\theta}\). We say that an instance of the CHILI model possesses homophily in the edge copy step of the growth process if \(\rho_+ > \rho_-\), and heterophily if \(\rho_+ < \rho_-\). We refer to homophily and heterophily in the extant and novel node addition steps in parallel fashion. We perform SEM on a set of hypergraphs that were generated with varying strengths of homophily across each of the three growth mechanisms: the edge copy step, the extant node addition step, and the novel node addition step. We also vary the expected number of nodes added to \(e\) in each of these steps. To measure the success of a parameter recovery, we compute the Kullback-Leibler (KL) divergence between the true and estimated distributions described by each of the six parameters, and add the results. Using standard results, our measure of error is given by \[\begin{align} \delta_{\mathrm{KL}} &= \phantom{+}\sum_{* \in \{+, -\}} c_*\left( \rho_{*} \log \frac{\rho_{*}}{\hat{\rho}_{*}} + (1 - \rho_{*}) \log \frac{1 - \rho_{*}}{1 - \hat{\rho}_{*}} \right) \nonumber\\ &\phantom{=} + \sum_{* \in \{+, -\}} \left( \gamma_{*} \log \frac{\gamma_{*}}{\hat{\gamma}_{*}} - \gamma_* + \hat{\gamma}_{*} \right) \label{eq:kl-divergence} \\ &\phantom{=} + \sum_{* \in \{+, -\}} \left( \eta_{*} \log \frac{\eta_{*}}{\hat{\eta}_{*}} - \eta_* + \hat{\eta}_{*} \right)\;, \nonumber \end{align}\tag{6}\] where \(c\) is a factor describing, in expectation, how many nodes of each label are selected in the copy step. For this analysis, we set \(c_* = 1\) for both labels.
As seen in 4, the quality of the learned parameter estimates depends on the seeded parameter vector. We observe especially low estimation errors when the seeded parameters reflect high copy rates \(\rho_{*}\) and strong homophily as reflected in the differences \(\gamma_+ - \gamma_-\) and \(\eta_+ - \eta_-\); a typical SEM run for these parameter settings is shown in the top row on the right of 4. That said, even in the least successful cases for estimation with high KL divergence, the learned estimates are still visually close to the true values, illustrated in the bottom row on the right of 4.
4.2 Parameter Estimation on Empirical Data↩︎
We also used our method to infer parameters for six empirical hypergraph data sets described in 1. Although many labeled hypergraph data sets are now available for community detection, relatively few have natively binary
node labels. The senate-bills and house-bills data sets [39], [40] are hypergraphs of cosponsorships of bills in the U.S. Senate and House of Representatives, respectively, with node labels corresponding to the political party of each legislator; we
use the version of these data sets used in [7]. The high-school and
primary-school data sets [3], [41], [42] are hypergraphs of social interactions between students in a high school and primary school, respectively, with node labels corresponding to the gender of
each student; these data sets were collected through the SocioPatterns project. In these data sets, we consider a sequence of interactions to correspond to a single hyperedge if those interactions are identical as sets and temporally contiguous. The
coauthorship data set [43] contains publishing records for 81 computer science conferences from 2000 to 2015, with an inferred
binary gender label for each author. Finally, we considered the enron-emails data set, which contains the email records for 145 users associated with the Enron Corporation [44].
| \(n\) | \(m\) | \(\langle d \rangle\) | \(\langle k \rangle\) | Nodes are | Edges are | Labels represent | |
|---|---|---|---|---|---|---|---|
| senate-bills [7], [39], [40] | 294 | 29,157 | 789.6 | 8.0 | Senators | Cosponsored bills | Political party |
| house-bills [7], [39], [40] | 1,494 | 60,987 | 835.8 | 20.5 | Congresspeople | Cosponsored bills | Political party |
| high-school [41], [42] | 320 | 60,544 | 484.6 | 2.6 | Students | Social interactions | Gender |
| primary-school [3] | 227 | 44,626 | 590.7 | 3.0 | Students | Social interactions | Gender |
| coauthorship [43] | 105,256 | 142,834 | 3.4 | 2.5 | Authors | Papers | Gender |
| enron-emails [44] | 14,901 | 83,372 | 18.5 | 3.3 | Email addresses | Email threads | Fringe vs. core |
5 shows our parameter estimates for each of the empirical data sets under the CHILI model. In each panel, we plot the same-label and opposite-label parameter estimates on opposing axes. Some data sets, like the
gendered primary-school contact network, appear on or near the line of equality in all three panels, suggesting that homophily with respect to the given binary labels is not a major driver of network evolution. The coauthorship
data, meanwhile, falls on the line of equality for the edge copy and extant node addition mechanisms, but there is a difference in the novel node addition parameters: \(\hat{\eta}_+ > \hat{\eta}_-\). Our model presents a
simplified picture of the process of collaboration formation as being initiated by a focal individual scientist, who then recruits from a previous collaboration, from other published scientists in the field, and from an unobserved pool of novices who have
never published before. Our findings suggest that of these mechanisms, only the latter (the recruitment of novices) displays substantial gender homophily, with novices who share the gender of the focal scientist more likely to be recruited to the
collaboration. This effect may be partially due to gender homophily between PhD students and their advisors, a known phenomenon in several scientific areas [45].
The senate-bills data set appears to have weak homophily in the estimates of the edge copy and extant node addition parameters. Because new senators are elected far less frequently than bills are cosponsored, the novel node addition
parameters are both nearly 0. The analogous house-bills data set has edge copy parameter estimates which are relatively similar to those of senate-bills, but which suggest a small tendency towards heterophily rather than homophily
in edge-copying. The slight heterophily in the copy step of house-bills may be related to reciprocity in bipartisan cosponsorship, a phenomenon that has been observed in U.S. legislative bodies [46], [47]. We illustrate the connection between reciprocity and
edge copy parameters through the following example. Consider some bill \(e\), where \(u\) is a cosponsor of both bill \(e\) and bill \(f\). If \(u\) is a Democrat, for instance, and the primary sponsor of \(f\) is a Republican, then \(u\) was a “bipartisan”
cosponsor of \(f\). The reciprocity of bipartisan cosponsorship suggests that the Republican cosponsors of \(f\) may also cosponsor \(e\) to return the
favor. However, the other Democratic cosponsors of \(f\) are not indebted in the same way to \(u\), possibly making them less likely to cosponsor \(e.\) The
mechanisms of and motivations for bipartisan reciprocity vary between the House and the Senate [47], which may in part explain why we
observe heterophily in edge copying for house-bills but not senate-bills. Further investigation of the relationship between reciprocity and edge copying in legislative cosponsorship is an interesting direction for future work.
The novel node estimates for both the congressional datasets are very small. The largest difference in our estimates for these two systems is in the extant node addition parameter estimates; the house-bills data set has estimates \(\hat{\gamma}_+\) and \(\hat{\gamma}_-\) which are each roughly three times larger than their values in the senate-bills data set. These estimates suggest the operation of similar
levels of party homophily in the cosponsorship of bills via extant node addition, while reflecting the higher mean edge size of house-bills, where a bill is co-sponsored by a mean of over 20 congresspersons, as opposed to the mean of 8 in
senate-bills.
The enron-emails hypergraph falls in the heterophilic region of the edge copy panel. This is a feature of the sampling strategy. The Enron data set contains the email data from 145 users associated with the Enron Corporation. These users
are the “core” nodes of the hypergraph. Any users who were included in enron-emails with core users, but whose complete email data was not sampled, are referred to as “fringe” nodes. Each hyperedge represents an email between a set of nodes,
or users. Therefore, each hyperedge must contain at least one core node. Since 73% of edges in the hypergraph have size 2, and less than one percent of the nodes in the hypergraph are core nodes, most edges (67%) join a core node to a fringe node,
generating disassortativity with respect to these labels. The high magnitude of the copy parameters, particularly \(\rho_-\), is a result of the nature of email — rarely does an email interaction between a group of users
occur only once. Often, an email is precipitated by an existing conversation between that set of users, making the edge copy mechanism well suited for modeling email network growth.
4.3 Community Detection via Simulated Annealing↩︎
We evaluate our simulated annealing algorithm on a range of small synthetic data sets under two conditions, shown in 6. Under the Nishimori condition [48], we assume that we have access to the exact parameter vector used to generate the synthetic data. In this condition, we therefore generate a synthetic hypergraph according to a parameter vector \(\boldsymbol{\theta}\) and then attempt to learn \(\mathbf{z}\) via simulated annealing using the same parameter vector \(\boldsymbol{\theta}\). Under the “fixed” condition, which simulates a more realistic data analysis scenario, we always use the same parameter vector \(\boldsymbol{\theta}'\) to learn \(\mathbf{z}\) via simulated annealing regardless of the true parameter vector \(\boldsymbol{\theta}\) used to generate the synthetic data. In this condition, we use a default value of \(\boldsymbol{\theta}' = (0.9, 0.1, 0.2, 0.1, 2.0, 0.3)\). We compare these methods against each other and against a baseline of Laplacian spectral clustering [49] applied to the weighted clique projection of the hypergraph. We assess the effectiveness of these techniques and our method with a measure of clustering effectiveness, the adjusted Rand index (ARI) [50]. In each panel of 6, we show the performance of each of these models as a single entry of \(\boldsymbol{\theta}\) is varied while the others are held fixed. Unsurprisingly, simulated annealing under the Nishimori condition outperforms simulated annealing under the fixed condition. For strong homophily (larger values of the varied parameter), simulated annealing under the Nishimori condition can nevertheless be outperformed by Laplacian spectral clustering. Notably, however, there are regions of parameter space, especially those with weak homophily (smaller values of the varied parameter), in which both simulated annealing under the fixed condition and under the Nishimori condition outperform Laplacian spectral clustering. This result is in some sense unsurprising, since we perform our experiment on synthetic samples drawn directly from the CHILI model. Nevertheless, these experiments suggest that there exist classes of hypergraphs in which community structure can be more reliably detected via likelihood maximization under CHILI than by at least this standard graph method.
Indeed, we observe relatively favorable performance of the CHILI model on empirical data as well. We evaluate our simulated annealing algorithm on a subset of our empirical data sets: senate-bills, primary-school, and
high-school; the poor scaling of simulated annealing with the number of nodes \(n\) made it computationally prohibitive to run on the larger house-bills, enron-emails, and
coauthorship data sets. We use the class membership of the students in primary-school and high-school as the ground truth labels in this section instead of gender because community structure has been previously
detected with these labels [7]. For our experiments we compare the binary labels
generated by our simulated annealing algorithm against the multinary ground-truth class labels. Since in the context of community detection the labels \(\mathbf{z}\) (and therefore parameter estimates \(\hat{\boldsymbol{\theta}}\)) are not available, we instead use a single, fixed parameter vector for each of these three experiments: \(\hat{\boldsymbol{\theta}} = (0.90, 0.10, 1.00, 0.25, 0.001,
0.001)\). The zero values for the novel node addition parameters reflect a choice of algorithmic convenience. All of these data sets are appropriately modeled with very small values of \(\eta_+\) and \(\eta_-\) under our model due to their large number of edges compared to the small number of nodes. To simplify computation, we instead simply assume that all nodes exist in the system at all times, even if they have not yet
participated in an edge.
We compare CHILI simulated annealing against several competitors. The all-or-nothing hypergraph modularity of [7] does not allow the user to specify a desired number of clusters, but because it returns a binary partition on senate-bills, we can compare directly on this data set. Belief-propagation nonbacktracking
spectral clustering [12] allows the user to specify the number of clusters, and the
authors report results for binary clustering. We therefore compare our method against the reported results from these two prior works.3 Notably, both hypergraph modularity
maximization and nonbacktracking spectral clustering are based on generative models which assume conditional edge independence, the assumption relaxed by CHILI. We also include as comparisons two dyadic methods, implemented on the weighted clique
projection of the hypergraph, that allow the user to specify the number of clusters: Laplacian spectral clustering [49] and greedy
modularity maximization [51], both of which also have connections to generative models that assume conditional (dyadic) edge
independence [37], [52].
![Figure 7: Community detection results on three selected empirical data sets. Box plots describe the distribution of the ARI of CHILI simulated annealing (SA) on the selected data set across 10 separate runs. Green filled circles give the ARI achieved by the single one of these runs to achieve the highest likelihood under simulated annealing. We compare these results with the ARI of the labels generated via several comparisons: hypergraph modularity [7] (only for senate-bills), belief-propagation nonbacktracking spectral clustering [12], greedy dyadic modularity maximization [51], and dyadic Laplacian spectral clustering [49]; the latter two dyadic methods are denoted by hollow markers. Values for hypergraph modularity maximization and nonbacktracking spectral clustering are reported in [7] and [12], respectively.](/article/2607.16046/community_detection_boxplot.png "fig:")
In 7, we observe that CHILI simulated annealing outperforms both of the dyadic methods—greedy-modularity and spectral clustering—on senate-bills in all 10 runs in ARI. On
primary-school, simulated annealing performs comparably in ARI to the competitor methods, with a slightly lower average ARI compared to spectral clustering and greedy modularity’s results. However, simulated annealing performs worse on
high-school in ARI compared to its competitors. The relatively consistent performance of CHILI simulated annealing across these three data sets suggests that explicitly modeling dependence between edges may provide promising future directions
in community detection, motivating a need for much more efficient heuristics with which to carry out the optimization.
5 Discussion↩︎
We have proposed Copying Hyperedges Influenced by Label Interactions (CHILI), a generative model for the formation of binary-labeled hypergraphs in which edges form as noisy copies of previous edges to which extant and novel nodes are also added. The copying, extant, and novel node addition steps are all tuned by node labels and parameters which control the strength of homophily in each of these steps. We have demonstrated efficient parameter inference in the CHILI model via stochastic expectation maximization (SEM). We have also demonstrated that, in some circumstances, maximum-likelihood inference of node labels \(\mathbf{z}\) via simulated annealing, though computationally expensive, can outperform both dyadic and recent polyadic community detection algorithms that proceed from the assumption of conditional edge independence. We interpret this latter result as proof-of-concept that explicitly modeling edge-dependence can be a useful framework for community detection in hypergraphs.
There are several directions for future work. Most obviously, although our results show that maximum-likelihood inference in the CHILI model is promising for detecting certain kinds of community structure, the simulated annealing method we have used is
computationally impractical on data sets of even relatively modest size. Developing more efficient heuristics for maximizing the likelihood over \(\mathbf{z}\) will be an important direction for scalable application of the
CHILI model for community detection. Scalability may also imply opportunities for improved performance. We note, for example, that belief-propagation nonbacktracking spectral clustering [12] achieves an ARI of approximately 0.25 on the senate-bills data set (7), higher than CHILI. That algorithm requires alternating between parameter estimation and community detection phases. In principle, we could attempt a similar alternating procedure between SEM and simulated
annealing in the CHILI model, but the computational cost of simulated annealing makes this approach infeasible in practice. Future work to develop more efficient heuristics for community detection under our model could, in principle, enable such an
alternating algorithm, with potentially improved performance.
There are also opportunities to extend CHILI to more realistic settings and more diverse data sets. One obvious extension is to multinary node labels, which would enable multiway community detection algorithms. Another approach is to relax the assumption that every edge is formed through the edge-copying mechanism; one could instead allow some edges to form spontaneously without copying an earlier edge. A final possible extension is to incorporate recency-bias [22] into the copy mechanism, so that more recent edges are more likely to be copied while edges that are sufficiently old may cease to be copied at all. We hope that CHILI can provide a useful jumping-off point for future work in generative modeling and algorithms for hypergraphs with edge dependence mediated by node labels.
PSC acknowledges support from the National Science Foundation under DMS-2407058. Compute was provided by Middlebury College through support from the National Science Foundation under award 1827373.
Data Availability↩︎
The code used to generate the results in this paper is available at https://github.com/Violet-Ross/labeled-growth. Data sets are publicly available through the papers cited in 1.
Use of Generative Language Models↩︎
Generative language models were used to assist the authors in implementing the model, methods and experiments shown in this paper. The authors carefully reviewed and revised all code and associated outputs. The authors take full responsibility for the correctness of all theory, methods and results presented in this paper.
6 Power-Law Degree Distribution↩︎
We now derive the asymptotic power law degree distribution described in 2.1.2. Our derivation follows that of He et al. [24], which in turn is based on a derivation by Mitzenmacher [34].
At stationarity, the mean degree \(\left\langle d \right\rangle\) of a node satisfies \[\begin{align} \left\langle d \right\rangle = \frac{m}{n}\left\langle k \right\rangle \;, \end{align}\] where \(\left\langle k \right\rangle = \left\langle k_0 \right\rangle + \left\langle k_1 \right\rangle\) is the mean edge size (computable from the stationary distribution \(q\left( k_0,k_1 \right)\)), \(m\) is the number of edges, and \(n\) is the number of nodes. In expectation, there are \(\eta{} = \eta{}_- + \eta{}_+\) new nodes added per new edge, which means that \(\frac{m}{n} \to \frac{1}{\eta{}}\) as \(m, n \to \infty\) in expectation. This implies the asymptotic relation \[\begin{align} \left\langle d \right\rangle = \frac{\left\langle k \right\rangle}{\eta{}} \;. \label{eq:average-degree} \end{align}\tag{7}\]
Define the moments \[\begin{align} \mu_{ij} &= \left\langle \frac{k_i k_j}{k} \right\rangle \;. \end{align}\] We will first calculate the expected number of nodes of label 0 sampled in the edge-sampling step. Suppose that an edge \(f\) is sampled with \(k_0\) nodes of label 0 and \(k_1\) nodes of label 1, for a total of \(k = k_0 + k_1\) total nodes. Conditional on this event, the probability that a node of label 0 is selected as the seed node \(u\) is \(\frac{k_0}{k}\), and when this occurs, the expected number of nodes of label 0 in \(e\) is \(1 + \rho{}_+ (k_0 - 1)\). If, on the other hand, a node of label 1 is selected as the seed node \(u\), which occurs with probability \(\frac{k_1}{k}\), then the expected number of nodes of label 0 in \(e\) is \(\rho{}_- k_0\). So, the expected number of nodes of label 0 copied into the new edge \(e\) from \(f\) is \[\begin{align} \nu{}_0 &= \left\langle \frac{k_0}{k}\left( 1 + \rho{}_+(k_0 - 1) \right) + \frac{k_1}{k} \rho{}_- k_0 \right\rangle \\ &= \left\langle \frac{k_0}{k} \right\rangle + \rho{}_+ \left[ \mu_{00} - \left\langle \frac{k_0}{k} \right\rangle \right] + \rho{}_- \mu_{01}\;. \end{align}\] Similarly, the expected number of nodes of label 1 copied into the new edge \(e\) from \(f\) is \[\begin{align} \nu{}_1 = \left\langle \frac{k_1}{k} \right\rangle + \rho{}_+ \left[ \mu_{11} - \left\langle \frac{k_1}{k} \right\rangle \right] + \rho{}_- \mu_{01}\;. \end{align}\] The expected total number of nodes copied into the new edge \(e\) from \(f\) is then \[\begin{align} \nu{} = \nu{}_0 + \nu{}_1 = 1 + \rho{}_+(\mu_{00} + \mu_{11} - 1) + 2\rho{}_- \mu_{01}\;. \end{align}\] Importantly, \(\nu{}\) depends only on the parameters and the moments of the joint distribution \(q\left( k_0,k_1 \right)\), which can be computed via linear algebra as described in 2.1.1.
To proceed, we make a mean-field assumption that the degree \(d\) of a node is independent of the ratio \(\frac{k_0}{k}\) of label-0 nodes in the edge \(f\) sampled to form the new edge \(e\). Under this assumption, the probability that a single node selected in the sampling step has degree \(d\) is \(dr_{t}(d)\), where \(r_{t}(d)\) is the fraction of nodes with degree \(d\) at time \(t\); the mean-field assumption implies that the node participates in \(d\) edges, independently of the ratio \(\frac{k_0}{k}\) of label-0 nodes in the edge \(f\) sampled to form the new edge \(e\). The expected number of nodes of degree \(d\) sampled in the edge-sampling step is then \(\nu{} \frac{d r_{t}(d)}{\left\langle d \right\rangle}\), and the expected change in the count of degree \(d\) nodes through the edge-sampling step is \[\begin{align} \frac{\nu{}}{\left\langle d \right\rangle} \left[ (d-1)r_{t}(d-1) - d r_{t}(d) \right] \end{align}\] Similarly, the expected number of nodes that are sampled in the extant node addition step is \(\gamma{} = \gamma{}_+ + \gamma{}_-\). Since this step does not involve any selection proportional to degree, the expected change in the count of degree \(d\) nodes through this step is
\[\begin{align} n_{t+1}r_{t+1}(d) - n_{t}r_{t}(d) = \frac{\nu{}}{\left\langle d \right\rangle} \left[ (d-1)r_{t}(d-1) - d r_{t}(d) \right] + \gamma{} \left[ r_{t}(d-1) - r_{t}(d) \right] \;. \end{align}\]
We know that at stationarity, \(n_{t+1} - n_{t}\) is in expectation equal to \(\eta{} = \eta{}_+ + \eta{}_-\), the expected number of novel nodes, while at stationarity \(r_{t+1}(d) = r_{t}(d) = r(d)\) independent of \(t\). Thus, we obtain \[\begin{align} r(d) \eta{} = \frac{\nu{}}{\left\langle d \right\rangle} \left[ (d-1)r(d-1) - d r(d) \right] + \gamma{} \left[ r(d-1) - r(d) \right] \;. \end{align}\] Rearranging gives \[\begin{align} \frac{r(d)}{r(d-1)} &= \frac{\frac{\nu{}}{\left\langle d \right\rangle}(d-1) + \gamma{}}{\frac{\nu{}}{\left\langle d \right\rangle} d + \eta{} + \gamma{}} \\ &= 1 - \frac{\frac{\nu{}}{\left\langle d \right\rangle} + \eta{}}{\frac{\nu{}}{\left\langle d \right\rangle} d + \eta{} + \gamma{}} \;. \end{align}\] Allowing \(d\) to become large (which describes the tail of the degree distribution), we have \[\begin{align} \frac{r(d)}{r(d-1)} &\approx 1 - \frac{1}{d} \frac{\frac{\nu{}}{\left\langle d \right\rangle} + \eta{}}{\frac{\nu{}}{\left\langle d \right\rangle}} \\ &= 1 - \frac{1}{d} \left( 1 + \frac{\eta{}}{\nu{} / \left\langle d \right\rangle} \right) \\ &= 1 - \frac{1}{d} \left( 1 + \frac{\left\langle k \right\rangle}{\nu{} } \right)\;. \end{align}\] a relation which corresponds to a power law with exponent \(\zeta = 1 + \frac{\left\langle k \right\rangle}{\nu{}}\). Expanding out \(\nu{}\) yields 2 in the main text. Explicitly calculating this exponent requires the moments \(\left\langle k \right\rangle\), \(\mu_{00}\), \(\mu_{11}\), and \(\mu_{01}\), which can be computed from the stationary joint distribution \(q_{\infty} \left( k_0,k_1 \right)\).
7 Additional Details on Stochastic Expectation-Maximization (SEM)↩︎
7.1 Sufficient Statistics↩︎
We describe the function \(\boldsymbol{\psi}\) from 3.1 which computes the conditional sufficient statistics for the parameters \(\boldsymbol{\theta}= (\rho{}_+, \rho{}_-, \gamma{}_+, \gamma{}_-, \eta{}_+, \eta{}_-)\) given an observed edge \(e\), a candidate parent edge \(f\), a candidate seed node \(u\), and the node labels \(\mathbf{z}\). Let \(z = z_u\) be the label of the seed node \(u\). Let \(e_z\) be the set of nodes in \(e\) with label \(z\), and let \(f_z\) be the set of nodes in \(f\) with label \(z\). Then, the sufficient statistics for the parameters \(\boldsymbol{\theta}\) [53] are given by \[\begin{align} \boldsymbol{\psi}(e, f, u, \mathbf{z}) = \begin{pmatrix} \left| e_z \cap f_z \right| - 1 \\ \left| f_z \right| - 1 \\ \left| e_{\bar{z}} \cap f_{\bar{z}} \right| \\ \left| f_{\bar{z}} \right| \\ \left| (e_z \setminus f_z) \cap \mathcal{N}_{t} \right| \\ \left| (e_{\bar{z}} \setminus f_{\bar{z}})\cap \mathcal{N}_{t} \right| \\ \left| e_z \setminus \mathcal{N}_{t} \right| \\ \left| e_{\bar{z}} \setminus \mathcal{N}_{t} \right| \end{pmatrix}^\top \;. \end{align}\] As described in 3.1, the vector \(\mathbf{s}\) of expected sufficient statistics is computed as a weighted average of \(\boldsymbol{\psi}(e, f, u, \mathbf{z})\) over the belief distribution \(p(f, u | e; \hat{\boldsymbol{\theta}}_{\ell})\) over the latent variables \(f\) and \(u\). Then, the function \(g\) which computes the maximum likelihood estimates of the parameters \(\boldsymbol{\theta}\) from \(\mathbf{s}\) is given by \[\begin{align} g(\mathbf{s}) = \left( \frac{s_1}{s_2}, \frac{s_3}{s_4}, s_5, s_6, s_7, s_8 \right) \;. \end{align}\]
7.2 Learning Rate↩︎
We use an exponentially decaying learning rate for SEM with form \(\alpha^{[\ell]} = \alpha^{[0]} \cdot e^{-c\ell}\). The initial learning rate is set to \(\alpha^{[0]} = 0.01\), and the decay rate \(c = 0.001\). We initialized the vector of sufficient statistics \(\mathbf{s}^{[0]} = (1, 2, 1, 2, 0.5, 0.5, 0.5, 0.5)\). Other reasonable choices of hyperparameters did not substantially change the results of our experiments.
8 Synthetic SEM Parameter Values↩︎
In this section we give more details on the parameter values to demonstrate the performance of SEM on synthetic data in 4.
9 Approximation Argument for Simulated Anealing↩︎
This section provides further detail on the approximation to the likelihood described in 3.2 and used in our simulated annealing implementation. The likelihood of an edge \(e\) being generated in timestep \(t\) can be written
\[\begin{align} p(e, \mathbf{z}_e'|\mathbf{z}_{t-1};\boldsymbol{\theta}) = \sum_{\substack{f\prec e\\ u\in f}} p(e, \mathbf{z}_e' | f, u, \mathbf{z}_{t-1}; \boldsymbol{\theta}) p(f,u| \mathcal{E}_{t-1}, \mathbf{z}_{t-1}; \boldsymbol{\theta}) = \frac{1}{m_{t}}\sum_{f\prec e} \frac{1}{|f|} \sum_{u \in f\cap e} p(e, \mathbf{z}_e' | f, u, \mathbf{z}_{t-1}; \boldsymbol{\theta}) \;, \end{align}\]
where \(m_{t}\) is the number of edges at time \(t\). It is helpful to reindex this sum by the size of the intersection \(|f\cap e|\) between \(f\) and \(e\). \[\begin{align} p(e, \mathbf{z}_e'|\mathbf{z};\boldsymbol{\theta}) &= \sum_{h=0}^{|e|} \sum_{\substack{f\prec e\\ |f\cap e| = k}} \frac{1}{m_{t}}\sum_{u \in f\cap e} \frac{1}{|f|} p(e,\mathbf{z}_e' | f, u, \mathbf{z}; \boldsymbol{\theta}) \;. \end{align}\] We will now argue that, when \(n\) is large, the terms in this sum corresponding to the largest intersections (for which \(h = \operatorname*{argmax}_{f\prec e} |f\cap e|\)) will dominate the likelihood, allowing us to neglect edges \(f\) whose intersection with \(e\) is smaller.
Fix edge \(e\) generated at time \(t\) and candidate parent edge \(f\). We subscript \(e_z\), \(f_z\), and \({\mathcal{N}_{t}}_z\) to denote the subsets of \(e\), \(f\), and \(\mathcal{N}_{t}\) with label \(z\). We let \({n_{t}}_z = |{\mathcal{N}_{t}}_z|\) be the number of nodes of label \(z\) at time \(t\) and note that \({n_{t}}_0 + {n_{t}}_1 = n_{t}\). Let \(h_z = |(e_z \setminus f_z )\cap \mathcal{N}_z^t|\) be the number of extant nodes of label \(z\) in \(e\). Then, the extant node addition process appears in \(p(e,\mathbf{z}_e' | f, u, \mathbf{z}; \boldsymbol{\theta})\) as a factor \[\begin{align} E(h_0, h_1) = \frac{\frac{h_{z_u}^{\gamma_+} e^{-h_{z_u}}}{\gamma_+!}\times\frac{h_{\bar{z}_u}^{\gamma_-} e^{-h_{\bar{z}_u}}}{\gamma_-!}}{\dbinom{{n_{t}}_0 - |f_{0}|}{h_{0}}\dbinom{{n_{t}}_1 - |f_{1}|}{h_{1}}} \;, \end{align}\] where the numerator describes the likelihood of adding fixed numbers of extant nodes of each label and the denominator describes the likelihood of selecting the specific set of extant nodes from the set of all candidates. Under the assumption that \(n_{t} \gg |f|, |e|\), we can estimate \[\begin{align} \dbinom{{n_{t}}_z - |f_{z}|}{h_{z}} \approx \frac{({n_{t}}_z)^{h_z}}{h_z!} \approx \frac{(c_z n_{t})^{h_z}}{h_z!} \;, \end{align}\] where we have defined \(c_z = \frac{{n_{t}}_z}{n_{t}}\) as the fraction of nodes with label \(z\). Then, absorbing terms that do not depend on \(n_{t}\) into a constant term \(C(h_0, h_1)\), we have \[\begin{align} E(h_0, h_1) &= C(h_0, h_1) n_{t}^{-h_0 - h_1} + O(n_{t}^{-h_0 - h_1 - 1}) \\ &= C(h_0, h_1) n_{t}^{-h} + O(n_{t}^{-h - 1}) \;, \end{align}\] where \(h = h_0 + h_1 = |(e\setminus f) \cap \mathcal{N}_{t}|\) is the total number of extant nodes added to \(e\) and \[\begin{align} C(h_0, h_1) &= \frac{\frac{h_{z_u}^{\gamma_+} e^{-h_{z_u}}}{\gamma_+!}\times\frac{h_{\bar{z}_u}^{\gamma_-} e^{-h_{\bar{z}_u}}}{\gamma_-!}}{\frac{(c_0)^{h_0}}{h_0!}\frac{(c_1)^{h_1}}{h_1!}} \;. \end{align}\] It follows that \(p(e,\mathbf{z}_e' | f, u, \mathbf{z}; \boldsymbol{\theta})\) is of order \(O(n_{t}^{-h})\) as \(n_{t}\) grows large. Terms corresponding to choices of \(f\) for which \(h\) is larger (of which there may be up to approximately \(m_{t}\)) decay quickly as \(n_{t}\) grows large. We therefore have
\[\begin{align} p(e,\mathbf{z}_e'|\mathbf{z};\boldsymbol{\theta}) &= \sum_{\substack{f\prec e\\ |f\cap e| = h'}} \frac{1}{m_{t}}\sum_{u \in f\cap e} \frac{1}{|f|} p(e,\mathbf{z}_e' | f, u, \mathbf{z}; \boldsymbol{\theta}) + m_{t}O(n_{t}^{-h'-1}) \label{eq:sa-approx} \end{align}\tag{8}\]
where \(h' = \min_{f\prec e} |(e\setminus f) \cap \mathcal{N}_{t}|\) is the minimum number of extant nodes added to \(e\) across all choices of \(f\).
This argument suggests that, when \(n_{t}\) is relatively large in relation to \(m_{t}\), the likelihood of \(e\) is dominated by terms corresponding to choices of \(f\) for which the number of extant nodes added to \(e\) is small. There are two limitations to this argument in the context of empirical data sets. First, it may be the case that only a small number of candidate parent edges \(f\) have large intersections with \(e\) , which may lead to a poor approximation of the likelihood when we evaluate only a small number of corresponding terms. Second, several of our data sets have \(m_{t}\) large in relation to \(n_{t}\), making it dubious to neglect the error term. We therefore introduce two variations to our approximation heuristic:
Top \(j\) heuristic: Rather than evaluate only terms for \(f\) for which the number of extant nodes \(h\) to be added is minimal, we evaluate terms for the top \(j\) choices of \(f\) with the smallest values of \(h\).
Batch normalization: Rather than normalizing by \(\frac{1}{m_{t}}\) in 8 , we normalize by the number of choices for \(f\) that we evaluate for a specific \(e\). This may be \(j\) under the top \(j\) heuristic, or it may be the number of \(f\) choices with \(h\) equal to the minimum \(h'\) if we do not use the top \(j\) heuristic. Letting \(\mathcal{F}_e\) be the set of seed edges \(f\) evaluated for a given \(e\) under the top \(j\) or similar heuristic, our approximation is \[\begin{align} \tilde{p}(e,\mathbf{z}_e'|\mathbf{z};\boldsymbol{\theta}) &=\frac{1}{|\mathcal{F}_e|}\sum_{\substack{f\in \mathcal{F}_e}} \sum_{u \in f\cap e} \frac{1}{|f|} p(e,\mathbf{z}_e' | f, u, \mathbf{z}; \boldsymbol{\theta})\;. \end{align}\] The approximate likelihood of the complete data set is then \[\begin{align} \tilde{p} (\mathcal{H}, \mathbf{z}; \boldsymbol{\theta})&= \prod_{e\in \mathcal{E}}\tilde{p}(e,\mathbf{z}_e'|\mathbf{z};\boldsymbol{\theta})\;. \label{eq:sa-approx-2} \end{align}\tag{9}\]
It is important to note that the batch normalization 9 tends to overestimate the likelihood, since we have restricted ourselves to an average of high-probability seed-edges. For the purposes of simulated annealing, this is not a major problem since we are only interested in the relative likelihood of different labelings \(\mathbf{z}\) and not the absolute scale of likelihood.
10 Simulated Annealing for Community Detection↩︎
We now describe the simulated annealing algorithm from 3.2 in greater detail. Our custom learning schedule uses the first epoch to set the absolute scale of the acceptance probability decay, and we therefore describe this epoch as a separate algorithm (10). We then proceed to the main simulated annealing algorithm (11), which uses the scale of acceptance probability decay \(\sigma_\mathcal{H}\) computed in the first epoch to set the acceptance probability for label flips in subsequent epochs. We only allow the probability to begin to decay after the first four epochs in order to encourage exploration of the label space.
Here, \(\sigma(\boldsymbol{\delta})\) is the square root of the empirical variance \[\begin{align} \sigma(\boldsymbol{\delta}) = \sqrt{\frac{\sum_{i = 1}^n \delta_i^2 - (\frac{1}{n}\sum_{i=1}^n \delta_i)^2}{n-1}}\;. \end{align}\]
11 Simulated Annealing Illustration↩︎
We provide simple illustration of a single run of simulated annealing as described in [algocf:simulated_annealing_first_epoch,algocf:simulated_annealing] for community detection in a hypergraph generated from our model. In this simulation, we generated a
simple hypergraph according to our model with 67 nodes and 200 edges. Running simulated annealing for 1340 steps, we obtained a maximum log-likelihood of -10.38 nats per edge with an Adjusted Rand Index (ARI) against ground-truth of 0.39. Laplacian
spectral clustering as implemented in the Python package networkx [54] obtained an ARI of -0.03. The parameters for the model used to generate
the hypergraph are \[\begin{align} \boldsymbol{\theta}= (\rho{}_+, \rho{}_-, \gamma{}_+, \gamma{}_-, \eta{}_+, \eta{}_-) = (0.90{}, 0.10{}, 1.00{}, 0.25{}, 0.20{}, 0.10{})\;.
\end{align}\] For the purposes of this illustration, we perform community detection under the Nishimori condition, using a guess \(\hat{\boldsymbol{\theta}}\) equal to the true parameters \(\boldsymbol{\theta}\) used to generate the hypergraph.