Generative Hypergraph Clustering: From Blockmodels to Modularity

There are many methods for inference in stochastic blockmodels and their relatives, including variational coordinate ascent airoldi2008mixed, variational belief propagation decelle2011asymptotic; zhang2014scalable, and Markov Chain Monte Carlo nowicki2001estimation; peixoto2020merge. We perform approximate maximum-likelihood inference via coordinate ascent. We do so in order to exploit a recent connection between maximum-likelihood inference in stochastic blockmodels and the popular modularity objective for graph clustering newman2016equivalence. Our coordinate ascent framework, in which we alternate between estimating parameters and node labels, is a close relative of expectation-maximization (EM) algorithms for blockmodel inference decelle2011asymptotic. Standard versions of EM construct “soft” clusters, in which each node is given a weighted assignment in every possible cluster. “The” cluster for a given node is often taken to be the cluster in which the node has largest weight. In contrast, our approach generates “hard” clusters in which each node belongs to exactly one cluster. Profile likelihood methods offer an alternative framework for maximum-likelihood inference in SBMs bickel2009nonparametric, and their development for hypergraphs is another promising avenue of future work.

In the maximum-likelihood framework, we learn estimates $\hat{\mathbf{z}}$ of the node labels, $\hat{\Omega}$ of the affinity function, and $\hat{\bm{\theta}}$ of the degree parameters by solving the optimization problem

where $\mathbf{A}$ is a given data set represented by a collection of (integer-weighted) hyperedges. As usual, it is easier to work with the log-likelihood, which has the same local optima. The log-likelihood is

where

The first term $Q(\mathbf{z},\Omega,\bm{\theta})$ is the only part of the log-likelihood that depends on the group assignments $\mathbf{z}$ and affinity function $\Omega$. The second term depends on $\bm{\theta}$, while the third term depends only on the data $\mathbf{A}$ and can be disregarded for inferential purposes.

In the coordinate ascent approach to maximum-likelihood, we alternate between two stages. In the first stage, we assume a current estimate $\hat{\mathbf{z}}$ and obtain new estimates of $\Omega$ and $\bm{\theta}$ by solving

The resulting $\hat{\Omega},\hat{\bm{\theta}}$ can be viewed as maximum-likelihood estimates, conditioned on the current estimate $\hat{\mathbf{z}}$ of the label vector $\mathbf{z}$. In the second stage, we assume current estimates $\hat{\Omega}$ and $\hat{\bm{\theta}}$ and obtain a new estimate of $\mathbf{z}$ by solving

We alternate between these two stages until convergence.

There are several identifiability issues that must be addressed. First, permuting the group labels in $\mathbf{z}$ and $\Omega$ does not alter the value of the likelihood. We therefore impose an arbitrary order on group labels. Second, the number of possible groups $\bar{\ell}$ can in principle be larger than the number of groups present in $\mathbf{z}$. Such a case would correspond to the presence of groups which are statistically possible but empty in the given data realization. While other treatments are possible, we choose to disregard empty groups and treat $\bar{\ell}$ as equal to the number of distinct labels in an estimate of $\mathbf{z}$. A final form of unidentifiability relates to the scales of $\bm{\theta}$ and $\Omega$. For a fixed $\bm{\theta}$ and $\Omega$, we can construct $\bm{\theta}^{\prime}\neq\bm{\theta}$ and $\Omega^{\prime}\neq\Omega$ such that $\mathcal{L}(\mathbf{z},\Omega,\bm{\theta})=\mathcal{L}(\mathbf{z},\Omega^{\prime},\bm{\theta}^{\prime})$ (LABEL:sec:unidentifiable). To enforce identifiability, we must therefore place a joint normalization condition on either $\bm{\theta}$ or $\Omega$. We choose to constrain $\bm{\theta}$ such that

where $\text{{vol}}(\ell)=\sum_{i=1}^{n}d_{i}\delta(z_{i},\ell)$ and $d_{i}$ is the (weighted) number of hyperedges in which node $i$ appears.³³3Here, and throughout the rest of the text, $\delta$ is an indicator function evaluating to 1 if all its inputs are equal, and is 0 otherwise.

