{'fraction_non_alphanumeric': 0.03947607549913026, 'fraction_numerical': 0.011799984422462809, 'mean_word_length': 4.3123922488104265, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nonbacktracking spectral clustering of nonuniform hypergraphs

Spectral methods offer a tractable, global framework for clustering in graphs via eigenvector computations on graph matrices. Hypergraph data, in which entities interact on edges of arbitrary size, poses challenges for matrix representations and therefore for spectral clustering. We study spectral clustering for nonuniform hypergraphs based on the hypergraph nonbacktracking operator. After reviewing the definition of this operator and its basic properties, we prove a theorem of Ihara-Bass type which allows eigenpair computations to take place on a smaller matrix, often enabling faster computation. We then propose an alternating algorithm for inference in a hypergraph stochastic blockmodel via linearized belief-propagation which involves a spectral clustering step again using nonbacktracking operators. We provide proofs related to this algorithm that both formalize and extend several previous results. We pose several conjectures about the limits of spectral methods and detectability in hypergraph stochastic blockmodels in general, supporting these with in-expectation analysis of the eigeinpairs of our studied operators. We perform experiments in real and synthetic data that demonstrate the benefits of hypergraph methods over graph-based ones when interactions of different sizes carry different information about cluster structure.

then be used to obtain clusters from the embedding space defined by these eigenvectors. There are many possibilities for the choice of matrix M. These include the graph adjacency matrix [57], various Laplacian matrices [65,72], the modularity matrix [58], and the nonbacktracking matrix [43]. While greedy methods such as the famous modularity-maximizing Louvain algorithm [14] are guaranteed only to find locally optimal solutions to cluster optimization problems, many spectral methods can be interpreted as approximations of globally optimal solutions to such problems.
In this paper we study extensions of nonbacktracking spectral clustering to the setting of hypergraphs. Several spectral methods, including nonbacktracking methods, exist for uniform hypergraphs-in which all edges contain the same number of nodes. Many interesting hypergraph data sets, however, are nonuniform, containing edges of multiple sizes [9]. Our extension focuses on the challenges and opportunities posed by such data. In addition to the development of new spectral algorithms, we also contribute conjectures related to the possibility of community detection in generative random hypergraph models.
The paper proceeds as follows. In Section 2, we survey related work on the dyadic nonbacktracking operator and its uses in graph data science, especially including its role in results on the theory of the detectability of communities in graphs. We also discuss spectral clustering techniques applied to hypergraphs, and discuss the hypergraph nonbacktracking operator B introduced by Storm [68]. In Section 3 we prove a theorem of Ihara-Bass type relating the spectrum of the nonbacktracking operator to that of a related matrix which is usually smaller. We also provide a first clustering algorithm based on this matrix. In Section 4 we describe a hypergraph stochastic blockmodel (HSBM) and state in-expectation eigenrelations for the nonbacktracking operator. Then, in Section 5, we study the relationship between the nonbacktracking operator and the belief-propagation algorithm in the HSBM. We prove a precise statement of the known heuristic [43,4] that the stability of an uninformative fixed point of the belief-propagation dynamics can be studied via the nonbacktracking operator. We also prove a relationship of Ihara-Bass type between the Jacobian governing the stability of this fixed point and a smaller matrix, often enabling faster computations. We use this relationship to propose an alternating spectral clustering algorithm based on belief-propagation. In Section 6, we pose several conjectures on the spectral properties of hypergraph nonbacktracking operators, and derive from them conjectured thresholds bounding the ability of our proposed algorithms to detect clusters. We support these conjectures with experiments on synthetic data. We move on to several empirical data sets in Section 7, finding that hypergraph spectral clustering using the belief-propagation Jacobian outperforms methods based on the projected graph when edges of different sizes play statistically distinct roles. We conclude in Section 8 with a discussion of limitations of our algorithm and suggestions for future work.
2. Related Work. We now introduce our primary notation and briefly discuss related work in clustering methods for graphs and hypergraphs. Table 1 gives a summary of our notation.

Notation and Preliminaries.
A hypergraph H is a tuple (N , E) consisting of a node set N and an edge set E. The node set N contains n ∶= N nodes. The edge set E consists of subsets of N . For each k in a set K of possible edge sizes, we let E k denote the set of edges of size k. We let ∂ k i ⊂ E k denote the set of edges of size k containing i, and we let ∂i = ∪ k∈K ∂ k i. We let m k = E k , and m = ∑ k∈K m k . We usually assume K = {2, 3, . . . ,k} for some maximum edge sizek, and set κ ∶= K . A hypergraph is k-uniform if we may take K = {k}, and nonuniform if it is not k-uniform for any k. A graph is a 2-uniform hypergraph; in this case we write G instead of H. We let ⟨k⟩ = 1 m ∑ k∈K km k be the empirical average edge size.
We will consider many linear maps defined on structures related to graphs and hypergraphs. Accordingly, if S is a finite set, we let V (S) be the vector space of formal sums of the form ∑ s∈S sa s for scalars a s . For integer j, we define the notation V (jS) ∶= V (S) ⊕ V (S)⋯ ⊕ V (S) with j direct summands. Given two finite sets S and T , we let L(S, T ) be the set of linear maps between V (S) and V (T ). We abbreviate L(S) ∶= L(S, S). We speak of elements of L(S, T ) interchangeably as either linear maps or matrices. In the latter case, we regard S and T as bases of V (S) and V (T ), respectively. So defined, V (S) is isomorphic to R S and L(S, T ) is isomorphic to R S × T . Our notation is intended to emphasize the relationship between the many linear maps we will encounter and the structures on which they act.
Elements of V (S) are written in lowerbase bold: v ∈ V (S). An entry of v is v s for some s ∈ S. Elements of L(S, T ) are written in uppercase bold: A ∈ L(S, T ). An entry of A is a s,t for s ∈ S and t ∈ T . In many cases we will need to consider multiply-indexed structures; for example, an element u ∈ V (S × T ) has elements of the form v st for s ∈ S and t ∈ T . An element C ∈ L(S × T, S ′ × T ′ ) has elements of the form c st,s ′ t ′ . When considering an indexed family of matrices such as {A 1 , . . . , A k }, we separate the family index with a semicolon when writing the entries, e.g., a k;i,j is the (i, j)-th entry of A k . When later considering objects indexed by community labels, we use upper indexing with parentheses. For example, a vector c has typical entry c (s) when s is a community label, and a matrix C has entries c (s,t) when t is also a community label. We use 1[P ] to denote the indicator function of the proposition P , and the shorthand δ i,j ∶= 1[i = j].
We now define the nonbacktracking operator on hypergraphs. This definition is due to Storm [68].
Definition 2.1 (Pointed Line Graph). A pointed edge iQ in a hypergraph H consists of an edge Q ∈ E along with a choice of point i ∈ Q. Define ⃗ E to be the set of pointed edges. Let ⃗ ∂ k i be the set of pointed edges of size k with point i, and let ⃗ ∂i ∶= ∪ k∈K ⃗ ∂ k i We say that pointed edge jR follows pointed edge iQ, written iQ → jR, if j ∈ Q ∖ i and Q ≠ R. The pointed line graph P of H is a directed graph whose nodes are elements of ⃗ E. There is a directed edge (iQ, jR) in P if iQ → jR. Definition 2.2 (Nonbacktracking Operator). The nonbacktracking operator B ∈ L( ⃗ E) associated to a hypergraph H is the directed adjacency operator of the pointed line graph P. Its entries are b iQ,jR ∶= 1[iQ → jR].

Nonbacktracking Methods in Network Data
Science. The nonbacktracking matrix B has found several applications in the the study of graphs (i.e., 2-uniform hypergraphs). Alon et al. [3] show that random walks on graphs governed by the nonbacktracking matrix B mix more rapidly than random walks governed by the adjacency matrix A, implying that these walks may be more efficient for graph exploration tasks. Martin et al. [49] propose an eigenvector centrality measure based on B, and show that this centrality avoids the pathological localization of classical adjacency-based eigenvector centrality in sparse graphs. Torres et al. [70] and Mellor and Grusovin [52] impose metrics on the space of point clouds in the complex plane. These metrics enable the comparison of the spectra of two nonbacktracking operators, and induce a pseudometric on simple graphs. The resulting pseudometrics can then be used for graph clustering tasks. Torres et al. [71] develop perturbation theory for the eigenvalues of B in order to identify influential nodes in spreading processes on graphs.
The nonbacktracking matrix plays an important role in the community detection problem on graphs. In a simple graph with two planted clusters, the eigenvectors of B can be used to assign labels to nodes correlated with the true clusters [43]. We refer to this approach to community detection as nonbacktracking spectral clustering (NBSC).
A standard way to study the theoretical behavior of many clustering algorithms, including NBSC, is to consider their behavior under generative models of random clustered graphs. A common such model is the sparse binary planted partition model. We begin with two groundtruth labeled communities of n 2 nodes. We then draw edges between nodes independently in such a way that each node has, in expectation, a edges joining it to other nodes in its own community and b edges joining it to nodes in the opposite community. 1 The expected performance of a given algorithm can then be studied as a function of the parameters a and b, usually as n → ∞.
Several spectral algorithms admit statistical guarantees for the recovery of planted clusters in the n → ∞ limit [73,44]. There are also important limitations. Heuristically, an algorithm detects communities in a generative model if that algorithm is able to reliably return labels that have better-than-random correlation with ground truth. A series of deep results [57,43,16] have shown that various spectral graph clustering methods, including NBSC, are able to detect communities in the binary planted partition model as n → ∞ if and only if The quantity φ may be viewed as a signal-to-noise ratio for the clustering problem. The general structure of such results is to show that, with high probability as n grows large, the graph matrix used for clustering possesses an automatically selectable eigenvalue whose eigenvector can be used to approximately recover the two planted communities. Formal concentration results along these lines for the nonbacktracking matrix B are recently available [16]. These results show that, when φ > 1, with high probability, the second eigenvalue of B is real and holds cluster information. The equation φ = 1 is often called the detectability threshold for the sparse binary planted partition model. It was conjectured [24] and later proven [53,54,50] that the condition φ > 1 is both necessary and sufficient for the existence of algorithms (of any kind) that detect communities in this model. The analysis of spectral methods can thus offer insight on the conditions under which the cluster detection task is possible at all.

Spectral Clustering Methods for Hypergraphs.
There is a wide range of algorithms for clustering and community detection in hypergraphs; see Chodrow et al. [20] for a recent overview. We focus here on spectral methods, i.e., methods which make use of the eigenvectors of some object associated to the hypergraph. One approach begins by replacing the hypergraph with a dyadic graph. This is frequently done by clique-projection, in which each edge e of size k is replaced by a clique of k 2 , possibly-weighted 2-edges between each pair of nodes contained in e. Graph spectral methods can then be applied [76]. Multiple versions of this approach possess asymptotic consistency guarantees under hypergraph planted partition models [35,29]. A challenge for such approaches is the need to treat edges of multiple sizes. In order to realize the desired statistical guarantees it is usually necessary to assume either a k-uniform hypergraph or a generative model in which edges of varying sizes carry similar information about community structure. Both assumptions can be restrictive in practice.
Tensor methods can provide explicit representations of polyadic relationships. A k-uniform hypergraph can be represented as a symmetric adjacency k-tensor A. In the k = 3 case, for example, we would have entry a i,j,h = 1 iff {i, j, h} ∈ E. There are several spectral methods for clustering k-uniform hypergraphs, many of which rely on concepts connected to tensor eigenvalues and eigenvectors [45]. Ke et al. [40] use a normalized tensor power iteration to compute eigenvectors, while Chang et al. [18] take an explicit optimization approach. Hu and Wang [36] derive a clustering method for uniform hypergraphs based on angular separability, and provide several statistical guarantees. The representation of hypergraphs using adjacency tensors is constrained by the need to represent all edges in a tensor of fixed dimensions; in particular, it is unclear how edges of multiple sizes should be represented in the same tensor. A proposal was recently made in this direction by Galuppi et al. [34], but the applications of the resulting tensor for data analysis problems remains to be explored. We also note recent work by Mulas and collaborators developing spectral theory for hypergraphs via tensors [34], random walks [56], and Laplace operators [55,38]. To our knowledge, the application of these techniques in the data scientific context has not yet been studied.
The use of the nonbacktracking operator B for clustering uniform hypergraphs was considered by Angelini et al. [4]. These authors studied properties of this operator in the context of uniform hypergraphs, including a weakened version of the Ihara-Bass theorem and some conjectures regarding the locations of the eigenvalues of the operator. The authors offered a conjecture on the detectability threshold in a uniform hypergraph planted partition model. Several of these conjectures were recently formalized and proved by Stephan and Zhu [67].
3. An Ihara-Bass Theorem for Nonuniform Hypergraphs. A theorem often attributed to both Ihara [37] and Bass [5], originally formulated for graphs, relates the spectrum of B to the spectrum of a different matrix B ′ which is usually smaller. We now extend this theorem to nonuniform hypergraphs. Our result generalizes a computation by Angelini et al. [4] and theorem by Stephan and Zhu [67].
It is useful to distinguish blocks of B by edge size. Let ⃗ E k be the set of pointed edges of size k, let ⃗ m k = ⃗ E k , and let ⃗ m = ∑ k∈K ⃗ m k .
The k-th nonbacktracking operator B k ∈ L( ⃗ E) associated to hypergraph H has the block form where 0 m 1 ×m 2 is the matrix of zeros of specified dimensions.
By construction, B = ∑ k∈K B k . Because the only indices in B k permitted to be nonzero in both rows and columns are the indices corresponding to B k→k , any eigenvalue of B k must either be zero or an eigenvalue of B k→k .
We now define adjacency and degree operators. We also distinguish these by edge size.
Definition 3.2 (Adjacency and Degree Operators). The k-th adjacency operator A k ∈ L(N ) associated to H has entries The k-th degree matrix D k ∈ L(N ) associated to H has entries

Define the block matrices
where A k and D k are as in Definition 3.2. Let K be the square matrix of size κ × κ with entries (2, 3, . . . ,k) along the diagonal, and zeroes everywhere else. Theorem 3.3 (Ihara-Bass for nonuniform hypergraphs). For any hypergraph H, we have and ⊗ is the Kronecker product.
We supply a proof of Theorem 3.3, as well as of Corollary 3.4 and Lemma 3.5 below, in Appendix A. The computational significance of Theorem 3.3 is that we can compute the interesting eigenvalues of B via a different matrix B ′ . This matrix is usually smaller than B, although it can also be less sparse. • For each k, if m k > n, then β = 1 − k is an eigenvalue of B with algebraic multiplicity at least m k − n. • If ∑ k∈K m k (k − 1) > κn, then β = 1 is an eigenvalue of B with algebraic multiplicity at least ∑ k∈K m k (k − 1) − κn.
• The remaining eigenvalues of B are eigenvalues of the matrix Corollary 3.4 expresses a relationship between the eigenvalues of B and B ′ . There is an associated relationship between their eigenvectors.
are doubly-indexed vectors defined entrywise by Then, B ′ x = βx. In particular, either x is an eigenvector of B ′ with eigenvalue β or x = 0.
Heuristically, Lemma 3.5 states that one can aggregate cluster information in u to the level of nodes by summing over all edges of size k with specified points. To complete a first spectral clustering algorithm, we further sum x 2 over edge sizes, obtaining vectorx with entries In the case of uniform hypergraphs, there exist recent probabilistic guarantees ensuring thatx is correlated with planted communities under a certain choice of data generating process [67]. Extending these guarantees to the nonuniform cases is an avenue of future work. We find experimentally that the sign ofx i is most informative, and it is therefore helpful to use the vectorx ∶= sgn(x), with the sign computed entrywise. 2 We now state our first spectral clustering algorithm, Nonbacktracking Hypergraph Spectral Clustering (NBHSC) in Algorithm 3.1. We compute a desired number of eigenvectors u ( ) , formx ( ) for each, and then cluster in the Euclidean embedding described by the set x ( ) . There are multiple choices for the Euclidean clustering subroutine. We assume that this subroutine accepts as input a matrix giving the coordinates of a point cloud in Euclidean space, and returns a vector of labels z. Throughout this paper we use the standard k-means algorithm [47,64], but other choices could in principle lead to superior performance.
Carrying out NBHSC requires the user to choose h, the number of eigenvectors with real eigenvalues to extract from B ′ . Analogy with the uniform hypergraph case [67] would suggest the use of h − 1 eigenvectors to cluster into h communities, provided that there are indeed h − 1 such eigenvectors with eigenvalues separated from the bulk other than an uninformative eigenvector corresponding to the largest real eigenvalue. We will later argue that NBHSC is limited in cases in which edges of different sizes carry different cluster information in H and that no choice of h can circumvent this limitation. To highlight these limitations and build a foundation for further development, we now discuss a generative model of random hypergraphs and study the spectral structure of B under this model. Computex ( ) via (3.3) and (3.4) 5:X ⋅ ← sgn(x ) 6: end for 7: z = Cluster(X) 8: return z 4. The Sparse Hypergraph Stochastic Blockmodel. In this section we briefly review the hypergraph stochastic blockmodel (HSBM), a generative model of clustered hypergraphs. Our choice of notation and formulation most closely resembles that of Ke et al. [40]. Many other related formulations exist in the literature [35,4,20,29,67,21]. We prove in-expectation results for eigenpairs of the matrices B k , and pose conjectures generalizing recent proofs by Stephan and Zhu [67] of eigenpair concentration results in the uniform case. These conjectures will also inform our development of belief-propagation hypergraph spectral clustering in Section 5.
Our blockmodel is a probability distribution over hypergraphs. We denote this distribution by η. To sample from this model, we first assign each node i ∈ N a label z i from a finite label alphabet Z of size . These labels are drawn independently from a probability vector q ∈ V (Z), so that q (s) gives the probability that z i = s for each i ∈ N . We collect the labels in a vector z. Let R give the set of possible edges, which we usually take to be sets of nodes with some specified possible sizes. Let R k denote the set of node subsets of size k. Define R(i) to be the set of subsets containing node i, and let R k (i) = R k ∩ R(i). To realize edges, we consider each set of nodes R ∈ R and add this set to the edge set E with probability η(R ∈ E z R ) depending on the labels in z R . The structure of the model is specified by the functional form of η. Here, we will focus on the sparse setting in which (4.1) η(R ∈ E z R ) = ω(z R ) n R −1 for some function ω that does not depend on n. This structure imposes sparsity on the hypergraph; the number of k-edges for a given node is asymptotically constant with respect to n. The overall probability to realize a given combination of label vector z and edge set E is The function ω controls the relative rates of edges within and between communities, and is an extension of the parameters a and b in the graph planted partition model described in Section 1. A common choice of the function ω is the "all-or-nothing" affinity, defined as When a k > b k , edges of size k form at higher rates within planted communities than between them. This choice of ω is also the partition used by Chodrow et al. [20] to derive a fast modularity maximization algorithm for hypergraphs. However, many other choices for ω are also possible, including ones that favor edge formation between communities rather than within them. Our work is independent of the choice of ω, subject only to the assumptions described in Subsection 4.1.
In our development below, it will be useful to reason asymptotically about many quantities. Define x ≐ y if there exists some constant r > 0 such that x = (1 + O(n −r ))y either deterministically or with high probability with respect to the blockmodel η as n grows large and the largest edge-size k remains fixed. Let This is the asymptotic expected number of k-edges attached to a given node i in cluster s, as can be verified via direct calculation: We also define which gives the expected number of cluster t neighbors of a node in cluster s through edges of size k. We have the identity which computes c (s) k by conditioning on the label of a second node in each possible edge. 4.1. Spectral Structure of B in the HSBM. We now aim to describe the eigenvectors of B k which are correlated with community structure under the blockmodel (4.2). We implicitly condition on a realization of the label vector z for simplicity, so that the only remaining randomness is in the edge generation process. We impose the following assumptions: • There are exactly two groups, labeled 1 and 2, of equal expected size. Thus, q (1) = q (2) = 1 2 . We further assume that the empirical distribution of group labels in the realized label vector z is close to q; this holds with high probability as n grows large.
• The group-specific expected degrees are equal: c for all k. Formal concentration results under these hypotheses are available for graphs [16] and kuniform hypergraphs [67]. While we anticipate that many of the techniques used for these cases will generalize to nonuniform hypergraphs, pursuing formal proofs is beyond our present scope. We instead provide informal results by reasoning in expectation.
Then, for all tuples Q and nodes i ∈ Q, we have The full proof of (4.7) is provided in Appendix C; the proof of (4.6) is ommitted but similar.
In expectation, Theorem 4.1 says that, B k has a Perron eigenvector with eigenvalue α k , and a community-correlated eigenvector with eigenvalue β k . To substantiate the claim that v is indeed community-correlated, we can sum over k−edges incident to i, obtaining This calculation is similar to the proof of Lemma 3.5. Indeed, this calculation shows that E[x 2;k ] ≐ β k σ with x 2;k as in Lemma 3.5, explicitly showing why the vector x 2 can be expected to be correlated with planted communities. In practice, we do not have access to the expectation and only to our realization of x 2 , and therefore our estimate of σ, is noisy. Formal concentration results along the lines of Bordenave et al. [16] and Stephan and Zhu [67] are necessary to provide probabilistic guarantees. From Theorem 4.1 we can also obtain in-expectation eigenpair relations for the full matrix B.
Proof. The corollary follows from the formula B = ∑ k∈K B k and the fact that the vectors u and v do not depend on k.

Implications for NBHSC.
In Algorithm 3.1, we form the nonbacktracking operator B, obtain a spectral embedding of the nodes, and then cluster in the embedding space. The spectral embedding is obtained by considering the h−1 eigenvectors of B with real eigenvalues largest in magnitude, excepting the uninformative eigenvector u. A necessary condition for NBHSC to perform better-than-random guessing is that at least one of these eigenvectors is correlated with ground-truth community labels. In expectation, 80% of both 2-edges and 3-edges contain nodes within a single cluster while 20% contain nodes from both clusters. The uninformative real eigenvalue with largest magnitude and community-informative eigenvalue are highlighted (far right). Ticks on the bottom margin give the theoretical predictions α and β for these two eigenvalues given by Corollary 4.2. In this case, spectral clustering using the eigenvector v associated to the second eigenvalue in magnitude achieves Adjusted Rand Index (ARI) 0.98 against ground truth. The solid circle has radius √ α. (b): As in (a), but only 10% of 2-edges and 50% of 3-edges are within-cluster. A community-correlated eigenvalue is again present, but lies inside the bulk of the spectrum. Clustering based on the associated eigenvector would achieve ARI 0.87 against ground truth, but it is not possible to automatically distinguish this informative eigenvalue from nearby uninformative ones.
It was argued informally by Angelini et al. [4] and recently proven by Stephan and Zhu [67] that, in the uniform case, eigenvalues other than α and β concentrate with high probability as n grows large in the disc of radius √ α. The informal argument generalizes smoothly. A generalized formal argument, alongside concentration results on the eigenvalues of B that strengthen Corollary 4.2, would be sufficient to prove the following conjecture: In the blockmodel η, with high probability as n → ∞, 1. The eigenvalue of largest magnitude of B approaches α.

2.
There exists a community-correlated eigenvector with an associated eigenvalue approaching β. 3. The bulk of the spectrum of B is confined to the disk of radius √ α. 4. NBHSC is able to recover labels correlated with ground-truth clusters iff These conjectures parallel known results for graphs [16] and uniform hypergraphs [67]. Conjecture 4.3 highlights a limitation of NBHSC in the setting in which edges of different sizes reflect cluster structure in different ways. Consider, for example, a setting with hyperedge sizes k and k ′ in which edges of size k tend to be within-cluster (β k > 0), while edges of size k ′ tend to be between-cluster (β k ′ < 0). In such a case, the approximate eigenvalue β = β k + β k ′ of B can be smaller in magnitude than either of β k and β k ′ , and could potentially be fully lost in the eigenvalue bulk, causing spectral clustering to fail. This can occur even when β k and β k ′ are sufficiently large in magnitude that clustering based on edges of size k or k ′ alone would succeed. This phenomenon is illustrated in Figure 1(b), in which an informative real eigenvalue lies hidden in the bulk of the spectrum and is hidden among other real, uninformative eigenvalues. The need to adaptively synthesize signals from hyperedges of differing sizes is one motivation of spectral methods based on the belief-propagation Jacobian, which we develop now.

Nonbacktracking
Operators and Belief Propagation. We will derive a connection between nonbacktracking operators and the belief-propagation algorithm (BP) [13] for community detection in sparse hypergraphs. Our derivation extends arguments by Krzakala et al. [43] in dyadic graphs and Angelini et al. [4] in uniform hypergraphs. We first sketch out the bird's-eye view of the argument. BP is known to be exact on graphical models which are trees, and often reliable on sparse models with local tree-like structure. Sparse stochastic blockmodels are an example of the latter case. Following standard arguments [43], we show that under certain additional symmetry assumptions, the approximate belief-propagation dynamics have a distinguished fixed point which contains no cluster information. Perturbations around this point are encoded in a Jacobian matrix-expressible in terms of the nonbacktracking operators B k -whose eigenvectors may therefore contain cluster information.

Belief-Propagation
Algorithm. We will work in the framework of the hypergraph stochastic blockmodel described in Section 4. In detection problems, we assume that we observe the edge set E but not the label vector z. We would like to obtain information about the conditional distribution η(z E) = η(E,z) η(E) of labels given the edge data E. An especially relevant summary is the marginal distribution of labels for each node, η(z i E).
We will estimate the marginals using belief-propagation [13]. We treat both subsets R ∈ R and labels z i ∈ Z as factors in a factor graph. There is a label factor for each node i. The message that each label factor sends to node i about its belief that z i = s is q (s) . The node subset factors are somewhat more complex. Let µ (s) iR denote the message that node i passes to the factor R expressing its belief that z i = s. Let ν (s) Ri denote the message that factor R passes to node i expressing its belief that z i = s. Then, the standard BP updates for this model read Here, Z iR and Z Ri are normalizing constants ensuring that ∑ s µ Here and below, sums over label vectors z should be assumed to run over Z subject to the explicitly specified constraints.
On factor graphs which are trees, the updates (5.1) and (5.2) converge, and the desired marginals can be obtained by computing the marginal message for node i: Our factor graph is admittedly not a tree. One possibility is to modify the belief propagation algorithm to account for loops [42]. We instead follow the standard argument that: • The sparsity assumption on η implies that the realized hypergraph is locally tree-like as n grows large. • The updates (5.1) and (5.2) can be approximated by updates that travel only along the edges of the realized hypergraph. Combined, these two points imply that BP can be approximated by an algorithm that takes place on a locally tree-like factor graph. The first point is discussed in the graph setting by [24]. The second has been argued heuristically in several papers [24,43,4]. 3 We provide below a heuristic statement (Claim 5.1) that expresses the thrust of this argument. We offer a rigorous statement with proof in Theorem B.1 (Appendix B).
Let F be the function with components This expression is obtained by substituting (5.2) into (5.1). We can then write the dynamics, restricted to the variable µ, as µ ← F(µ).
We are now prepared to relate the belief-propagation algorithm to the nonbacktracking matrices B k . A precise statement and proof of the relationship requires large amounts of additional notation so we defer them to Appendix B. Here, we state a heuristic version of the result.
k for all s, t ∈ Z and k ∈ K. Then, with high probability as n → ∞: • The pointμ with coordinatesμ (s) iR = q (s) for all i ∈ N , R ∈ R(i), and s ∈ Z is approximately a fixed point of F, in the sense that F(μ) ≐μ.
• The Jacobian of F atμ, restricted to a space of appropriately normalized perturbations, has entries of order O(n −1 ) except in a block corresponding to the set of realized edges. The restriction of the Jacobian to this block is J + O(n −1 ), where Claim 5.1 states that the uninformative pointμ is approximately a fixed point of the beliefpropagation update F. The stability of this fixed point is governed by the Jacobian matrix, which is dominated by a block which can be approximated by J. Because J approximates the BP dynamics in a neighborhood ofμ, and because the BP dynamics by design attempt to find community structure, we expect that the eigenvectors of J with large real eigenvalues to carry information about community structure in the hypergraph H.
m can be very large, and computing its eigenpairs can be costly. We therefore ask whether it is possible to compute on a smaller matrix. In the case of a k-uniform hypergraph, the relevant eigenvalues of J = G k ⊗ B k are of the form γβ, where γ is an eigenvalue of G k and β an eigenvalue of B k . This allows Angelini et al. [4] to compute only on B k , or indeed on the smaller B ′ k . In nonuniform hypergraphs, however, such a simple reduction is not available. We therefore offer a partial generalization in Theorem 5.2 that enables us to compute on a smaller matrix on nonuniform hypergraphs.
Theorem 5.2. There exists a matrix J ′ ∈ L(Z × K × N ) which can be expressed in terms of the parameter matrices {G k } k∈K and the adjacency matrices {A k } k∈K such that, if λu = Ju, then λx = J ′ x. In particular, either x = 0 or x is an eigenvector of J ′ with eigenvalue λ.
Writing J ′ and proving Theorem 5.2 explicitly requires the introduction of some cumbersome notation. We defer the detailed statement and proof to Theorem D.1 in Appendix D.
We focus on the eigenvectors u such that Lu ≠ 0 for clustering algorithms. In analogy to the use of x 2 in nonbacktracking hypergraph spectral clustering, we again use x 2 here. One intuition for this choice is that x (s) 2;ik sums beliefs that node i belongs to cluster s across edges of size k. The vector x 2 can be obtained either by computing eigenvectors of J and applying the transformation L or directly by computing eigenvectors of J ′ .

Alternating
Belief-Propagation Spectral Clustering. Claim 5.1 and Theorem 5.2 jointly suggest a modified algorithm based on the eigenvectors of J or J ′ rather than the eigenvectors of B or B ′ . Since J depends on the blockmodel parameters through the matrices {G k }, we alternate between spectral clustering steps and updates to these parameters. This alternating structure is reminiscent of expectation-maximization [27] and other coordinateascent algorithms. However, our alternating algorithm is not literally a form of coordinateascent because the spectral clustering step does not maximize a likelihood objective.
Indeed, one can carry out belief-propagation hypergraph spectral clustering without even specifying a likelihood objective. While the absence of a likelihood makes certain tasks harder, there is also an important computational benefit. In general, fully specifying a stochastic blockmodel requires specifying η(R ∈ E z R ) for every possible combination of labels z R . In a hypergraph with edges up to size k and group labels, there are k k! such combinations. Calculating a likelihood under the blockmodel therefore requires the estimation of a potentially very large number of parameters.
In contrast, belief-propagation hypergraph spectral clustering does not require specification of η(R ∈ E z R ), but only the entries of G k for each k. To do this, we need to estimate the label proportions q (s) and the pairwise edge counts m k (s, t) for each k, s, and t. This is a total of + 1 2 k 2 parameters, a number which scales much more favorably than the k k! parameters required for likelihood maximization. As a result, BPHSC is both less sensitive to the fine details of the HSBM parameters and more robust against overfitting concerns. We show how to estimate the necessary parameters given an estimated label vector in Appendix F.
The eigenvectors of J ′ are elements of V (2Z × K × N ). In order to carry out the clustering step, we need to obtain from these eigenvectors a set of feature vectors in V (N ) that carry information on the level of nodes. To do so, we use (5.3), which suggests that the linearized perturbations in the marginal label distributions can be obtained by summing over all messages incoming to node i. From the eigenvector x j we form the matrixỸ j ∈ L(N , Z) with entries y j;i,s = sign(∑ k∈K x (s) 2;ik ), with x (s) 2;ik as defined in (5.6). The extraction of the sign is again useful for noise reduction in small instances. 4 The columns of the matrices {Ỹ j } form the embedding coordinates on which we cluster. This clustering step gives a new label vector z. We use z to compute updated estimates of the model parameters, and repeat. A single stage of the clustering step is illustrated in Figure 2. It is also possible to carry out BPHSC using the full Jacobian matrix J rather than the reduced matrix J ′ . In this case, one computes eigenvectors {u } h =1 of J and forms from them the vectors {x } h =1 via (5.6). One then proceeds with the remainder of the algorithm. Whether it is best to use J or J ′ depends on the size of the hypergraph under consideration. The cost of BPHSC is dominated by (a) the cost of allocating the matrix used and (b) the cost of computing eigenpairs. Experimentally, we have found that the allocation cost of J ′ can exceed that of J on small instances. When performing many small experiments, such as in Figure 3, for example, we therefore use J directly. For larger instances, we found faster overall performance using J ′ .

Number of Eigenvectors.
In Algorithm 5.1, we extract from J ′ (or, alternatively, J) the h eigenvectors with real eigenvectors larger than unity in magnitude. This is very different from the standard approach in the case of sparse uniform hypergraphs and graphs. There, in a blockmodel with detectable communities, there are real eigenvalues outside the bulk of the spectrum with high probability. The real eigenvalue of largest magnitude is uninformative, and the remaining −1 real eigenvalues have eigenvectors from which the communities can be detected. In moderately nonsparse regimes, there are also real eigenvalues with community- Finally, we study the nodes in the space spanned by the columnsỹj;⋅,s for each aggregate eigenvectorỹj. In our example, this space has 8 × 3 = 24 dimensions. In a full run of Algorithm 5.1, we would then perform a clustering algorithm such as k-means to obtain labels. Here, we instead visualize the true clusters in two dimensions using t-SNE. Nodes in the same cluster appear close in the embedding space.
informative eigenvectors inside the bulk of the spectrum [22]; these eigenvalues can also be used in algorithms [23]. In the nonuniform case, the need to consider the Jacobian J (or J ′ ) rather than than the nonbacktracking matrix B complicates the situation considerably. The Jacobian J is a sum of terms (5.5), each of which contains a different size-specific nonbacktracking operator B k . The informative eigenvectors of an individual term are, unfortunately, not guaranteed to be even correlated with informative eigenvectors of the complete sum J. As a result, simply counting eigenvalues outside the bulk of the spectrum of J is not a reliable guide to the total number of detectable communities. Relating the number of detectable communities to the number of eigenvectors to extract from J is a nontrivial avenue of future work. Instead, we justify the choice of eigenvectors in Algorithm 5.1 from a dynamical perspective. Each eigenvalue of J larger than unity in magnitude corresponds to an unstable direction. Since messages correspond to probability distributions, only real perturbations to messages are interpretable. A reasonable heuristic choice is therefore to extract the eigenvectors of J or J ′ whose eigenvalues are both real and larger than unity in magnitude. Unlike in the case of uniform hypergraphs and graphs, informative real eigenvalues can lie inside the bulk of spectrum of J, as illustrated in Figure 2(a). The development of automated guidance for the location of informative eigenpairs in this setting in terms of the desired number of clusters and current parameter estimate would be a useful avenue of future study. The need to compute more eigenvectors on the relatively large matrices J or J ′ is an important practical limitation of our proposed algorithms.

Number of Clusters.
It is also necessary to specify the number of clusters in Algorithm 5.1, as this number is not determined by the number of informative eigenvectors. Many approaches to this problem take a model-selection perspective, choosing a number of clusters to optimize an information criterion or summarize a Bayesian posterior. Unfortunately, such approaches are not available here because BPHSC does not optimize a likelihood or other objective function.
In order to select the number of clusters and determine which clustering to finally accept, we therefore use a surrogate objective function. We use k-means for the Cluster step, and we use as an objective function the proportion of variance in the embedding space explained by the returned clusters. This enables direct comparisons between candidate clusterings with the same numbers of cluster labels, while scree plots can assist choices about the correct number of cluster labels to use.
6. Conjectured Thresholds for BPHSC. We now consider the performance of our spectral clustering algorithms NBHSC and BPHSC on sparse synthetic data generated by a simple hypergraph stochastic blockmodel. Our development is motivated in part by known behavior of spectral methods in sparse random graphs [57,43]. We first restate a standard definition of detection in clustering problems. For each n, let η n (θ, z n ) be a probability distribution over graphs on n nodes, parameterized by some vector θ. Each such distribution possesses the same shared set of parameters θ, as well as a planted partition z n of nodes. Let G n ∼ η n (θ, z n ). Let A be a clustering algorithm, which we view as a map G n ↦ẑ n from the data to an estimate of the true labels. Let ρ be a correlation function measuring label agreement with the property that if z n andẑ n are length-n labels sampled independently, then ρ(z,ẑ) → 0 with high probability as n → ∞. Examples of measures include mutual information, the Adjusted Rand Index, and suitably adjusted versions of the overlap [24,1]. Definition 6.1. Algorithm A detects communities in the sequence {η n , θ, z n } with respect to the correlation function ρ if there exists > 0 such that, with high probability as n grows large, ρ(z n , A(G n )) > .
As briefly described in Section 1, a notable setting is the graph stochastic blockmodel with equal-sized communities in which the parameter θ = (a, b) governs the rates of edges within and between communities. In this case, the size of the signal-to-noise ratio φ from (2.1) determines whether any algorithm exists that detects communities in {η n , θ, z n } as n → ∞. Definition 6.1 generalizes to the setting of hypergraphs with planted cluster structure. It is necessary only to allow η n (θ, z n ) to be a probability distribution over hypergraphs. Indeed, Angelini et al. [4] have offered conjectures concerning the ability of nonbacktracking spectral methods to detect planted clusters in the setting of uniform hypergraphs. Several of these conjectures were recently proven by Stephan and Zhu [67]. Our purpose in this section is to extend these conjectures to the setting of nonuniform hypergraphs. We will offer experimental support of these conjectures and leave their proofs to future work.
In BPHSC, it is necessary to form an estimate of the parameter matrices {G k } k∈K . In Algorithm 5.1, we form new estimates of these matrices in each iteration. While realistic, the need for re-estimation considerably complicates detectability analysis. We therefore consider an idealized setting in which the correct values of the parameter matrices {G k } k∈K are known a priori. The assumption that the parameters are known exactly is sometimes called the Nishimori condition in statistical physics [24]. 5 For each k, direct calculation shows that the vector 1 = (1, 1) is an eigenvector of G k with eigenvalue β k α k . Let u ∈ V ( ⃗ E) be as in Theorem 4.1. Letû = 1 ⊗ u andv = 1 ⊗ v. Let λ = ∑ k∈K β 2 k α k . Equation (5.5) and Theorem 4.1 now imply the following result: Corollary 6.2. Under the two-group blockmodel η, for = 1, 2, we have Importantly, it is not guaranteed that λ ≤ β. Unlike in the case of uniform hypergraphs, therefore, it may be the real eigenvalue of largest magnitude that carries community information. Figure 8 in Appendix G computes the two real eigenvalues of J of largest magnitude in a nonuniform hypergraph, finding excellent agreement with Corollary 6.2 and illustrating this phenomenon. Motivated by these results, we pose the following conjecture: Conjecture 6.3. In the same setting as Theorem 4.1, with high probability as n grows large, • J possesses a real, community-correlated eigenvector with eigenvalue λ + o(1).
• BPHSC, when initialized with knowledge of the true parameters c in k and c out k , is able to detect ground-truth clusters if λ > 1.
We now verify that BPHSC, initialized with knowledge of c in k and c out k , is able to detect clusters in broader parameter regimes than NBHSC. After algebraic rearrangement, we can write the respective detectability conditions from Conjectures 4.3 and 6.3 as Noting that α k ≥ 0 and β k = 0 if α k = 0, an application of Bergström's inequality [10,8] yields φ 2 ≥ φ 1 . Since φ 2 ≥ φ 1 , we conclude that, under Conjectures 4.3 and 6.3, BPHSC with true parameters succeeds in detecting communities in all cases for which NBHSC succeeds. We also expect to find cases in which φ 2 > 1 ≥ φ 1 . In such cases cases, our conjectures imply that BPHSC succeeds when NBHSC fails. We support this expectation experimentally in Figure 3 below.
6.1. Parameterized Thresholds. We now illustrate Conjectures 4.3 and 6.3 computationally. To do so, we first need to express conditions on the eigenvalues in terms of the parameters of an HSBM. We specify the affinity function implicitly via a ball-dropping process [61]. We first generate a number of k-edges. With probability p k , a k-edge is sampled uniformly at random from the set of all node k-subsets R such that all nodes in R have the same cluster label. With probability 1 − p k , the k-edge is sampled uniformly at random from the set of all k-subsets in which at least two nodes have differing labels. We refer to the former type of edge as within-cluster and the latter as between-cluster.
We first consider NBHSC. In this HSBM, (4.11) in Conjecture 4.3 defines a pair of hyperplanes in the coordinates {p k } k∈K . This follows from two facts. First, by construction, α does not depend on {p k }. Second, β is an affine function of {p k }: The proof is a direct calculation and provided in Appendix E. Lemma 6.4 in conjunction with Conjecture 4.3 define a pair of hyperplanes in the coordinates {p k }. The region between these hyperplanes is, under these conjectures, the region in which NBHSC fails to detect clusters. Figure 3(a). shows an experiment on a blockmodel on 400 nodes with 2-and 3-edges with varying p 2 and p 3 . In this panel, we run NBHSC 20 times for each parameter value, and compute the average Adjusted Rand Index (ARI) of the retrieved clustering against the planted clustering. White lines are the boundaries given by β 2 = α, with α and β given by Lemma 6.4. Under Conjecture 4.3, as n → ∞, spectral clustering returns labels with ARI bounded above 0 iff (p 2 , p 3 ) does not lie between the two white lines. The experimental results shown are consistent with this conjecture.
We now consider detectability thresholds for belief-propagation hypergraph spectral clustering under the Nishimori condition. We can compute the entries of the matrix G k in terms of c in k using (5.4), obtaining The eigenvalue of G k that appears in Corollary 6.2 is then Direct computation now shows that the condition λ = ∑ k∈K β 2 k α k = 1 defines an axis-aligned ellipsoid with coordinates (p k 1 , p k 2 , . . .), centroid (x k 1 , x k 2 , . . .) and radii (a k 1 , a k 2 , . . .), where (6.6) Conjecture 6.3 claims that BPHSC succeeds outside this ellipse (λ > 1) and fails inside (λ < 1). Panels (b)-(d). of Figure 3 show a sequence of cluster recovery experiments. As before, each pixel gives the average Adjusted Rand Index of the recovered cluster against ground truth across 20 runs of BPHSC in a hypergraph stochastic blockmodel on 400 nodes, with varying parameters {p k }. In each experiment, the ellipse defined by (6.6) is shown in white. An implication of our conjectures is that, as n → ∞, with high probability, belief-propagation spectral clustering returns a clustering with ARI approaching 0 iff the point {p k } lies in the interior of the ellipse. Careful examination shows that the algorithm occasionally succeeds within the ellipse, and occasionally fails outside it. We attribute these deviations from conjectured theory to finite-size effects. Recalling that we have treated the true matrices {G k } k∈K as known, the observed performance and thresholds should be regarded as idealizations of the more realistic case in which these matrices must be inferred along the way.
From Corollary 6.2 we know that there is also an approximate eigenpair (β,û) which is uncorrelated with planted cluster structure. The equation λ = β again describes an ellipsoid in parameter space at which the two eigenvalues collide. This collision induces noise in the corresponding eigenvectors, resulting in observably lower-quality cluster recovery along this ellipsoid (Figure 3(b-d), dashed black curves).
Our conjectured thresholds are asymptotic, while the synthetic hypergraphs on which we compute in Figure 3 are, at 400 nodes, relatively small. Even in this regime, the agreement with theory is quite strong. In Figure 9, we show an experiment in a smaller region of parameter space on hypergraphs of 10,000 nodes, again finding close agreement with conjectured theory.

Experiments on Data.
We first study several data sets in which ground-truth labels are available. The contact-primary-school [66,9] data set logs close-proximity human contact interactions detected by wearable sensors. Nodes are students or teachers. A hyperedge exists between a set of nodes that were jointly in proximity to each other within a short time window. Each student is assigned to a unique classroom, which we use as a ground-truth label. The data contain timestamps associated to each interaction, although we do not use these timestamps here. There are n = 242 nodes and m = 12, 704 hyperedges in contact-primary-school. Hyperedges range from size k = 2 to size k = 5. Figure 4 shows a suite of clustering experiments on contact-primary-school, which possesses 11 ground-truth clusters, including 10 homeroom classes and one cluster containing all teachers. We ran BPHSC multiple times, choosing from among runs the clustering that resulted in the smallest value of the k-means within-group sum-of-squares objective. We also varied the number of clusters to be learned. We repeated this process for both the original hypergraph data and the projected (clique-expansion) graph obtained by replacing each khyperedge with a k-clique. We refer to this algorithm as belief-propagation projected graph spectral clustering (BPPGSC). The lefthand plot shows that the k-means objective is able to give some guidance as to the appropriate number of groups to infer, with the objective function leveling off close to the true number of groups for both BPPGSC and BPHSC. At center, we observe that for most values of , including all those close to the correct value, BPHSC considerably outperforms dyadic spectral clustering in retrieving labels correlated with ground truth. One explanation for this phenomenon may be observed at right, where we plot the diagonal entries of the matrix C k for each edge size k, computed using the true cluster labels. These diagonal entries measure the rate at which nodes connect to other nodes in their same group, and may therefore be interpreted as a measure of affinity or assortativity. These affinities vary considerably according to edge size, suggesting that edges of differing sizes play meaningfully distinct roles in this data set. BPHSC again outperforms spectral clustering on the projected graph, with the gap likely due to the different connection structure across edges of varying sizes. While these results suggest that belief-propagation hypergraph spectral clustering is preferable to BPPGSC on these school contact data sets, Figure 4. Cluster recovery in the contact-primary-school data set [66,9]. We ran BPHSC on the data for 10 rounds, using the 30 eigenvectors of the belief-propagation Jacobian with largest real eigenvalues and with a varying number of clusters to be estimated. In each round, we update the estimate of the labelsẑ by choosing the best of 20 runs of k-means according to the within-group sum-of-squares objective. We repeat this experiment on the projected (clique-expansion) graph. (Left): scree plot of the mean within-group sum-ofsquares obtained by the k-means step as a function of the number of groups to be estimated. The vertical grey line gives the true number of labels in the data. (Center): Adjusted Rand Index of the clustering with lowest k-means objective against ground truth. (Right): The diagonal entries of the matrix C k for varying edge size k. A similar experiment for the contact-high-school data set [51,9] is given in Figure 10.  Figure 4, using the senate-bills data set [32,33,9].
neither algorithm achieves the perfect cluster recovery obtainable via other methods [20].
The senate-bills data set [32,33] provides a contrasting case. Nodes are U.S. senators. An edge exists on a set of senators for each bill that that set of senators cosponsored. The data reflects bills from the 103rd through 115th U.S. Congresses, which span the years 1993-2017. There are 293 senators and 20,006 bills represented. We consider bills cosponsored by between k = 2 and k = 25 senators, although a small number of bills sponsored by larger groups exist. Labels correspond to the two major U.S. political parties. Unlike in the contact networks, Figure 6. (Best viewed in color): Clustering the tags-math-sx data set [9] using belief-propagation spectral clustering on the projected graph (BPPGSC) and original hypergraph data (BPHSC). Only tags that appeared in at least 20 questions are included. In each panel, we performed 50 steps of Algorithm 5.1 using = 4 groups and h = 15 eigenvectors. We repeated this process 50 times, resulting in 50 candidate clusterings. For each algorithm, the clustering shown is the one which achieved the lowest k-means objective (total within-group sum-of-squares) of these 50 candidates. Colors give the learned clustering, while coordinates in the plane give the 2-dimensional principal components projection of the eigenvector embedding from which that clustering was learned. The 30 most frequently-used tags in the data set are labeled.
the scree plot entirely fails to track the number of true clusters. Both belief-propagation hypergraph spectral clustering and belief-propagation projected graph spectral clustering recover partitions with Adjusted Rand Index near 0.25 when instructed to search for exactly two groups, and do much more poorly otherwise. These results are qualitatively aligned with those of Chodrow et al. [20], who found that greedy modularity-based methods also struggle with the label recovery task on senate-bills. One explanation for why belief-propagation hypergraph spectral clustering does not achieve noticeable improvement over belief-propagation projected graph spectral clustering in this case is given by the plot of within-group affinities (far right). Here, the variability of these affinities with edge size is much smaller than it is for contact-primary-school, suggesting that the ability of BPHSC to distinguish roles for different edge sizes may not confer a useful advantage over BPPGSC on senate-bills.
It is also possible to use BPHSC to cluster data in which no labels are natively present. Figure 6 shows an experiment in which we use Algorithm 5.1 to cluster tags on the forum Math StackExchange [9]. Nodes are tags that usually express certain concepts or areas within mathematics. Edges correspond to questions asked on the forum. An edge exists between a set of nodes for each question asked to which the corresponding tags were applied. The data contains edges of size k = 2, . . . , 5, with an average total degree of 720 edges per node. There are a total of 1,419 nodes in the original data. We removed nodes with degree less than 20, resulting in a hypergraph containing 1,038 nodes. We show the result of belief-propagation spectral clustering on both the projected graph and the original hypergraph in Figure 6.
The large-scale output structure is relatively similar between the two algorithms. In both cases, for example, there is a clearly separated cluster (purple) associated to calculus , including the topics limits, infinity, and radicals. There is also a cluster (gold) associated with topics in discrete mathematics, including combinatorics, fibonacci-numbers, and elementary-set-theory. A third cluster (pink) focuses on linear algebra, with tags including orthogonality, matrices, and inner-product-space. A final cluster (cyan) includes a number of tags that do not fit neatly in a single mathematical subfield, such as definition, terminology, and proof-writing.
There are, however, several notable differences in the clusterings produced. The linear algebra cluster in panel (b) contains only six words each with clear connections to the field, while the cluster in panel (a) is much larger and includes several superficially unrelated tags, such as soft-question and algebra-precalculus. In addition, inner-product-space is separated in a different cluster. As another example, in panel (b) decimal-expansion is grouped with other topics in discrete mathematics, and is close in embedding space to number-systems. In contrast, decimal-expansion is far from any related topics in panel (a).
We note that cluster assessment is a challenging task and our approach is admittedly ad hoc. That said, we do find it clear that clusters found by NBHSC are preferable to those found by NBSC in groupings of mathematical topic tags.
8. Discussion. We have proposed and analyzed two spectral algorithms for clustering hypergraphs. We have focused on the distinctive challenges posed by nonuniform hypergraphs, which comprise many of the most interesting hypergraph data sets. To address the challenge of representing edges of multiple sizes in compatible data structures, we have employed the hypergraph nonbacktracking matrix B proposed by Storm [68]. We have considered both a simple spectral algorithm (NBHSC) that relies only on the eigenvectors of B, as well as a more complex spectral algorithm (BPHSC) based on the eigenvectors of the belief-propagation Jacobian J evaluated at the uninformative fixed point. In each case, we have provided theorems that allow us to compute on alternative matrices that are usually smaller, thereby enabling faster computation in some cases. We have demonstrated the latter algorithm on several hypergraph data sets, showing that it enjoys superior performance over nonbacktracking methods on projected graphs due to its ability to explicitly represent distinct statistical roles for edges depending on size. That said, we emphasize that belief-propagation hypergraph spectral clustering is not state-of-the-art for label recovery in hypergraph clustering tasks as measured by accuracy and especially by scalability. Even the use of Theorem 5.2 yields a matrix J ′ of size 2nk , which can become computationally challenging when either k or are large. In practice, the use of BPHSC is considerably constrained by both the memory requirements of forming J or J ′ and the performance of sparse eigensolvers [62].
Part of the significance of our approach is that it admits closed-form analysis of the parameter regime in which BPHSC fails to recover clusters correlated with ground truth. In the parameterization we use in Figure 3, this regime is ellipsoidal in coordinates {p k }. Our description of this region relies on several conjectures related to the location of the eigenvalues of the matrices {B k } and their associated eigenvectors. These conjectures are inspired by known results for the graph case [16]. Proofs of these conjectures would provide conclusive characterizations of the success and failure regimes for BPHSC. We expect such proofs to require tools from random matrix theory similar to those used in Stephan and Zhu [67].
In the dyadic stochastic blockmodel, it has been proven that, for the stochastic blockmodel on graphs with two clusters, the regime in which nonbacktracking spectral clustering fails coincides precisely with the regime in which the regime in which no algorithm can detect clusters [54,50]. We pose a similar conjecture for the hypergraph blockmodel: within the ellipsoid described by (6.6), the cluster detection problem cannot be solved by any algorithm. We propose a proof of this conjecture as a direction of future work.
Acknowledgements. Much of PSC's work on this paper was completed during his time in the Department of Mathematics of the University of California, Los Angeles. We are grateful to Mason Porter for donating computational resources to support this research.
Software and Data Availability. Code and data sufficient to reproduce experiments in this paper may be found at the GitHub repository jamiehadd/HypergraphSpectralClustering Primary computations were performed in the Julia language using a custom-written package [11]. Visualizations were constructed using ggplot2 for the R programming language [60,75].
Gender Representation in Cited Works. Recent work in several fields of science has identified gender bias in citation practices-papers by women and other gender-minoritized scientists are systematically under-cited in their fields [28,17,48,31,69,30,74,46].
In the spirit of Zurn et al.
[77], we performed an analysis of gender representation in the references cited in the main text of this manuscript. We manually gender-coded the first and last authors in the works cited according to personal acquaintance, instances of pronoun usage online, or first name. We focused on the first and last authors because typically, though not always, the former is the leading researcher and the latter the senior author in the disciplines included in our references. Our method of coding is limited in several ways. Gender is fundamentally nonbinary. Names and pronouns may not be indicative of gender. Gender may change over time. Manual coding is inherently flawed and subject to error. Furthermore, the heuristic that the first and last authors correspond to those which make the most important contributions to a manuscript is of varying validity in different areas of science, especially in mathematics.
Of the works cited in the main text (excluding this statement), we estimate that 18% had a non-male first author and 19% had a non-male last author. Of those with at least two authors, 25% had either a non-male first author or a non-male last author. Our proof approach extends Kempton's proof of the Ihara-Bass formula for graphs [41]. The same approach was used by Stephan and Zhu [67] for the case of uniform hypergraphs.
For each k, define operators S k ∈ L( ⃗ E, N ), T k ∈ L(N , ⃗ E), and W k ∈ L( ⃗ E) with entries: These operators satisfy several important relations. We begin with entrywise calculations:

Direct multiplication and use of (A.1)-(A.3) gives the relations
We are now prepared for the main computation. The push-through identity states that provided that X is invertible and all matrix products are well-defined. Kempton [41] provides an elementary proof. Using (A.6) to write I − µB = I − µST + µW and applying (A.7) to the righthand side gives det(I − µB) = det(I + µW) det(I − µT(I + µW) −1 S) .
The derivation of this inverse uses the fact that I + µW is a block-diagonal matrix with one block for each edge. Each block for an edge of size k has the form I + (µ − 1)E, where E is a k × k matrix of ones. The expressions p k (µ), q k (µ), and h k (µ) can be derived by assuming (I + (µ − 1)E) p k (µ)I + q ′ k (µ)E = I and solving. Then, q k (µ) = q ′ k (µ) − p k (µ).
Performing n row multiplications by h k (µ) for each k yields This gives the second factor in the statement of Theorem 3.3, so our final step is to address the factor det(I + µW). We have We find in turn det(I + µW) This completes the computation and the proof.
A.1. Proof of Corollary 3.4. Recall that ⃗ m is the total number of pointed edges. We make the substitution µ = 1 β in (3.1). From the result, we extract copies of the characteristic polynomials p B of B and p B ′ of B ′ . We obtain A.2. Proof of Lemma 3.5. Our proof closely follows that of Stephan and Zhu [67]. Define the matrices We then have Let x h;k ∈ V (N ) be the vector containing only the entries of x h corresponding to edges of size k for h = 1, 2. The calculation 10) shows that x 2;k =T k u. A similar calculation shows that x 1;k =T k W −1 u. We also make use of the following identities, which can be verified through calculations similar to those shown in Appendix A.
Definex h;k to be the vector with components [x h;k ] k ′ ,i = δ k,k ′ x h;k,i for h = 1, 2. Then, the relationT k W −1 u implies thatT k W −1ū k =x 1;k . On the other hand, we computē We have used the fact that, sinceū k is nonzero only in the entries in whichW k is, The identity (A.13) implies that W −1 S = diag(T 1 , . . . , Tk). The identity (A.14) then gives We have ∑ k∈KDk = D. We have thus far shown that Summing over k gives which establishes our first needed relation.
Let us now instead premultiply the relation βū k = Bū k byT k . The lefthand side becomes βx 2;k . The righthand side becomes We haveT k S =Ā k , whereĀ k is the matrix Using this and (A.11), we obtain Summing over edge sizes k gives Let Z be the alphabet of cluster labels, with Z = . Let M be the space of possible messages µ; we can identify M with a product of probability ( − 1)-simplices containing one factor for each node-tuple pair. Letμ be the vector of messages with entriesμ (s) iR = q (s) for all nodes i, subsets R ∈ R(i), and labels s ∈ Z.
We will consider perturbations to the belief-propagation dynamics. The normalization condition ∑ s∈Z µ (s) iR = 1 on elements of M requires that perturbations to a message vector µ must satisfy ∑ s∈Z (s) iR = 0. Letting Π denote the projection operator onto the subspace defined by this relation, we have For a given hypergraph realization, we can separate µ into components (µ 0 , µ 1 ), where entries of µ 0 correspond to unrealized edges and entries of µ 1 correspond to realized edges. We can similarly separate the components of the function F. This allows us to write the BP update dynamics in the form We are now prepared to state a precise analog to to the heuristic Claim 5.1.
Theorem B.1. Let H be sampled from the sparse Bernoulli blockmodel. Then, as n grows large, F(μ) ≐μ. Furthermore, Proof. Throughout this proof, sums indexed by z R are assumed to run through Z R subject to stated constraints.
We'll first compute several approximations describing how messages propagate along unrealized edges, i.e. subsets R such that a R = 0. Since η(a R = 0 z R ) = 1 − ω(z R )n 1− R , (5.2) becomes The last line follows from the normalization of the messages µ (s) jR , since the sum ranges over all possible labelings of the nodes in R ∖ i. Since the smallest possible edge size is k = 2, we have shown that ν (s) Ri is approximately constant with respect to s. We can also approximate µ Here, we are able to identify the normalizing constant Z i = Z iR Z Ri independent of R because it normalizes an expression independent of R.
Let's now consider how messages are passed along edges R such that a R = 1. This corresponds to the consideration of F 1 . Substituting (5.2) into (5.1) and absorbing normalizing constants allows us to eliminate the messages ν (s) Ri entirely, obtaining an explicit form for F 1 : The updates in µ 1 are now µ 1 ← F 1 (µ 0 , µ 1 ).
Let us write F 1 in the form where M where we have used the normalization of the messages µ (z j ) jQ in the last line. We will now approximate N (s) iR by a "field term" N (s) which does not depend on i or R. First, since a Q = 0 for each Q appearing in the product defining N (s) iR , we approximate µ (z j ) The number of factors in the denominator is equal to the degree of node i, which is binomial and therefore concentrates about its mean c (s) k . We therefore have that, with high probability as n grows large, the entire denominator is then also 1 − O(n −1 ). With high probability, then, to be a constant "field term" which does not depend on i or R. Thus, with high probability as n grows large, the message passing update (B.5) satisfies Here we have absorbed factors that do not depend on s into Z iR . Importantly, F 1 depends on µ 0 only through the field term N (s) .
We next consider the behavior of (B.7) near the pointμ with entriesμ (s) iR = q (s) . Let be a perturbation vector assumed small. We'll first consider the field term. Let ∆ refer to terms of order O(n −1 + ). The perturbed field term N (s) now reads where we have defined N to absorb the products. We have shown that, nearμ, the field term N (s) approximately does not depend on s. Paralleling the partition µ = (µ 0 , µ 1 ), we can partition the entries of as = ( 0 , 1 ), again corresponding to unrealized and realized edges. From (B.7), By conditioning on the label of j, we can simplify the second term in the factor: Since the messages must normalize in s, we have that, up to errors that can be absorbed into the first factor, N Z iR = 1. This shows that A similar calculation shows that F 0 (μ 0 ,μ 1 ) ≐μ 0 , which relations jointly give F(μ) ≐μ. This proves the first clause of the theorem.
Furthermore, since F 1 depends onμ 0 only through the factor N Z iR = 1 + ∆ = 1 + O(n −1 + ), we have that any derivative of F 1 in a direction corresponding to µ 0 is of order O(n −1 ). Thus, the Jacobian ∂F 1 ∂µ 0 , evaluated atμ, has entries of order O(n −1 ). The projected matrix is also of order O(n −1 ), proving the first equation in (B.2). It remains to compute J 11 atμ. Expanding the product in (B.8) to first order in and separating the arguments of F 1 gives Using (B.9) gives We now apply the projection Π onto the subspace of admissible perturbations, yielding We can identify R with an edge e, and the pair (i, R) as a pointed edge ⃗ e with p(⃗ e) = i. Doing the same for Q and j, we can recognize the two rightmost sums as the action of B k on the perturbation vector . Using the definition of G k , we can write this relation as Ignoring the error term, this relation would define the Jacobian J 11 ∶= Π ∂F 1 ∂µ 1 as equal to the righthand side. Allowing → 0, we conclude that J 11 satisfies which establishes the second clause of (B.2) and completes the proof.
Appendix C. Proof of Theorem 4.1.
We will prove (4.7); the proof of (4.6) is similar but somewhat shorter. Letṽ be a vector indexed by tuples Q ∈ [n] k and nodes i ∈ Q with entries otherwise.
We are going to show that E (B kṽ ) iQ − β kṽiQ Q ∈ E ≐ 0; sinceB k andṽ agree with B k and v conditioned on the event Q ∈ E, this will imply (4.7). We proceed via direct computation. Expanding the expectation, we can write The third line follows from the fact that, conditioned on the event Q, R ∈ E,b iQ,jR = b iQ,jR andṽ jR = v jR . Proceeding from the fourth line, we can evaluate b iQ,jR and rearrange the sums: For the third line, we have used the fact that there are approximately n 2 terms in each sum. This completes the proof.
In the second line we have used the mixed product property. The third follows from direct multiplication, finding that the products involving mixed edge sizes zero out. The same argument as in Appendix A.1 shows thatT k W −1 B kū (s) k = D k − I κn x (s) 2;k , wherex (s) 2;k ∈ V (N ) has entries equal to x 2 in edge size k and zero otherwise, with all entries corresponding to group s. We thus have The second line is again the mixed product property, while for the third we have definedx 2;k to have entries that agree with x 2 on edges of size k and which are zero otherwise. Summing over k yields, after some algebra, our first reduced relation: where the matrices G, H, andD are defined in the statement of Theorem D.1.
We now premultiply both sides of the eigenvector relation ξu = Ju by the matrix I ⊗ ∑ k∈KTk . The lefthand side becomes ξx 2 . For the righthand side, we compute The steps of the calculation so far precisely parallel the calculation of (D.2). Definingx Performing the sum over k yields, after some further algebra, Defining the matrix inside the braces as J ′ completes the proof.
where r k = 1−2 2−k 2−2 2−k . Combining this with our previous results, we have In turn, we have  On a small technical note,m s,t k should be computed counting multiplicities; for example, a 5edge with labels (a, a, b, b, c) would make four contributions tom k (a, b) and two contributions to bothm a,c k andm b,c k . Appendix G. Additional Experiments. Figure 7 supports the use of sign-based thresholding in Algorithm 5.1, finding overall improved recovery as measured by the Adjusted Rand Index when using thresholding. Figure 8 offers support of our conjectures for the locations of informative eigenvalues in binary detection experiments for both of the matrices B and J. This figure also illustrates that the informative eigenvalue for the matrix J may be the largest real eigenvalue in magnitude, rather than the second-largest as is true for B. Figure 9 supports the accuracy of the thresholds predicted by Conjecture 6.3 in a much larger synthetic hypergraph of 10, 000 nodes. The parameter space explored corresponds to a vertical slice of Figure 3(a-b) with p 2 = 0.35.
We also include Figure 10 below, which replicates Figure 4 for the contact-high-school data set. Figure 8. Locations of the two real eigenvalues with largest magnitude for the nonbacktracking matrix B and the belief-propagation Jacobian matrix J (with true parameters) for 100 realizations of a random hypergraph as described in Section 6.1. The parameters are p2 = 0.2, p3 = 0.9, and c2 = c3 = 5, with 200 nodes in each of the two groups. The vertical axis gives the Adjusted Rand Index of the corresponding signs of the eigenvector against the planted community labels. Ticks on the horizontal axis give the predictions according to Corollary 4.2 and Corollary 6.2. The predictions for the second eigenvalues of these two matrices overlap exactly. The real eigenvalue of B with largest magnitude is uncorrelated with community structure, while the eigenvalue of second-largest magnitude is correlated. In contrast, the situation is reversed for J. Figure 9. Binary detection experiment similar to those in Figure 3, with p2 = 0.35 and varying p3. In this experiment, c2 = c3 = 5 and c4 = 0, with n = 10, 000 nodes. We show performance as measured by ARI for both BPHSC and NBHSC. Dotted lines give detectability thresholds for each algorithm (the two lower thresholds nearly overlap) from Conjecture 6.3. The dashed line gives the estimated location collision of the two eigenvalues described by Corollary 6.2. Figure 10. Experiment on the contact-high-school data set [51,9], analogous to the experiment shown in Figure 4. There are n = 327 nodes and m = 7, 818 hyperedges. We ran BPHSC on the data for 10 rounds, using the 30 eigenvectors of the belief-propagation Jacobian with real eigenvalues of greatest magnitude and with a varying number of clusters to be estimated. In each round, we update the estimate of the labelsẑ by choosing the best of 20 runs of k-means according to the within-group sum-of-squares objective. We repeat this experiment on the projected (clique-expansion) graph. (Left): scree plot of the mean within-group sum-of-squares obtained by the k-means step as a function of the number of groups to be estimated. The vertical grey line gives the true number of labels in the data. (Center): Adjusted Rand Index of the clustering with lowest k-means objective against ground truth. (Right): The diagonal entries of the matrix C k for varying edge size k.