The usefulness of (9) is that, when $\mathbf{z}$ is known or estimated, the conditional maximum-likelihood estimates $\hat{\bm{\theta}}$ and $\hat{\Omega}$ in (7) take simple, closed forms. First, for a fixed label vector $\mathbf{z}$, when using the normalization (9), the maximum-likelihood estimate for $\bm{\theta}$ is

(See LABEL:sec:proof_constant.) Second, conditioned on $\mathbf{z}$, if $\Omega$ takes constant value $\omega$ on some set $Y$ of unordered tuples of labels,⁴⁴4In full generality, we can have one such $\omega$ for every possible label arrangement for each hyperedge size in the data. Later, we will make natural restrictions on $\ fOmega$. the maximum likelihood estimate for $\omega$ is

(See LABEL:sec:omega_estimate.) Although (10) assumes that $\mathbf{z}$ was fixed, it is not necessary to know $\mathbf{z}$ in order to form the estimate $\hat{\bm{\theta}}$. However, forming the estimate $\hat{\omega}$ via (11) requires that we know or estimate $\mathbf{z}$. It is therefore important to remember that $\hat{\omega}$ is not a globally optimal estimate, but rather a locally optimal estimate conditioned on the currently estimated group labels.

The formula (11) could also be inserted directly into the full likelihood maximization problem (2), eliminating the parameters corresponding to $\Omega$ and producing a lower-dimensional profile likelihood, which could then in principle be optimized directly. This approach has been successful for dyadic blockmodels bickel2009nonparametric, and the development of similar methods for hypergraph blockmodels would be of considerable interest. The advantage of our coordinate ascent framework is that we are able to develop fast heuristics for solving (8), by generalizing widely-used algorithms for graph clustering (Section 4). Solving problem (7) in our framework can also be interpreted as evaluating the profile likelihood for a fixed cluster vector $\mathbf{z}$, highlighting the relationship between these approaches.

We now turn to the problem of inferring the label vector $\mathbf{z}$. This problem leads naturally to a class of modularity-type objectives for hypergraph clustering.

We obtain an important class of objective functions by stipulating that $\Omega$ is symmetric with respect to permutations of node labels. In this case, $\Omega(\mathbf{z}_{R})$ depends not on the specific labels $\mathbf{z}_{R}$ in a given node tuple $R$, but only on the number of repetitions of each. Statistically, the corresponding DCHSBM generates hypergraphs in which all groups are statistically identical, conditioned on the degrees of their constituent nodes. Symmetric affinity functions thus give a flexible generalization of the planted partition stochastic blockmodel jerrum1998metropolis; condon2001algorithms to the setting of hypergraphs.

Define the function $\phi(\mathbf{z})=\mathbf{p}$, where $p_{j}$ is the number of entries of $\mathbf{z}$ in the $j$th largest group in $\mathbf{z}$, with ties broken arbitrarily. For example, if $\mathbf{z}=(1,1,4,1,2,3,2)$, then $\mathbf{p}=(3,2,1,1)$. We call $\mathbf{p}$ a partition vector. The symmetry assumption implies that $\Omega$ is a function of $\mathbf{z}_{R}$ only through $\mathbf{p}=\phi(\mathbf{z}_{R})$. Accordingly, we abuse notation by writing $\Omega(\mathbf{p})\equiv\Omega(\mathbf{z})$ when $\mathbf{p}=\phi(\mathbf{z})$.

We now define generalized cuts and volumes corresponding to a possible partition vector $\mathbf{p}$ for tuples of $k$ nodes:

where $\mathcal{R}^{k}$ is the subset of tuples in $\mathcal{R}$ consisting of $k$ nodes. The function $\text{{cut}}_{\mathbf{p}}(\mathbf{z})$ counts the number of edges that are split by $\mathbf{z}$ into the specified partition $\mathbf{p}$, while the function $\text{{vol}}_{\mathbf{p}}(\mathbf{z})$ is a sum-product of volumes over all grouping vectors $\mathbf{y}$ that induce partition $\mathbf{p}$. Let $\mathcal{P}$ be the set of partition vectors on sets up to size $\bar{k}$, the maximum size of hyperedges. We show in LABEL:sec:mod_cut_vol_deriv that the symmetric modularity objective can then be written as

For a partition vector $\mathbf{p}$ for tuples of $k$ nodes, direct calculation of $\text{{vol}}_{\mathbf{p}}(\mathbf{z})$ is a summation of $\overline{\ell}^{k}$ elements, which can be impractical when either $\overline{\ell}$ or $k$ are large. We show in LABEL:sec:eval_sums, however, that it is possible to evaluate these sums efficiently via a combinatorial identity. We also give a formula for updating volume terms $\mathrm{\text{{vol}}}_{\mathbf{p}}(\mathbf{z})$ when a candidate labeling is modified.

The objective function (14) is related to the multiway hypergraph cut problem studied by veldt2020hypergraph. They formulate the hypergraph cut objective in terms of splitting functions, which associate penalties when edges are split between two or more clusters. One then aims to minimize the sum of penalties subject to constraints that certain nodes must not lie in the same cluster. Symmetric affinity functions in our framework correspond to signature-based splitting functions in their terminology. Table 1 lists four of many families of affinity functions.

The All-Or-Nothing affinity function distinguishes only whether or not a given edge is contained entirely within a single cluster. This affinity function is especially important for scalable computation, and we discuss it further below. The Group Number affinity depends only on the number of distinct groups represented in an edge, regardless of the number of incident nodes in each one. Special cases of the Group Number affinity arise frequently in applications. When $f\left(\left\lVert\mathbf{p}\right\rVert_{0},k\right)=\text{exp}(\left\lVert\mathbf{p}\right\rVert_{0}-1)$, the first term of the modularity objective corresponds to a hyperedge cut penalty that is known in the scientific computing literature as “$\text{connectivity}-1$” deveci2015hypergraph, the $K-1$ metric karypis2000multilevel, or the boundary cut hendrickson2000graph. It has also been called fanout in the database literature kabiljo2017social. The related Sum of External Degrees penalty karypis2000multilevel is also a special case of the Group Number affinity. The Relative Plurality affinity considers only the relative difference between the size of the largest group represented in an edge and the next largest group. This rather specialized affinity function is especially appropriate in contexts where groups are expected to be roughly balanced, as we find, for example, in party affiliations in Congressional committees. Finally, the Pairwise affinity counts the number of pairs of nodes within the edge whose clusters differ. While this affinity function uses similar information to that used in dyadic graph models, there is no immediately apparent equivalence between any dyadic random graph and a DCHSBM with the Pairwise affinity function. There are many more symmetric affinity functions; see Table 3 of veldt2020hypergraph for several other splitting functions which can be used to define affinities.

An important subtlety was recently raised by zhang2020statistical concerning the relationship between blockmodels and modularity in newman2016equivalence, which also applies to our derivation of (14) above and (15) below. We derived the conditional maximum-likelihood estimates (10) and (11) of $\bm{\theta}$ and $\Omega$ under the assumption of a general, unconstrained affinity function $\Omega$. It is not guaranteed that these same estimates maximize the likelihood when additional constraints—such as the symmetry constraint $\Omega(\mathbf{z})=\Omega(\phi(\mathbf{z}))$—are imposed. Indeed, as zhang2020statistical show for the case of dyadic graphs, (10) and (11) for estimating $\hat{\bm{\theta}}$ and $\hat{\Omega}$ are only exact under the symmetry assumption on $\Omega$ when $\text{{vol}}(\ell)$ is constant for each $\ell\in[\hat{\ell}]$. When the sizes of groups vary, as is typical in most data sets, (10) and (11) are instead approximations of the exact conditional maximum-likelihood estimates. The situation is reminiscent of the tendency of the graph modularity objective to generate clusters of approximately equal sizes gleich2016mining. The objectives and algorithms we develop below should therefore be understood as approximations to coordinate-ascent maximum likelihood inference, which are exact only in the case that all clusters have equal volumes. See zhang2020statistical for a detailed discussion of these issues in the context of dyadic graphs.

All-Or-Nothing affinity function is of special interest for modeling and computation. This affinity function is a natural choice for systems in which the occurrence of an interaction or relationship depends strongly on group homogeneity.

Inserting the All-Or-Nothing affinity function from Table 1 into (14) yields, after some algebra (LABEL:sec:mod_cut_vol_deriv_aon), the objective

where $\beta_{k}=\log\omega_{k1}-\log\omega_{k0}$, $\gamma_{k}=\beta_{k}^{-1}(\omega_{k1}-\omega_{k0})$, and $J(\bm{\omega})$ collects terms that do not depend on the partition $\mathbf{z}$. We collect $\{\beta_{k}\}$ and $\{\gamma_{k}\}$ into vectors $\bm{\beta},\bm{\gamma}\in\mathbb{R}^{\bar{k}}$. We have also defined

In this expression, $m_{k}$ is the (weighted) number of hyperedges of size $k$, i.e., $m_{k}=\sum_{R\in\mathcal{R}^{k}}a_{R}$. The cut terms $\text{{cut}}_{k}(\mathbf{z})$ thus count the number of hyperedges of size $k$ which contain nodes in two or more distinct clusters. This calculation is a direct generalization of that from newman2016equivalence for graph modularity. Indeed, we recover the standard dyadic modularity objective by restricting to $k=2$. We call (15) the All-Or-Nothing (AON) hypergraph modularity.

Recently, Kami2018 proposed a “strict modularity” objective for hypergraphs. This strict modularity is a special case of (15), obtained by choosing $\omega_{k0}$ and $\omega_{k1}$ such that $\beta_{k}=1$ and $\gamma_{k}=\frac{m_{k}}{(\text{{vol}}(H))^{k}}$, where $\text{{vol}}(H)=\sum_{i=1}^{n}d_{i}$ is the sum of all node degrees in hypergraph $H$. However, leaving these parameters free lends important flexibility to our proposed AON objective (15). Tuning $\bm{\beta}$ allows one to specify which hyperedge sizes are considered to be most relevant for clustering. In email communications, for example, a very large list of recipients may carry minimal information about their social relationships, and it may be desirable to down-weight large hyperedges by tuning $\bm{\beta}$. Tuning $\bm{\gamma}$ has the effect of modifying the sizes of clusters favored by the objective, in a direct generalization of the resolution parameter in dyadic modularity reichardt2006statistical; veldt2018cc. Importantly, it is not necessary to specify the values of these parameters a priori; instead, they can be adaptively estimated via (11).

The first phase of our Symetric HMLL algorithm mirrors standard graph Louvain: nodes begin in singleton clusters and in turn greedily move to adjacent clusters until no more improvement is possible. Phase two of graph Louvain reduces edges between clusters into weighted edges between supernodes in a structure-preserving way. However, in the hypergraph case, naively collapsing clusters and hyperedges would discard important information about hyperedge sizes and the way each hyperedge is partitioned across clusters. Therefore, in subsequent stages of our algorithm, we greedily improve the objective by moving entire sets of nodes in the original hypergraph, rather than greedily moving individual nodes. In this way, a set of nodes that was assigned to the same cluster in a previous iteration is essentially treated as a supernode and moved as a unit, without collapsing the hypergraph and losing needed information about hyperedge structure.

Our overall procedure is formalized in Algorithms 1 and 2. Algorithm 1 visits in turn each set of nodes $S_{c}$ that represents a cluster $c$ from a previous iteration. The algorithm evaluates the change $\Delta Q$ in the objective function $Q$ associated with moving all of $S_{c}$ to an adjacent cluster, and then carries out the move that gives the largest positive change to the objective. This is repeated until moving a set $S_{c}$ can no longer improve the objective. The entire Symetric HMLL method is summarized by running Algorithm 2, which starts by placing every node in a singleton cluster before calling Algorithm 1 for the first time.

Algorithm 1 relies on a function $\Delta Q$ which computes the change in the objective $Q$ associated with moving all nodes from $S_{c}$ to an adjacent cluster. Changes to the second (volume) term in the objective can be computed efficiently using combinatorial identities (LABEL:sec:eval_sums, LABEL:prop:recursion). Changes to the first (cut) term require summing across all hyperedges incident to a node or set of nodes. At each hyperedge, we must evaluate the affinity $\Omega(\mathbf{p})$ on the current partition $\mathbf{p}$, as well as the affinity $\Omega(\mathbf{p}^{\prime})$ associated with the candidate updated partition $\mathbf{p}^{\prime}$. This situation contrasts with the case of the graph Louvain algorithm, in which it is sufficient to check whether a given edge joins nodes in the same or different clusters. The fact that we need to store and update the partition vector $\mathbf{p}$ for each hyperedge is what prevents us from collapsing a cluster of nodes into a monolithic supernode and recursively applying Algorithm 2 on a reduced data structure, as customary in graph Louvain. Thus, while clusters of nodes move as a unit in Algorithm 1 as well, it is necessary in this case to operate on the full adjacency data $\mathcal{A}$ at all stages of the algorithm. This can make Algorithm 2 slow on hypergraphs of even modest size. Developing efficient algorithms for optimizing the general symmetric modularity objective or various special cases is an important avenue of future work.

When $\Omega$ is the All-Or-Nothing affinity function, considerable simplification is possible. For each edge we need not compute the full partition vector $\mathbf{p}$ but only check whether or not $\left\lVert\mathbf{p}\right\rVert_{0}=1$. Rather than a general affinity function $\Omega$, we instead supply the parameter vectors $\bm{\beta}$ and $\bm{\gamma}$ appearing in (15). This allows us to compute on considerably simplified data structures. In particular, we are able to follow the classical Louvain strategy of collapsing clusters into single, consolidated “supernodes,” and restrict attention to hyperedges that span multiple supernodes. Because we do not need to track the precise way in which the hyperedges span multiple supernodes, we can forget much of the original adjacency data $\mathcal{A}$ and instead simply store the edge sizes of the hypergraph. These simplifications enable both significant memory savings and very rapid evaluation of the objective update function $\Delta Q$. We provide pseudocode for exploiting these simplifications in LABEL:sec:AON_louvain.

Like most Louvain-style modularity methods, the user cannot directly control the number of clusters returned by HMLL. One approach to influence the number of clusters is to manually set values for the affinity function $\Omega$ or the parameters $\bm{\beta}$ and $\bm{\gamma}$ (if using the All-Or-Nothing affinity). Rather than inferring these parameters from data, one can set them a priori and perform a single round of optimization over $\mathbf{z}$. This approach generalizes the common treatment of the resolution parameter in dyadic modularity maximization as a tuning hyperparameter rather than a number to be estimated from data reichardt2006statistical. Considerable experimentation may be required in order to obtain the desired number of clusters, and retrieving an exact number may not be possible.

Another approach to influencing the number of clusters is to impose a Bayesian prior on the community labels. In the simplest version of a Bayesian approach, one assumes that each node is independently assigned one of $\bar{\ell}$ labels with equal probability, prior to sampling edges. The probability of realizing a given label vector $\mathbf{z}$ is then $\bar{\ell}^{-n}$, which generates a term of the form $-n\log\bar{\ell}$ in the log-likelihood $\mathcal{L}$. This term may then be incorporated into Louvain implementations, with the result that greedy moves which reduce the total number of clusters $\bar{\ell}$ are strongly incentivized. The resulting regularized algorithm may then label vectors $\mathbf{z}$ with slightly smaller numbers of distinct clusters. This can be useful when it is known a priori that the true number of clusters in the data is small. Our implementation of AON HMLL incorporates this optional regularization term. We use this term only in the synthetic detectability experiments presented in LABEL:fig:detectability.

Dyadic Louvain algorithms are known for being highly efficient in large graphs. Here, we show that AON HMLL can achieve similar performance on synthetic data to Graph MLL (GMLL), a variant of the standard dyadic Louvain algorithm in which we return the combination of resolution parameter and partition that yield the highest dyadic likelihood. For a fixed number of nodes $n$, we consider a DCHSBM-like hypergraph model with $\bar{\ell}=n/200$ clusters and $m=10n$ hyperedges with size $k$ uniformly distributed between 2 and 4. Each $k$-edge is, with probability $p_{k}$, placed uniformly at random on any $k$ nodes within the same cluster. Otherwise, with probability $1-p_{k}$, the edge is instead placed uniformly at random on any set of $k$ nodes. We use this model rather than a direct DCHSBM to avoid the computational burden of sampling edges at each $k$-tuple of nodes, which is prohibitive when $n$ is large. For the purpose of performance testing, we compute estimates of the parameter vectors $\bm{\beta}$ and $\bm{\gamma}$ (in the case of AON HMLL) and the resolution parameter $\gamma$ (in the case of GMLL) using ground truth cluster labels. We emphasize that this is typically not possible in practical applications, since the ground truth labels are not known. We make this choice in order to focus on a direct comparison of runtimes of each algorithm in a situation in which both can succeed. In later sections, we study the ability of HMLL and GMLL to recover known groups in synthetic and empirical data when affinities and resolution parameters are not known.

Section 5.1 shows runtime, Adjusted Rand (ARI) and number of clusters returned on synthetic hypergraphs when $p_{2}=0.6$, $p_{3}=1/n^{3}$ and $p_{4}=1/n^{4}$. These parameter are chosen so that hyperedges of size three and four are rarely (if ever) contained completely inside clusters. Thus, hyperedges of different sizes provide different signal regarding ground truth clusters. For this experiment, we implemented Graph MLL by computing a normalized clique projection, in which nodes are joined by weighted dyadic edges with weights

We also performed experiments on an unnormalized clique projection with $w_{ij}=\left|e:i,j\in e\right|$, but do not show these results because, in this experiment, the associated MLL algorithm consistently fails to recover labels correlated with the planted clusters.

On smaller instances, HMLL outperforms Graph MLL in recovering planted clusters, as measured by the ARI. For larger instances, the recovery results are comparable. Interestingly, although GMLL and HMLL obtain similar accuracy in this experiment, they do so in different ways, with HMLL tending to generate more, smaller clusters than GMLL. Importantly, the runtimes are nearly indistinguishable, indicating that dyadic clique projections are necessary neither for accuracy nor for performance. We observed other choices of the parameters $p_{2}$, $p_{3}$, and $p_{4}$ in which HMLL substantially outperformed GMLL in cluster recovery and vice versa; however, in each case the algorithms’ respective runtimes tended to differ by only a small constant factor.

Interestingly, in this synthetic experiment, a combination of the two algorithms leads to the strongest recovery results. In addition to running each algorithm independently, we also ran a two-stage algorithm in which GMLL is used to generate an intermediate partition, and then HMLL is used to refine it. This procedure is similar to the warmstart approach from kaminski2020community. We emphasize again that these results are obtained on synthetic hypergraphs with pre-optimized affinity parameters, and so the effectiveness of the refinement strategy may not generalize to real data sets. Indeed, in the experiments on empirical data shown in LABEL:sec:empirical, we do not show results from the refinement procedure because the output partition was in each case essentially indistinguishable from the output of the dyadic partition. This may reflect the fact that we did not allow the algorithms to learn a priori the affinity parameters associated with the true data labels. Further investigation into the performance of hybrid strategies would be of considerable practical importance.

Informally, an algorithm is able to detect communities in a random graph model with fixed labels $\mathbf{z}$ when the output labeling $\hat{\mathbf{z}}$ of that algorithm is, with probability bounded above zero, correlated with $\mathbf{z}$. Using arguments from statistical physics, decelle2011asymptotic conjectured the existence of a regime in the graph stochastic blockmodel in which no algorithm can successfully detect communities. This conjecture has since been refined and proven in various special cases; see abbe2017community for a survey. In the dyadic stochastic blockmodel with two equal-sized planted communities, a necessary condition for detectability in the large-graph limit is

where $c_{i}$ is the mean number of within-cluster edges attached to a node, and $c_{o}$ is the mean number of between-cluster edges attached to a node. If this condition is not satisfied, no algorithm can reliably detect communities in the associated graph stochastic blockmodel, even though the communities are statistically distinct. This bound limits direct inferential methods, such as Bayesian or maximum-likelihood techniques, as well as methods based on maximization of modularity or other graph objectives nadakuditi2012graph. Several recent papers have considered the detectability problem in the case of uniform hypergraphs ghoshdastidar2014consistency; Ghoshdastidar2017; Angelini2016. In our model, the presence of edges of multiple sizes complicates analysis. Here, we limit ourselves to an experimental demonstration that the regimes of detectability for the graph SBM and our DCHSBM can differ significantly.