remarkRemark hypothesisHypothesis claimClaim Hypergraph Spectral ClusteringP.
S. Chodrow, N. Eikmeier, and J. L. Haddock conjConjecture corCorollary
supplement

# Nonbacktracking Spectral Clustering of Nonuniform Hypergraphs

Philip S. Chodrow Department of Mathematics, University of California, Los
Angeles<EMAIL_ADDRESS>Nicole Eikmeier Department of Computer Science,
Grinnell College, Grinnell, IA Jamie Haddock Department of Mathematics, Harvey
Mudd College, Claremont, CAJH is partially supported by NSF DMS #2211318

###### Abstract

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 to enable faster computation
of eigenpairs. We then propose an alternating algorithm for inference in a
hypergraph stochastic blockmodel via linearized belief-propagation, offering
proofs that both formalize and extend several previous results. We perform
experiments in real and synthetic data that underscore the benefits of
hypergraph methods over graph-based ones when interactions of different sizes
carry different information about cluster structure. Through an analysis of
our algorithm, we pose several conjectures about the limits of spectral
methods and detectability in hypergraph stochastic blockmodels writ large.

###### keywords:

networks, hypergraphs, spectral methods, community detection, nonbacktracking
matrix

05C50, 05C65, 15A18, 62H30, 62R07, 91D30

##  1 Introduction

Graphs provide a classical representation for systems with pairwise
relationships: components are modeled by nodes, and relationships are modeled
by edges. In many systems however, relationships may simultaneously involve
more than two components. Examples include three scholars writing a paper
together, or three or more genes interacting to influence phenotypic traits.
While certain graph techniques can be used in such cases, there is often value
in keeping multiway interactions intact [2, 13]. In such cases, _polyadic_ or
_higher-order_ representations are useful, and analysis based on such
representations can lead to qualitatively and qualitatively different
conclusions [17]. There has been a wealth of recent work using hypergraphs,
simplicial complexes, and related structures for modeling complex polyadic
systems [10, 6, 7].

In this paper, we focus on the _community detection_ problem for hypergraphs.
The community detection problem asks us to partition the nodes of a hypergraph
into useful or insightful subsets, which are often called _communities_ ,
_clusters_ , or simply _groups_. Many algorithms exist for dyadic graphs [54].
There is also a growing body of algorithms for hypergraph community detection;
Chodrow et al. [18] give a brief survey of extant approaches and applications.

In the context of graphs, _spectral methods_ are a well-studied family of
clustering algorithms. A spectral method for clustering a graph $\mathcal{G}$
proceeds by computing a distinguished set of eigenvectors of some matrix
$\mathbf{M}$ associated to $\mathcal{G}$. A point-cloud clustering algorithm
may then be used to obtain clusters from the embedding space defined by these
eigenvectors. There are several possibilities for the choice of matrix
$\mathbf{M}$, including the graph adjacency matrix [51], various Laplacian
matrices [59, 66], the modularity matrix [52], and the nonbacktracking matrix
[37]. Unlike greedy methods such as the famous Louvain algorithm [12], many
spectral methods can be interpreted as approximate global solutions to cluster
optimization problems. Many spectral algorithms admit statistical guarantees
for the recovery of planted clusters in the large-sample limit [67, 38].
Additionally, certain methods admit asymptotic analysis of the conditions
under which they _fail_ to detect planted clusters [51, 37]; in some cases,
these conditions coincide with detectability thresholds that constrain the
performance of any algorithm [47, 48, 44]. The analysis of spectral methods
can thus offer insight on the conditions under which cluster detection is
possible. In practice, spectral methods are constrained by the performance of
sparse eigensolvers [57] and by general theoretical limitations inherited by
any clustering algorithm [53].

In this paper, we study spectral methods based on the hypergraph
nonbacktracking operator [62]. Several spectral clustering algorithms exist
for _uniform_ hypergraphs—hypergraphs in which all edges contain the same
number of nodes—including one based on the nonbacktracking operator [4, 61].
Many interesting hypergraph data sets, however, are nonuniform. We focus here
on the extension, analysis, and application of nonbacktracking methods to
nonuniform hypergraph data.

The paper proceeds as follows. In Section 2, we survey related work on the
dyadic nonbacktracking operator and its uses in graph data science. We discuss
spectral clustering techniques applied to hypergraphs, and discuss the
hypergraph nonbacktracking operator $\mathbf{B}$ introduced by Storm [62]. 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 [37, 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 between the Jacobian governing the
stability of this fixed point and a smaller matrix, 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 nonbacktracking spectral clustering
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.

###  2.1 Preliminaries and Notation

A hypergraph $\mathcal{H}$ is a tuple $(\mathcal{N},\mathcal{E})$ consisting
of a node set $\mathcal{N}$ and an edge set $\mathcal{E}$. The node set
$\mathcal{N}$ contains $n\coloneqq\lvert\mathcal{N}\rvert$ nodes. The edge set
$\mathcal{E}$ consists of subsets of $\mathcal{N}$. For each $k$ in a set $K$
of possible edge sizes, we let $\mathcal{E}_{k}$ denote the set of edges of
size $k$. We let $m_{k}=\lvert\mathcal{E}_{k}\rvert$, and $m=\sum_{k\in
K}m_{k}$. We usually assume $K=\\{2,3,\ldots,\bar{k}\\}$ for some maximum edge
size $\bar{k}$, and set $\kappa\coloneqq\lvert K\rvert$. 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 $\mathcal{G}$ instead of $\mathcal{H}$. We let $\langle
k\rangle=\frac{1}{m}\sum_{k\in 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 $\sum_{s\in S}sa_{s}$ for scalars $a_{s}$.
For integer $j$, we define the notation $V(jS)\coloneqq V(S)\oplus
V(S)\cdots\oplus 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)\coloneqq L(S,S)$. We speak of elements of $L(S,T)$
interchangeably as either linear maps or matrices with respect to $S$ and $T$
regarded as bases of $V(S)$ and $V(T)$, respectively. So defined, $V(S)$ is
isomorphic to $\mathbb{R}^{\lvert S\rvert}$ and $L(S,T)$ is isomorphic to
$\mathbb{R}^{\lvert S\rvert\times\lvert T\rvert}$. 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: $\mathbf{v}\in V(S)$. An
entry of $\mathbf{v}$ is $v_{s}$ for some $s\in S$. Elements of $L(S,T)$ are
written in uppercase bold: $\mathbf{A}\in L(S,T)$. An entry of $\mathbf{A}$ is
$a_{s,t}$ for $s\in S$ and $t\in T$. In many cases we will need to consider
multiply-indexed structures; for example, an element $\mathbf{u}\in V(S\times
T)$ has elements of the form $v_{st}$ for $s\in S$ and $t\in T$. An element
$\mathbf{C}\in L(S\times T,S^{\prime}\times T^{\prime})$ has elements of the
form $c_{st,s^{\prime}t^{\prime}}$. When considering an indexed family of
matrices such as $\\{\mathbf{A}_{1},\ldots,\mathbf{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 $\mathbf{A}_{k}$. We use $\mathbbm{1}[P]$ to denote
the indicator function of the proposition $P$, and the shorthand
$\delta_{i,j}\coloneqq\mathbbm{1}[i=j]$.

We now define the nonbacktracking operator on hypergraphs as given by Storm
[62].

######  Definition 2.1 (Pointed Line Graph).

A _pointed edge_ $iQ$ in a hypergraph $\mathcal{H}$ consists of an edge
$Q\in\mathcal{E}$ along with a choice of _point_ $i\in Q$. Define
$\vec{\mathcal{E}}$ to be the set of pointed edges. Let $\vec{\mathcal{E}}(i)$
be the set of pointed edges with specified point $i$. We say that pointed edge
$jR$ _follows_ pointed edge $iQ$, written $iQ\rightarrow jR$, if $j\in
Q\setminus i$ and $Q\neq R$. The _pointed line graph_ $\mathcal{P}$ of
$\mathcal{H}$ is a directed graph whose nodes are elements of
$\vec{\mathcal{E}}$. There is a directed edge $(iQ,jR)$ in $\mathcal{P}$ if
$iQ\rightarrow jR$.

######  Definition 2.2 (Nonbacktracking Operator).

The nonbacktracking operator $\mathbf{B}\in L(\vec{\mathcal{E}})$ associated
to a hypergraph $\mathcal{H}$ is the directed adjacency operator of the
pointed line graph $\mathcal{P}$. Its entries are
$b_{iQ,jR}\coloneqq\mathbbm{1}[iQ\rightarrow jR]$.

####  2.1.1 Nonbacktracking Methods in Network Data Science

The nonbacktracking matrix $\mathbf{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 $\mathbf{B}$ mix
more rapidly than random walks governed by the adjacency matrix $\mathbf{A}$,
implying that these walks may be more efficient for graph exploration tasks.
Martin et al. [43] propose an eigenvector centrality measure based on
$\mathbf{B}$, and show that this centrality avoids the pathological
localization of classical adjacency-based eigenvector centrality in sparse
graphs. Torres et al. [64] and Mellor and Grusovin [46] 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. [65] develop perturbation theory for
the leading eigenvalues of $\mathbf{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 second
eigenvector of $\mathbf{B}$ can be used to assign labels to nodes correlated
with the true clusters [37]. We refer to this algorithm and variations as
_nonbacktracking spectral clustering_ (NBSC). The behavior of NBSC under
random graph models is governed by the behavior of the leading eigenpairs of
$\mathbf{B}$, for which concentration results are recently available [14].
Unlike spectral methods based on the adjacency or Laplacian matrices, NBSC
succeeds in detecting communities all the way down to the _detectability
threshold_ , which constrains the ability of any algorithm to return clusters
correlated with the true partition [37]. This threshold was first conjectured
by Decelle et al. [20] using insights from statistical physics. It was later
proved in a series of subsequent papers [47, 48, 44]. Angelini et al. [4]
posed a conjectured detectability threshold for the clustering problem in
uniform hypergraphs.

###  2.2 Spectral Clustering Methods for Hypergraphs

There is a wide range of algorithms for clustering and community detection in
hypergraphs; see Chodrow et al. [18] 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. Perhaps the most direct approach is to
replace 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
$\binom{k}{2}$, possibly-weighted 2-edges between each pair of nodes contained
in $e$. Graph spectral methods can then be applied [70]. Multiple versions of
this approach possess asymptotic consistency guarantees under hypergraph
planted partition models [29, 23]. 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 $k$-tensor
adjacency tensor $\mathbf{A}$. In the $k=3$ case, for example, we would have
entry $a_{i,j,k}=1$ iff $\left\\{i,j,k\right\\}\in\mathcal{E}$. There are
several spectral methods for clustering $k$-uniform hypergraphs, many of which
rely on concepts connected to tensor eigenvalues and eigenvectors [39]. Ke et
al. [34] use a normalized tensor power iteration to compute eigenvectors,
while Chang et al. [16] take an explicit optimization approach. Hu and Wang
[30] 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 represented in the same tensor. A proposal was
recently made in this direction by Galuppi et al. [28], 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 [28], random walks [50], and Laplace operators
[49, 32]. To our knowledge, the application of these techniques in the data
scientific context has not yet been studied.

The use of $\mathbf{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
hypergraph planted partition model. Several of the conjectures for uniform
hypergraphs were recently formalized and proved by Stephan and Zhu [61].

##  3 An Ihara-Bass Theorem for Nonuniform Hypergraphs

A theorem often attributed to both Ihara [31] and Bass [5] relates the
spectrum of $\mathbf{B}$ to the spectrum of a smaller matrix in the case of
graphs. We now extend this theorem to nonuniform hypergraphs. Our result
generalizes a computation by Angelini et al. [4] and theorem by Stephan and
Zhu [61].

It is useful to distinguish blocks of $\mathbf{B}$ by edge size. Let
$\vec{\mathcal{E}}_{k}$ be the set of pointed edges of size $k$, let
$\vec{m}_{k}=\lvert\vec{\mathcal{E}}_{k}\rvert$, and let $\vec{m}=\sum_{k\in
K}\vec{m}_{k}$.

######  Definition 3.1 (Size-Restricted Nonbacktracking Operators).

Let $\mathbf{B}_{k^{\prime}\rightarrow k}\in
L(\vec{\mathcal{E}}_{k},\vec{\mathcal{E}}_{k^{\prime}})$ have entries
$b_{k^{\prime}\rightarrow k;iQ,jR}\coloneqq\mathbbm{1}[iQ\rightarrow jR]\;,$
where $iQ\in\vec{\mathcal{E}}_{k^{\prime}}$ and $jR\in\vec{\mathcal{E}}_{k}$.
The _$k$ -th nonbacktracking operator_ $\mathbf{B}_{k}\in
L(\vec{\mathcal{E}})$ associated to hypergraph $\mathcal{H}$ has the block
form

$\displaystyle\mathbf{B}_{k}=\left[\begin{matrix}\mathbf{0}_{2m_{2}\times
2m_{2}}&\mathbf{0}_{2m_{2}\times 3m_{3}}&\cdots&\mathbf{B}_{2\rightarrow
k}&\cdots&\mathbf{0}_{2m_{2}\times\bar{k}m_{\bar{k}}}\\\
\mathbf{0}_{3m_{3}\times 2m_{2}}&\mathbf{0}_{3m_{3}\times
3m_{3}}&\cdots&\mathbf{B}_{3\rightarrow
k}&\cdots&\mathbf{0}_{3m_{3}\times\bar{k}m_{\bar{k}}}\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\mathbf{0}_{\bar{k}m_{\bar{k}}\times
2m_{2}}&\mathbf{0}_{\bar{k}m_{\bar{k}}\times
3m_{3}}&\cdots&\mathbf{B}_{\bar{k}\rightarrow
k}&\cdots&\mathbf{0}_{\bar{k}m_{\bar{k}}\times\bar{k}m_{\bar{k}}}\end{matrix}\right]\;,$

where $\mathbf{0}_{m_{1}\times m_{2}}$ is the matrix of zeros of specified
dimensions.

By construction, $\mathbf{B}=\sum_{k\in K}\mathbf{B}_{k}$. Because the only
indices in $\mathbf{B}_{k}$ permitted to be nonzero in both rows and columns
are the indices corresponding to $\mathbf{B}_{k\rightarrow k}$, any eigenvalue
of $\mathbf{B}_{k}$ must either be zero or an eigenvalue of
$\mathbf{B}_{k\rightarrow k}$.

Table 1: Table of selected major notation used. Symbol |  Meaning | Introduced In
---|---|---
$\mathcal{H}=(\mathcal{N},\mathcal{E})$ |  a hypergraph with nodes $\mathcal{N}$ and edges $\mathcal{E}$ | Section 2.1
$n=\lvert\mathcal{N}\rvert$ |  the number of nodes | Section 2.1
$K$, $\kappa$ |  set of possible edge sizes of a hypergraph, $\kappa=\lvert K\rvert$ | Section 2.1
$\mathcal{E}_{k}$, $m_{k}$ |  the set of edges of size $k$, $m_{k}=\lvert\mathcal{E}_{k}\rvert$ | Definition 2.2
$iQ$ |  a pointed edge, $i\in\mathcal{N}$ and $Q\in\mathcal{E}$ | Definition 2.1
$\vec{\mathcal{E}}$, $\vec{m}$ |  set of pointed edges of $\mathcal{H}$, $\vec{m}=\lvert\vec{\mathcal{E}}\rvert$ | Definition 2.2
$\vec{\mathcal{E}_{k}}$, $\vec{m}_{k}$ |  set of pointed edges of size $k$, $\vec{m}_{k}=\lvert\vec{\mathcal{E}}_{k}\rvert$ | Definition 2.2
$\vec{\mathcal{E}}_{k}(i)$ |  set of pointed edges with point $i$ | Definition 2.1
$\mathcal{P}$ |  pointed line graph of a hypergraph | Definition 2.2
$V(S)$ |  finite vector space with orthonormal basis indexed by elements of set $S$ | Section 2.1
$L(S,T,$) |  Space of linear maps $V(S)\rightarrow V(T)$ | Section 2.1
$\mathbf{B}_{k\rightarrow k^{\prime}}$ |  nonbacktracking operator from $k$-edges to $k^{\prime}$-edges | Definition 3.1
$\mathbf{B}_{k}$ |  nonbacktracking operator from $k$-edges to all edges | Definition 3.1
$\mathbf{B}$ |  hypergraph nonbacktracking operator | Definition 3.1
$\mathbf{A}_{k},\mathbf{D}_{k}$ |  $k$-th adjacency and degree operators for a hypergraph | Definition 3.2
$\mathbf{K}$ |  square matrix with diagonal entries in $K$ | Section 3
$\mathcal{R}$, $\mathcal{R}_{k}$ |  set of all subsets or $k$-subsets of nodes | Section 4
$\mathcal{R}_{k}(i)$, $\mathcal{R}_{k}(i)$ |  set of all subsets of nodes containing node $i$ | Section 4
$x\doteq y$ |  $x=(1+O(n^{-r}))y$ w.h.p. for some $r>0$. | Section 4
$c_{k}^{(s)}$ |  Mean $k$-degree of nodes in cluster $s$. | Equation 6
$\mathbf{J}$ |  Non-vanishing block of the Jacobian matrix of belief-propagation evaluated at the uninformative fixed point. | 2
NBHSC |  Nonbacktracking Hypergraph Spectral Clustering | Algorithm 1
BPHSC |  Belief-Propagation Hypergraph Spectral Clustering | Algorithm 2

We now define adjacency and degree operators, which we also distinguish by
edge size.

######  Definition 3.2 (Adjacency and Degree Operators).

The _$k$ -th adjacency operator_ $\mathbf{A}_{k}\in L(\mathcal{N})$ associated
to $\mathcal{H}$ has entries

$a_{k;i,j}\coloneqq\sum_{R\in\mathcal{E}_{k}}\mathbbm{1}[\left\\{i,j\right\\}\subseteq
R]\;.$

The _$k$ -th degree matrix_ $\mathbf{D}_{k}\in L(\mathcal{N})$ associated to
$\mathcal{H}$ has entries

$d_{k;i,j}\coloneqq\delta_{i,j}\sum_{R\in\mathcal{E}_{k}}\mathbbm{1}[i\in
R]\;.$

Define the block matrices

$\displaystyle\mathbf{A}\coloneqq\left[\begin{matrix}\mathbf{A}_{2}&\cdots&\mathbf{A}_{\bar{k}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{A}_{2}&\cdots&\mathbf{A}_{\bar{k}}\end{matrix}\right]\in
L(K\times\mathcal{N})\quad\text{and}\quad\mathbf{D}\coloneqq\left[\begin{matrix}\mathbf{D}_{2}&\cdots&\mathbf{D}_{\bar{k}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{D}_{2}&\cdots&\mathbf{D}_{\bar{k}}\end{matrix}\right]\in
L(K\times\mathcal{N})\;,$

where $\mathbf{A}_{k}$ and $\mathbf{D}_{k}$ are as in Definition 3.2. Let
$\mathbf{K}$ be the square matrix with entries $(2,3,\ldots,\bar{k})$ along
the diagonal, and zeroes everywhere else.

######  Theorem 3.3 (Ihara-Bass for nonuniform hypergraphs).

For any hypergraph $\mathcal{H}$, we have

$\det(\mathbf{I}-\mu\mathbf{B})=f_{\mathcal{H}}(\mu)\det\left(\mathbf{I}_{\kappa
n}+\mu((\mathbf{K}-2\mathbf{I}_{\kappa})\otimes\mathbf{I}_{n}-\mathbf{A})+\mu^{2}((\mathbf{K}-\mathbf{I}_{\kappa})\otimes\mathbf{I}_{n})(\mathbf{D}-\mathbf{I}_{\kappa
n})\right)\;,$

where

$f_{\mathcal{H}}(\mu)=\prod_{k\in
K}(1-\mu)^{m_{k}(k-1)-n}(1+\mu(k-1))^{m_{k}-n}\;,$

and $\otimes$ is the Kronecker product.

We supply a proof of Theorem 3.3, as well as of 1 and 3.4 below, in Appendix
A.

The computational significance of Theorem 3.3 is that we can compute the
interesting eigenvalues of $\mathbf{B}$ via a smaller matrix
$\mathbf{B}^{\prime}$.

######  Corollary 1.

Let $\mathcal{H}$ be a hypergraph with $m_{k}$ edges of size $k$ for each
$k\in K$. Then:

  * •

For each $k$, if $m_{k}>n$, then $\beta=1-k$ is an eigenvalue of $\mathbf{B}$
with algebraic multiplicity at least $m_{k}-n$.

  * •

If $\sum_{k\in K}m_{k}(k-1)>\kappa n$, then $\beta=1$ is an eigenvalue of
$\mathbf{B}$ with algebraic multiplicity at least $\sum_{k\in
K}m_{k}(k-1)-\kappa n$.

  * •

The remaining eigenvalues of $\mathbf{B}$ are eigenvalues of the matrix

(1)
$\displaystyle\mathbf{B}^{\prime}\coloneqq\left[\begin{matrix}\mathbf{0}_{\kappa
n}&\mathbf{D}-\mathbf{I}_{\kappa n}\\\
(\mathbf{I}_{\kappa}-\mathbf{K})\otimes\mathbf{I}_{n}&\mathbf{A}+(2\mathbf{I}_{\kappa}-\mathbf{K})\otimes\mathbf{I}_{n}\end{matrix}\right]\in
L(2K\times\mathcal{N})\;.$

1 expresses a relationship between the eigenvalues of $\mathbf{B}$ and
$\mathbf{B}^{\prime}$. There is an associated relationship between their
eigenvectors.

######  Lemma 3.4.

Let $\mathbf{u}\in V(\vec{\mathcal{E}})$ be an eigenvector of $\mathbf{B}$
with eigenvalue $\beta$. Let
$\mathbf{x}=(\mathbf{x}_{1},\mathbf{x}_{2})^{T}\in V(2K\times\mathcal{N})$,
where $\mathbf{x}_{1},\mathbf{x}_{2}\in V(K\times\mathcal{N})$ are doubly-
indexed vectors defined entrywise by

(2) $\displaystyle
x_{1;k,i}\coloneqq\sum_{\begin{subarray}{c}iQ\in\vec{\mathcal{E}}_{k}(i)\end{subarray}}u_{iQ}\quad\text{and}\quad
x_{2;k,i}\coloneqq\sum_{\begin{subarray}{c}\begin{subarray}{c}iQ\in\vec{\mathcal{E}}_{k}(i)\end{subarray}\end{subarray}}\sum_{\begin{subarray}{c}j\in
Q\setminus i\end{subarray}}u_{jQ}\;.$

Then, either $\mathbf{x}$ is an eigenvector of $\mathbf{B}^{\prime}$ with
eigenvalue $\beta$, or $\mathbf{x}=\mathbf{0}$.

Heuristically, Lemma 3.4 states that one can aggregate cluster information in
$\mathbf{u}$ to the level of nodes by summing over all edges of size $k$ with
specified points. To obtain a first spectral clustering algorithm, we further
sum $\mathbf{x}_{1}$ over edge sizes, obtaining vector $\bar{\mathbf{x}}$ with
entries

(3) $\bar{x}_{i}\coloneqq\sum_{k\in K}x_{1;k,i}\;.$

This particular choice will be further supported when considering the
relationship of $\mathbf{B}$ to the belief-propagation algorithm on hypergraph
stochastic blockmodels, which we do in Section 5. We find experimentally that
the sign of $\bar{x}_{i}$ is most informative, and it is therefore helpful to
use the vector $\tilde{\mathbf{x}}\coloneqq\mathrm{sgn}(\bar{\mathbf{x}})$,
with the sign computed entrywise.

In the case of uniform hypergraphs, there exist recent probabilistic
guarantees ensuring that $\mathbf{x}$ is correlated with planted communities
under a certain choice of data generating process [61]. Extending these
guarantees to the nonuniform cases is an avenue of future work.

We now state our first spectral clustering algorithm, Nonbacktracking
Hypergraph Spectral Clustering (NBHSC) in Algorithm 1. We compute a desired
number of eigenvectors $\left\\{\mathbf{u}^{(\ell)}\right\\}$, form
$\tilde{\mathbf{x}}^{(\ell)}$ for each, and then cluster in the Euclidean
embedding described by the set
$\left\\{\tilde{\mathbf{x}}^{(\ell)}\right\\}_{\ell}$. 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 $\mathbf{z}$. Throughout
this paper we use the standard $k$-means algorithm [41, 58], but other choices
could in principle lead to superior performance.

Algorithm 1 Nonbacktracking Hypergraph Spectral Clustering (NBHSC)

1: $\left\\{\mathbf{u}^{(\ell)}\right\\}_{\ell=1}^{h}\leftarrow$ leading $h$
eigenvectors of $\mathbf{B}^{\prime}$

2: Initialize $\tilde{\mathbf{X}}$

3: for $\ell=1,\ldots,h$ do

4: Compute $\bar{\mathbf{x}}^{(\ell)}$ via Eqs. 2 and 3

5:
$\tilde{\mathbf{X}}_{\cdot\ell}\leftarrow\mathrm{sgn}(\bar{\mathbf{x}}^{\ell})$

6: end for

7: $\mathbf{z}=$ Cluster$(\tilde{\mathbf{X}})$

8: return $\mathbf{z}$

We will argue in Section 6.1 that NBHSC is limited in cases in which edges of
different sizes carry different cluster information in $\mathcal{H}$. This
limitation motivates our development of spectral methods based on belief-
propagation in a hypergraph stochastic blockmodel, a development to which we
turn now.

##  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. [34]. Many other related
formulations exist in the literature [29, 4, 18, 23, 61, 19]. We prove in-
expectation results for eigenpairs of the matrices $\mathbf{B}_{k}$, and pose
conjectures generalizing recent proofs by Stephan and Zhu [61] of eigenpair
concentration results in the uniform case. These conjectures will also inform
our development of a spectral method based on belief-propagation in Section 5.

Our blockmodel is a probability distribution over hypergraphs, which we denote
$\eta$. To sample from this model, we first independently assign each node
$i\in\mathcal{N}$ a label $z_{i}$ from a finite label alphabet $\mathcal{Z}$
of size $\ell$. We let $q^{(s)}$ give the probability that $z_{i}=s$. We
collect the labels in a vector $\mathbf{z}$. Let $\mathcal{R}$ give the set of
possible edges, which we usually take to be sets of nodes with some specified
possible sizes. Let $\mathcal{R}_{k}$ denote the set of node subsets of size
$k$. Define $\mathcal{R}(i)$ to be the set of subsets containing node $i$, and
let $\mathcal{R}_{k}(i)=\mathcal{R}_{k}\cap\mathcal{R}(i)$. To realize edges,
we consider each set of nodes $R\in\mathcal{R}$ and add this set to the edge
set $\mathcal{E}$ with probability $\eta(R\in\mathcal{E}|\mathbf{z}_{R})$
depending on the labels in $\mathbf{z}_{R}$. Here, we will focus on the sparse
setting in which

(4)
$\eta(R\in\mathcal{E}|\mathbf{z}_{R})=\frac{\omega(\mathbf{z}_{R})}{n^{\lvert
R\rvert-1}}$

for some function $\omega$ 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 $\mathbf{z}$ and edge set
$\mathcal{E}$ is

(5)
$\eta(\mathcal{E},\mathbf{z})=\left(\prod_{i\in\mathcal{N}}q^{(z_{i})}\right)\left(\prod_{R\in\mathcal{R}}\eta(R\in\mathcal{E}|\mathbf{z}_{R})\right)\;.$

In our development below, it will be useful to reason asymptotically about
many quantities. Define $x\doteq 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 $\eta$ as $n$ grows large and the largest edge-size
$k$ remains fixed.

Let

(6)
$c_{k}^{(s)}\coloneqq\frac{1}{(k-1)!}\sum_{\mathbf{z}\in\mathcal{Z}^{k-1}}\omega(\mathbf{z},s)\prod_{t\in\mathbf{z}}q^{(t)}\;.$

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:

$\sum_{\begin{subarray}{c}R\in\mathcal{R}_{k}(i)\end{subarray}}\eta(R\in\mathcal{E}|\mathbf{z}_{R})\prod_{j\in
R\setminus
i}q^{(z_{j})}\doteq\binom{n-1}{k-1}\sum_{\mathbf{z}\in\mathcal{Z}^{k-1}}n^{1-k}\omega(\mathbf{z},s)\prod_{t\in\mathbf{z}}q^{(t)}=(1+O(n^{-1}))c_{k}^{(s)}\doteq
c_{k}^{(s)}\;.$

We also define

(7)
$c_{k}^{(s,t)}\coloneqq\frac{1}{(k-2)!}\sum_{\mathbf{z}\in\mathcal{Z}^{k-2}}\omega(\mathbf{z},s,t)\prod_{t\in\mathbf{z}}q^{(t)}\;,$

which gives the expected number of cluster $t$ neighbors of a node in cluster
$s$ through edges of size $k$. We have the identity

(8) $c_{k}^{(s)}=\frac{1}{k-1}\sum_{t\in\mathcal{Z}}q^{(t)}c_{k}^{(s,t)}\;,$

which computes $c_{k}^{(s)}$ by conditioning on the label of a second node in
each possible edge.

###  4.1 Community-Correlated Eigenvectors

We now aim to describe eigenvectors of $\mathbf{B}_{k}$ under the blockmodel
Eq. 5 which are correlated with community structure. We implicitly condition
on a realization of the label vector $\mathbf{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)}=\frac{1}{2}$. We further assume that the empirical
distribution of group labels in the realized label vector $\mathbf{z}$ is
close to $q$; this holds with high probability as $n$ grows large.

  * •

The group-specific expected degrees are equal:
$c_{k}^{(1)}=c_{k}^{(2)}\coloneqq c_{k}$ for each $k$.

  * •

We have $c_{k}^{(1,1)}=c_{k}^{(2,2)}\coloneqq c_{k}^{\mathrm{in}}$ and
$c_{k}^{(1,2)}=c_{k}^{(2,1)}\coloneqq c_{k}^{\mathrm{out}}$ for all $k$.

Formal concentration results under these hypotheses are available for graphs
[14] and $k$-uniform hypergraphs [61]. 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 therefore provide
informal results by reasoning in expectation.

Let $\mathbf{u}\in V(\vec{\mathcal{E}})$ have entries $u_{iQ}=\lvert
Q\rvert-1$. Let $\boldsymbol{\sigma}\in V(\mathcal{N})$ be the vector with
entries

$\sigma_{i}\coloneqq\begin{cases}+1&\quad z_{i}=1\\\ -1&\quad
z_{i}=2\;.\end{cases}$

Let $\mathbf{v}\in V(\vec{\mathcal{E}})$ have entries
$v_{iQ}\coloneqq\sum_{j\in Q\setminus i}\sigma_{j}$.

######  Theorem 4.1.

Consider a hypergraph sampled from $\eta$. Let $\alpha_{k}=c_{k}(k-1)$. Let
$\beta_{k}=\frac{c_{k}^{\mathrm{in}}-c_{k}^{\mathrm{out}}}{2}$. Then, for all
tuples $Q$ and nodes $i\in Q$, we have

(9)
$\displaystyle\mathbb{E}\left[(\mathbf{B}_{k}\mathbf{u})_{iQ}-\alpha_{k}u_{iQ}|Q\in\mathcal{E}\right]\doteq
0\;,$ (10)
$\displaystyle\mathbb{E}\left[(\mathbf{B}_{k}\mathbf{v})_{iQ}-\beta_{k}v_{iQ}|Q\in\mathcal{E}\right]\doteq
0\;.$

The full proof of Eq. 10 is provided in Appendix C; the proof of Eq. 9 is
similar.

In expectation, Theorem 4.1 says that, $\mathbf{B}_{k}$ has a Perron
eigenvector with eigenvalue $\alpha_{k}$, and a community-correlated
eigenvector with eigenvalue $\beta_{k}$. To substantiate the claim that
$\mathbf{v}$ is indeed community-correlated, we can sum over $k-$edges
incident to $i$, obtaining

(11)
$\mathbb{E}\left[\sum_{Q\in\mathcal{E}_{k}(i)}v_{iQ}\right]=\sum_{Q\in\mathcal{R}_{k}(i)}\eta(Q\in\mathcal{E}_{k})\sum_{j\in
Q\setminus i}\sigma(j)\doteq\frac{1}{n}\sum_{j\neq
i}\sigma_{j}c_{k}^{(z_{i},z_{j})}\doteq\beta_{k}\sigma_{i}\;.$

This calculation is similar to the proof of Lemma 3.4. The expectation of
$\mathbf{v}$, therefore, can be used to recover the vectors
$\boldsymbol{\sigma}$ and $\mathbf{z}$. In practice, we do not have access to
the expectation and our estimates of $\boldsymbol{\sigma}$ are noisy. Formal
concentration results along the lines of Bordenave et al. [14] and Stephan and
Zhu [61] are necessary to provide more guarantees in this case.

##  5 Nonbacktracking Operators and Belief Propagation

We now derive a connection between nonbacktracking operators and the belief-
propagation algorithm (BP) [11] for community detection in sparse hypergraphs.
Our derivation extends arguments by Krzakala et al. [37] 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, of
which sparse stochastic blockmodels are an example. Following standard
arguments [37], 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 $\mathbf{B}_{k}$—whose eigenvectors may therefore contain cluster
information.

###  5.1 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 $\mathcal{E}$ but not the label vector $\mathbf{z}$. We would like to
obtain information about the conditional distribution
$\eta(\mathbf{z}|\mathcal{E})=\frac{\eta(\mathcal{E},\mathbf{z})}{\eta(\mathcal{E})}$
of labels given the edge data $\mathcal{E}$. An especially relevant summary is
the marginal distribution of labels for each node, $\eta(z_{i}|\mathcal{E})$.

We will estimate the marginals using belief-propagation [11]. We treat both
subsets $R\in\mathcal{R}$ and labels $z_{i}\in\mathcal{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
$\zeta_{i}^{(s)}=q^{(s)}$. The tuple factors are somewhat more complex. Let
$\mu_{iR}^{(s)}$ denote the message that node $i$ passes to the factor $R$
expressing its belief that $z_{i}=s$. Let $\nu_{Ri}^{(s)}$ 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

(12) $\displaystyle\mu_{iR}^{(s)}$
$\displaystyle\leftarrow\frac{1}{Z_{iR}}q^{(s)}\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus
R\end{subarray}}\nu_{Qi}^{(s)}$ (13) $\displaystyle\nu_{Ri}^{(s)}$
$\displaystyle\leftarrow\frac{1}{Z_{Ri}}\sum_{\mathbf{z}_{R}:z_{i}=s}\eta(R\in\mathcal{E}|\mathbf{z}_{R})\prod_{j\in
R\setminus i}\mu_{jR}^{(z_{j})}\;.$

Here, $Z_{iR}$ and $Z_{Ri}$ are normalizing constants ensuring that
$\sum_{s}\mu_{iR}^{(s)}=\sum_{s}\nu_{Ri}^{(s)}=1$.

On factor graphs which are trees, the updates Eqs. 12 and 13 converge, and the
desired marginals can be obtained by computing the marginal message for node
$i$:

(14)
$\mu_{i}^{(s)}\coloneqq\frac{1}{Z_{i}}q^{(s)}\prod_{Q\in\mathcal{R}(i)}\nu_{Qi}^{(s)}\;.$

Our factor graph is admittedly not a tree. One possibility is to modify the
belief propagation algorithm to account for loops [36]. We instead follow the
standard argument that:

  * •

The sparsity assumption on $\eta$ implies that the realized hypergraph is
locally tree-like as $n$ grows large.

  * •

The updates Eqs. 12 and 13 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-treelike factor graph. The first point is
discussed in the graph setting by [20]. The second has been argued
heuristically in several papers [20, 37, 4]. We provide below a heuristic
statement (2) that expresses the thrust of this argument. In Theorem B.1
(Appendix B) we offer a rigorous statement with proof.

###  5.2 The BP Jacobian

For each $k$, define the matrix $\mathbf{G}_{k}\in L(\mathcal{Z},\mathcal{Z})$
with entries

(15) $g_{k;s,t}\coloneqq
q^{(s)}\left(\frac{c_{k}^{(s,t)}}{(k-1)c_{k}^{(s)}}-1\right)\;.$

Finally, let $\mathbf{F}$ be the function with components

$\mathbf{F}(\boldsymbol{\mu})_{iR}^{(s)}\coloneqq\frac{1}{Z_{iR}}q^{(s)}\prod_{Q\in\mathcal{R}(i)\setminus
R}\sum_{\mathbf{z}_{Q}:z_{i}=s}\eta(a_{Q}|\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{jQ}^{(z_{j})}\;.$

This expression is obtained by substituting Eq. 13 into Eq. 12. We can then
write the dynamics, restricted to the variable $\boldsymbol{\mu}$, as
$\boldsymbol{\mu}\leftarrow\mathbf{F}(\boldsymbol{\mu})$.

We are now prepared to relate the belief-propagation algorithm to the
nonbacktracking matrices $\mathbf{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.

######  Claim 2.

Let $\mathcal{H}$ be sampled from a sparse Bernoulli stochastic blockmodel
$\eta$ in which $c_{k}^{(s)}=c_{k}^{(t)}$ for all $s,t\in\mathcal{Z}$ and
$k\in K$. Then, with high probability as $n\rightarrow\infty$:

  * •

The point $\bar{\boldsymbol{\mu}}$ with coordinates
$\bar{\mu}_{iR}^{(s)}=q^{(s)}$ for all $i\in\mathcal{N}$,
$R\in\mathcal{R}(i)$, and $s\in\mathcal{Z}$ is approximately a fixed point of
$\mathbf{F}$, in the sense that
$\mathbf{F}(\bar{\boldsymbol{\mu}})\doteq\bar{\boldsymbol{\mu}}$.

  * •

The Jacobian of $\mathbf{F}$ at $\bar{\boldsymbol{\mu}}$, 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 $\mathbf{J}+O(n^{-1})$, where

(16) $\mathbf{J}\coloneqq\sum_{k\in K}\mathbf{G}_{k}\otimes\mathbf{B}_{k}\in
L(\mathcal{Z}\times\vec{\mathcal{E}})\;.$

2 tells us that the uninformative point $\bar{\boldsymbol{\mu}}$ is
approximately a fixed point of the belief-propagation update $\mathbf{F}$. In
turn, the stability of this fixed point is governed by the Jacobian matrix,
which is dominated by a block which can be approximated by $\mathbf{J}$.
Because $\mathbf{J}$ approximates the BP dynamics in a neighborhood of
$\bar{\boldsymbol{\mu}}$, and because the BP dynamics by design attempt to
find community structure, we expect that the leading eigenvectors of
$\mathbf{J}$ carry information about community structure in the hypergraph
$\mathcal{H}$.

In general, $\mathbf{J}\in
L(\mathcal{Z}\times\vec{\mathcal{E}})\simeq\mathbb{R}^{\ell\vec{m}\times\ell\vec{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 leading eigenvalues of
$\mathbf{J}=\mathbf{G}_{k}\otimes\mathbf{B}_{k}$ are of the form
$\gamma\beta$, where $\gamma$ is an eigenvalue of $\mathbf{G}_{k}$ and $\beta$
an eigenvalue of $\mathbf{B}_{k}$. This allows Angelini et al. [4] to compute
only on $\mathbf{B}_{k}$, or indeed on the smaller $\mathbf{B}_{k}^{\prime}$.
In nonuniform hypergraphs, however, such a simple reduction is not available.
We therefore offer a partial generalization in Theorem 5.1 that enables us to
compute on a smaller matrix on nonuniform hypergraphs.

Let $\mathbf{u}\in V(\mathcal{Z}\times\vec{\mathcal{E}})$. Index the entries
of $\mathbf{u}$ as $u_{iQ}^{(s)}$, where $iQ\in\vec{\mathcal{E}}$ is a pointed
edge and $s\in\mathcal{Z}$ is a group label. Let $\mathbf{L}\in
L(\mathcal{Z}\times\vec{\mathcal{E}},2\mathcal{Z}\times K\times\mathcal{N})$
be $\mathbf{L}:\mathbf{u}\mapsto(\mathbf{x}_{1},\mathbf{x}_{2})^{T}\in
V(2\mathcal{Z}\times K\times\mathcal{N})$, where

(17)
$x_{1;ik}^{(s)}\coloneqq\sum_{Q\in\mathcal{E}_{k}(i)}u_{iQ}^{(s)}\quad\text{and}\quad
x_{2;ik}^{(s)}\coloneqq\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in Q\setminus
i}u_{jQ}^{(s)}\;.$

Let $\mathbf{M}\in
L(\mathcal{Z}\times\vec{\mathcal{E}},\mathcal{Z}\times\mathcal{E})$ send
$\mathbf{x}\mapsto\mathbf{t}$, where $t_{Q}^{(s)}=\sum_{i\in Q}u_{iQ}^{(s)}$.
Finally, let $\mathbf{N}\in L(\mathcal{Z}\times\mathcal{E})$ send
$\mathbf{t}\mapsto\bar{\mathbf{t}}$, where

$\bar{t}_{Q}^{(s)}=(1-\lvert Q\rvert)\sum_{s^{\prime}\in\mathcal{Z}}g_{\lvert
Q\rvert}^{(s,s^{\prime})}t_{Q}^{(s^{\prime})}\;.$

We can now partition the eigenpairs of $\mathbf{J}$ into two subsets in
fashion similar to that of Theorem 3.3, although the proof is somewhat
different.

######  Theorem 5.1.

There exists a matrix $\mathbf{J}^{\prime}\in L(2\mathcal{Z}\times
K\times\mathcal{N})$ which depends only on $\mathcal{H}$ and the parameter
matrices $\\{\mathbf{G}_{k}\\}$, such that any eigenpair $(\xi,\mathbf{u})$ of
$\mathbf{J}$ satisfies exactly one of the following two alternatives:

  1. 1.

$(\xi,\mathbf{L}\mathbf{u})$ is an eigenpair of $\mathbf{J}^{\prime}$.

  2. 2.

$(\xi,\mathbf{M}\mathbf{u})$ is an eigenpair of $\mathbf{N}$.

The proof is similar to that of Lemma 3.4. We supply it, as well as an
explicit formula for $\mathbf{J}^{\prime}$, in Appendix D. In Algorithm 2
below, we focus on the eigenpairs in Case 1 of Theorem 5.1. These eigenpairs
are directly interpretable, since $u_{ik}^{(s)}$ sums beliefs that node $i$
belongs to cluster $s$ across edges of size $k$. One reason for this choice is
that, $\bar{\mathbf{t}}$ carries no information at the level of nodes, and
therefore cannot be directly used for node clustering. Since we have
$\bar{\mathbf{t}}=\mathbf{N}\mathbf{t}=\xi\mathbf{t}$ in Case 2, it follows
that $\mathbf{t}$ also cannot be directly applied for node-level clustering
tasks. The use of $\mathbf{t}$ for edge-level clustering tasks might be an
interesting direction of future work. Generalized results along the lines of
Theorem 2 of Stephan and Zhu [61] could provide probabilistic guarantees on
the eigenvectors satisfying Case 1.

###  5.3 Alternating Belief-Propagation Spectral Clustering

2 and 5.1 jointly suggest a modified algorithm based on the eigenvectors of
$\mathbf{J}$ or $\mathbf{J}^{\prime}$ rather than the eigenvectors of
$\mathbf{B}$ or $\mathbf{B}^{\prime}$. Since $\mathbf{J}^{\prime}$ depends on
the blockmodel parameters through the matrices $\\{\mathbf{G}_{k}\\}$, we
alternate between spectral clustering steps and updates to these parameters.
This alternating structure is reminiscent of expectation-maximization [21] and
other coordinate-ascent algorithms. However, our alternating algorithm is not
literally a form of coordinate-ascent 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 $\eta(R\in\mathcal{E}|\mathbf{z}_{R})$ for every possible
combination of labels $\mathbf{z}_{R}$. In a hypergraph with edges up to size
$k$ and $\ell$ group labels, there are $\frac{\ell^{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 $\eta(R\in\mathcal{E}|\mathbf{z}_{R})$, but only the
entries of $\mathbf{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 $\ell+\frac{1}{2}k\ell^{2}$ parameters,
a number which scales much more favorably than $\frac{\ell^{k}}{k!}$. We show
how to estimate these parameters given an estimated label vector in Appendix
F. As a result, BPHSC is both less sensitive to the fine details of the HSBM
parameters and more robust against overfitting concerns.

The eigenvectors of $\mathbf{J}^{\prime}$ are elements of
$V(2\mathcal{Z}\times K\times\mathcal{N})$. In order to carry out the
clustering step, we need to obtain from these eigenvectors a set of feature
vectors in $V(\mathcal{N})$ that carry information on the level of nodes. To
do so, we use Eq. 14, 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 $\mathbf{u}_{j}$ we form the matrix
$\tilde{\mathbf{X}}_{j}\in L(\mathcal{N},\mathcal{Z})$ with entries
$\tilde{x}_{j;i,s}=\mathrm{sign}(\sum_{k}x_{1;ik}^{(s)})$, with
$x_{1;ik}^{(s)}$ as in Eq. 17. The columns of the matrices
$\\{\tilde{\mathbf{X}}_{j}\\}$ form the embedding coordinates on which we
cluster. This clustering step gives a new label vector $\mathbf{z}$. We use
$\mathbf{z}$ to compute updated estimates of the model parameters, and repeat.
A single stage of the clustering step is illustrated in Figure 1.

Algorithm 2 Step of Alternating BP Hypergraph Spectral Clustering (BPHSC)

0: Hypergraph $\mathcal{H}$, current clustering $\mathbf{z}_{0}$

1:
$\\{\mathbf{G}_{k}\\}\leftarrow\mathrm{estimateParameters}(\mathcal{H},\mathbf{z}_{0})$

2: Form $\mathbf{J}^{\prime}$ as in Theorem 5.1.

3: $\\{\mathbf{u}^{(\ell)}\\}_{\ell=1}^{h}\leftarrow$ leading  $h$
eigenvectors of $\mathbf{J}^{\prime}$ with real eigenvalues.

4: Initialize $\tilde{\mathbf{X}}$.

5: for $j=1,\ldots,h\;,s=1,\ldots,\ell\;,i=1,\ldots,n$ do

6: $x_{1;ik}^{(s)}\leftarrow\sum_{Q\in\mathcal{E}_{k}(i)}u_{iQ}^{(s)}$

7: $\tilde{x}_{j;1;i,s}\leftarrow\mathrm{sign}(\sum_{k\in K}x_{1;ik}^{(s)})$

8: end for

9: $\mathbf{z}=$ Cluster$(\tilde{\mathbf{X}})$

9: $\mathbf{z}$.

In practice, one must choose both how many clusters to form and how many
iterations to perform before accepting the given clustering. This decision is
complicated by the fact that this clustering algorithm does not optimize a
likelihood or other objective function. Indeed, BPHSC does not even require us
to fully define a likelihood, even though we derived it from a hypergraph
stochastic blockmodel. This is advantageous, as computing a likelihood
requires the specification of a generative model with explicit affinity
function $\omega$ describing the rates of edges on node subsets containing
different cluster labels. Such a specification poses problems, as the number
of parameters required to specify $\omega$ grows rapidly with the number of
clusters and sizes of possible hyperedges. Importantly, BPHSC does not require
$\omega$ itself, but only the pairwise relationships encoded by
$\\{\mathbf{G}_{k}\\}$.

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 $\mathrm{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 screen plots can assist choices
about the correct number of cluster labels to use. The
$\mathrm{estimateParameters}(\mathcal{H},\mathbf{z}_{0})$ step can be
performed by maximum-likelihood estimation. We describe this step in Appendix
F.

Figure 1:  Illustration of the stages of BPHSC (Algorithm 2). In this example,
we generated a synthetic hypergraph with three clusters of 50 nodes each, with
$c_{2}=c_{3}=5$. 90% of 2-edges and 10% of 3-edges are within-cluster. (a):
First, we compute the spectrum of the matrix $\mathbf{J}^{\prime}$ described
in Theorem 5.1. There are eight real eigenvalues larger than unity in
magnitude. (b): Next, we form aggregated eigenvectors of $\mathbf{J}$, which
we visualize here as matrices. Shown are the signs of the entries of
aggregated eigenvectors for highlighted real eigenvalues in (A). The entries
are $\tilde{x}_{j;i,s}=\mathrm{sign}\left(\sum_{k}x_{ik}^{(s)}\right)$, with
$\mathbf{x}$ calculated from $\mathbf{u}_{j}$ as in Eq. 17. Thick horizontal
lines separate nodes in different ground-truth clusters. (c): Finally, we
study the nodes in the space spanned by the columns $\mathbf{x}_{j;\cdot,s}$
for each aggregate eigenvector $\mathbf{x}_{j}$. In our example, this space
has $8\times 3=24$ dimensions. In a full run of Algorithm 2, 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.

##  6 Thresholds for Spectral Clustering

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 [51, 37]. We first restate a standard
definition of _detection_ in clustering problems. For each $n$, let
$\eta_{n}(\theta,\mathbf{z}_{n})$ be a probability distribution over graphs on
$n$ nodes. Each such distribution possesses the same shared set of parameters
$\theta$, as well as a planted partition $\mathbf{z}_{n}$ of nodes. Let
$\mathcal{G}_{n}\sim\eta_{n}(\theta,\mathbf{z}_{n})$. Let $\mathtt{A}$ be a
clustering algorithm, which we view as a map
$\mathcal{G}_{n}\mapsto\hat{\mathbf{z}}_{n}$ from the data to an estimate of
the true labels. Let $\rho$ be a correlation function measuring label
agreement with the property that if $\mathbf{z}_{n}$ and
$\hat{\mathbf{z}}_{n}$ are length-$n$ labels sampled independently, then
$\rho(\mathbf{z},\hat{\mathbf{z}})\rightarrow 0$ as $n\rightarrow\infty$.
Examples of measures include mutual information, Adjusted Rand, and suitably
adjusted versions of the overlap [20, 1].

######  Definition 6.1.

Algorithm $\mathtt{A}$ _detects communities_ in the sequence
$\\{\eta_{n},\theta,\mathbf{z}_{n}\\}$ with respect to the correlation
function $\rho$ if there exists $\epsilon>0$ such that, with high probability
as $n$ grows large,
$\rho(\mathbf{z}_{n},\mathtt{A}(\mathcal{G}_{n}))>\epsilon$.

A particularly notable case is the graph stochastic blockmodel with two
equally-sized communities, where $\theta=(a,b)$ contains $a$, the expected
number of edges joining a node to others in its own community, and $b$, the
expected number of edges joining a node to nodes in the opposite community.
For some algorithms $\mathtt{A}$, $\mathtt{A}$ is able to detect clusters in
$\\{\eta_{n},\theta,\mathbf{z}_{n}\\}$ if

(18) $\phi\coloneqq\frac{(a-b)^{2}}{2(a+b)}>1\;.$

and unable if $\phi<1$ [20, 51, 37]. Mossel et al. [48] and Massoulié [44]
later proved that _no_ choice of algorithm $\mathtt{A}$ will suffice for
$\phi<1$, but that some choices will indeed succeed for $\phi>1$.

Definition 6.1 generalizes to the setting of hypergraphs with planted cluster
structure. It is necessary only to allow $\eta_{n}(\theta,\mathbf{z}_{n})$ to
be a probability distribution over hypergraphs. Indeed, Angelini et al. [4]
have offered conjectures concerning the ability of non-backtracking spectral
methods to detect planted clusters in the setting of uniform hypergraphs.
Several of these conjectures were recently proven by Stephan and Zhu [61]. 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.

###  6.1 Thresholds for NBHSC

In Algorithm 1, we form the nonbacktracking operator $\mathbf{B}$, obtain a
spectral embedding of the nodes, and then cluster in the embedding space. The
spectral embedding is obtained by considering the $h$ leading eigenvectors of
$\mathbf{B}$. A necessary condition for detection is that at least one of
these eigenvectors is correlated with ground-truth community labels. Theorem
4.1 implies in-expectation eigenrelations for this case.

######  Corollary 3.

Under the hypotheses as Theorem 4.1, and with $\mathbf{u}$ and $\mathbf{v}$ as
defined there,

(19) $\displaystyle\mathbb{E}\left[(\mathbf{B}\mathbf{u})_{iQ}-\alpha
u_{iQ}|Q\in\mathcal{E}\right]$ $\displaystyle\doteq 0\;$ (20)
$\displaystyle\mathbb{E}\left[(\mathbf{B}\mathbf{v})_{iQ}-\beta
v_{iQ}|Q\in\mathcal{E}\right]$ $\displaystyle\doteq 0\;,$

where $\alpha\coloneqq\sum_{k\in K}\alpha_{k}$ and $\beta\coloneqq\sum_{k\in
K}\beta_{k}$.

It was argued informally by Angelini et al. [4] and recently proven by Stephan
and Zhu [61] that, in the uniform case, eigenvalues other than $\alpha$ and
$\beta$ concentrate with high probability in the disc of radius
$\sqrt{\alpha}$. The informal argument generalizes smoothly. A formal
generalization alongside concentration results would be sufficient to prove
the following conjecture:

######  Conjecture 4.

In the blockmodel $\eta$, with high probability as $n\rightarrow\infty$,

  1. 1.

The leading eigenvalue of $\mathbf{B}$ is $\alpha\doteq\sum_{k\in K}c_{k}$.

  2. 2.

There exists a community-correlated eigenvector with an associated eigenvalue
$\beta\doteq\frac{1}{2}\sum_{k\in
K}\left(c_{k}^{\mathrm{in}}-c_{k}^{\mathrm{out}}\right)$.

  3. 3.

The bulk of the spectrum of $\mathbf{B}$ is confined to the disk of radius
$\sqrt{\beta}$.

  4. 4.

Spectral clustering using $\mathbf{B}$ is able to detect ground-truth clusters
iff

(21) $\beta^{2}>\alpha\;.$

These conjectures parallel known results for graphs [14] and uniform
hypergraphs [61].

Conjecture 4 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^{\prime}$ in which edges of
size $k$ tend to be within-cluster ($\beta_{k}>0$), while edges of size
$k^{\prime}$ tend to be between-cluster $(\beta_{k^{\prime}}<0$). In such a
case, the approximate eigenvalue $\beta=\beta_{k}+\beta_{k^{\prime}}$ of
$\mathbf{B}$ can be smaller in magnitude than either of $\beta_{k}$ and
$\beta_{k^{\prime}}$, and could potentially be fully lost in the eigenvalue
bulk, causing spectral detection to fail. This can occur even when $\beta_{k}$
and $\beta_{k^{\prime}}$ are sufficiently large in magnitude that detection
based on edges of size $k$ or $k^{\prime}$ alone would succeed. The need to
adaptively synthesize signals from hyperedges of differing sizes is one
motivation of spectral methods based on the belief-propagation Jacobian.

###  6.2 Thresholds for BPHSC

We can also use Theorem 4.1 to motivate a conjecture about the spectral
structure of the Jacobian matrix $\mathbf{J}$. For each $k$, the vector
$\mathbf{1}=(1,1)$ is an eigenvector of $\mathbf{G}_{k}$ with eigenvalue
$\gamma_{k}=\frac{\beta_{k}}{(k-1)c_{k}}$. Let
$\hat{\mathbf{u}}=\mathbf{1}\otimes\mathbf{u}$ and
$\hat{\mathbf{v}}=\mathbf{1}\otimes\mathbf{v}$. Let $\nu=\sum_{k\in
K}\gamma_{k}\alpha_{k}$ and $\lambda=\sum_{k\in K}\gamma_{k}\beta_{k}$.
Equations 16 and 4.1 now imply the following result:

######  Corollary 5.

Under the two-group blockmodel $\eta$, for $\ell=1,2$, we have

(22)
$\displaystyle\mathbb{E}\left[(\mathbf{J}\hat{\mathbf{u}})_{iQ}^{(\ell)}-\nu\hat{u}_{iQ}^{(\ell)}|Q\in\mathcal{E}\right]$
$\displaystyle\doteq 0$ (23)
$\displaystyle\mathbb{E}\left[(\mathbf{J}\hat{\mathbf{v}})_{iQ}^{(\ell)}-\lambda\hat{v}_{iQ}^{(\ell)}|Q\in\mathcal{E}\right]$
$\displaystyle\doteq 0\;.$

Motivated by this result, we pose the following conjecture:

######  Conjecture 6.

In the same setting as Theorem 4.1, with high probability as $n$ grows large,

  * •

$\mathbf{J}$ possesses a real, community-correlated eigenvector with
eigenvalue $\lambda+o(1)$.

  * •

BPHSC, when initialized with knowledge of the true parameters
$c_{k}^{\mathrm{in}}$ and $c_{k}^{\mathrm{out}}$, is able to detect ground-
truth clusters if $\lvert\lambda\rvert>1$.

###  6.3 Parameterized Thresholds

We now illustrate Conjectures 4 and 6 through some computational experiments.
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 [56]. 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 will now show that, in this model, Eq. 21 in Conjecture 4 defines a pair of
hyperplanes in the coordinates $\\{p_{k}\\}_{k\in K}$. This will follow from
two facts. First, by construction, $\alpha$ does not depend on $\\{p_{k}\\}$.
Second, $\beta$ is an affine function of $\\{p_{k}\\}$:

######  Lemma 6.2.

In this model, we have the following with
$r_{k}\coloneqq\frac{1-2^{2-k}}{2-2^{2-k}}$

(24) $\beta=\sum_{k\in K}(k-1)c_{k}\left[2(1-r_{k})p_{k}+2r_{k}-1\right]\;.$

The proof is a direct calculation and provided in Appendix E. Lemma 6.2 in
conjunction with Conjecture 4 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 2(a). shows an experiment on a blockmodel on $200$ 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 Eqs. 21 and 6.2. Under Conjecture 4, as
$n\rightarrow\infty$, spectral clustering returns a clustering 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.

Figure 2: Experiments and analytical boundaries for nonbacktracking spectral
clustering in synthetic hypergraphs of 200 nodes. Each pixel is the mean
Adjusted Rand Index (ARI) of the retrieved clustering against ground truth
across 20 trials. (a). Spectral clustering via the nonbacktracking operator
$\mathbf{B}$ according to NBHSC (Algorithm 1), $c_{2}=c_{3}=5$. The white
lines are given by Eqs. 21 and 24. (b)-(d). Spectral clustering in one round
of BPHSC (Algorithm 2) using the true values of the parameter matrices
$\left\\{\mathbf{G}_{k}\right\\}_{k\in K}$. In each panel, the white ellipse
has parameters given by Eq. 26. The dashed black curves trace the ellipse
described by the collision of the community-correlated eigenvalue $\lambda$
with the eigenvalue $\nu$ described in 5. We show hypergraphs with edges of
size $2$, $3$, and $4$. The mean $k$-degrees for each edge-size $k$ vary in
each panel. (b). $c_{2}=c_{3}=5$, $c_{4}=0$. (c). $c_{3}=c_{4}=5$, $c_{2}=0$.
(d). $c_{2}=5$, $c_{4}=50$, $c_{3}=0$.

###  6.4 Thresholds for BPHSC

We now consider belief-propagation hypergraph spectral clustering. In this
algorithm, it is necessary to form an estimate of the parameter matrices
$\left\\{\mathbf{G}_{k}\right\\}_{k\in K}$. In Algorithm 2, 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 $\left\\{\mathbf{G}_{k}\right\\}_{k\in K}$ are known _a priori_. The
assumption that the parameters are known exactly is sometimes called the
Nishimori condition in statistical physics [20]. In graph SBMs, the relaxation
of the Nishimori condition leads to a much more complicated analysis with
qualitatively different conclusions [33]. We now proceed with analysis under
the Nishimori condition.

We can compute the entries of the matrix $\mathbf{G}_{k}$ in terms of
$c_{k}^{\mathrm{in}}$ using Eq. 15, obtaining

$g_{k}^{\mathrm{in}}\coloneqq
g_{k}^{(s,s)}=\frac{1}{2}\left(\frac{c_{k}^{\mathrm{in}}}{(k-1)c_{k}}-1\right)\quad\text{and}\quad
g_{k}^{\mathrm{out}}\coloneqq g_{k}^{(s,t)}=-g_{k}^{\mathrm{in}}\;$

The eigenvalue of $\mathbf{G}_{k}$ that appears in 5 is then

(25)
$\gamma_{k}=g_{k}^{\mathrm{in}}-g_{k}^{\mathrm{out}}=2g_{k}^{\mathrm{in}}=\frac{c_{k}^{\mathrm{in}}}{(k-1)c_{k}}-1\;.$

Direct computation now shows that the condition $\lambda=\sum_{k\in
K}\gamma_{k}\beta_{k}=1$ defines an axis-aligned ellipsoid with coordinates
$(p_{k_{1}},p_{k_{2}},\ldots)$, centroid $(x_{k_{1}},x_{k_{2}},\ldots)$ and
radii $(a_{k_{1}},a_{k_{2}},\ldots)$, where

(26) $x_{k}=\frac{1-2r_{k}}{2-2r_{k}}\quad\text{and}\quad
a_{k}=\frac{\sqrt{(k-1)c_{k}}}{2-2r_{k}}\quad\text{with}\quad
r_{k}\coloneqq\frac{1-2^{2-k}}{2-2^{2-k}}\;.$

Conjecture 6 claims that BPHSC succeeds outside this ellipse ($\lambda>1$) and
fails inside ($\lambda<1$).

Panels (b)-(d). of Figure 2 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 200 nodes, with varying parameters $\\{p_{k}\\}$. In
each experiment, the ellipse defined by Eq. 26 is shown in white. An
implication of our conjectures is that, as $n\rightarrow\infty$, 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 $\left\\{\mathbf{G}_{k}\right\\}_{k\in 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 5 we know that there is also an approximate eigenpair
$(\nu,\hat{\mathbf{u}})$ which is uncorrelated with planted cluster structure.
The equation $\lambda=\nu$ 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 2(b-d), dashed black curves).

##  7 Experiments on Data

We first study several data sets in which ground-truth labels are available.
The contact-primary-school [60, 8] 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 3:  Cluster recovery in the contact-primary-school data set [60, 8]. We
ran BPHSC on the data for 10 rounds, using the 30 leading eigenvectors of the
belief-propagation Jacobian and with a varying number of clusters to be
estimated. In each round, we update the estimate of the labels
$\hat{\mathbf{z}}$ 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 $\mathbf{C}_{k}$ for varying edge size $k$. A similar experiment
for the contact-high-school data set [45, 8] is given in Figure 6.

Figure 3 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 $\ell$ to be learned. We repeated this process for both the original
hypergraph data and the projected (clique-expansion) graph obtained by
replacing each $k$-hyperedge 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 $\ell$, 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 $\mathbf{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, neither algorithm
achieves the perfect cluster recovery obtainable via other methods [18].

Figure 4:  As in Figure 3, using the senate-bills data set [26, 27, 8].

The senate-bills data set 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, 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. [18], 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.

Figure 5:  (Best viewed in color): Clustering the tags-math-sx data set [8]
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 2
using $\ell=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.

It is also possible to use BPHSC to cluster data in which no labels are
natively present. Figure 5 shows an experiment in which we use Algorithm 2 to
cluster tags on the forum Math StackExchange [8]. 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,\ldots,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 5.

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 $\mathbf{B}$ proposed by Storm [62]. We have considered
both a simple spectral algorithm (NBHSC) that relies only on the eigenvectors
of $\mathbf{B}$, as well as a more complex spectral algorithm (BPHSC) based on
the eigenvectors of the belief-propagation Jacobian $\mathbf{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. 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.1 yields
a matrix $\mathbf{J}^{\prime}$ of size $2n\bar{k}\ell$, which can become
computational challenging when either $k$ or $\ell$ are large.

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 2, 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 $\\{\mathbf{B}_{k}\\}$ and their associated
eigenvectors. These conjectures are inspired by known results for the graph
case [14]. 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 [61].

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 [48, 44]. We
pose a similar conjecture for the hypergraph blockmodel: within the ellipsoid
described by Eq. 26, the cluster detection problem cannot be solved by any
algorithm. We propose a proof of this conjecture as a direction of future
work.

## 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 [9]. Visualizations were constructed using ggplot2 for the R
programming language [55, 69].

## 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 [22, 15, 42, 25, 63, 24, 68, 40].

In the spirit of Zurn et al. [71], 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
subjective and subject to error. Furthermore, the heuristic that the first and
last authors correspond 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 20% had a non-male first author and 18% had a non-male last author. Of
those with at least two authors, 26% had either a non-male first author or a
non-male last author.

## References

  * Abbe [2017] E. Abbe.  Community detection and stochastic block models: Recent developments.  _The Journal of Machine Learning Research_ , 18(1):6446–6531, 2017.
  * Aktas et al. [2021] M. E. Aktas, T. Nguyen, J. Sidra, R. Rakin, and A. Esra.  Identifying critical higher-order interactions in complex networks.  _Scientific Reports_ , 11, 2021.  doi:http://dx.doi.org/10.1038/s41598-021-00017-y.
  * Alon et al. [2007] N. Alon, I. Benjamini, E. Lubetzky, and S. Sodin.  Non-backtracking random walks mix faster.  _Communications in Contemporary Mathematics_ , 9(04):585–603, 2007.
  * Angelini et al. [2015] M. C. Angelini, F. Caltagirone, F. Krzakala, and L. Zdeborová.  Spectral detection on sparse hypergraphs.  In _2015 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton)_ , pages 66–73. IEEE, 2015.
  * Bass [1992] H. Bass.  The Ihara-Selberg zeta function of a tree lattice.  _International Journal of Mathematics_ , 3(06):717–797, 1992.
  * Battiston et al. [2020] F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri.  Networks beyond pairwise interactions: Structure and dynamics.  _Physics Reports_ , 874:1–92, Aug. 2020.  ISSN 03701573.  10.1016/j.physrep.2020.05.004.
  * Battiston et al. [2021] F. Battiston, E. Amico, A. Barrat, G. Bianconi, G. Ferraz de Arruda, B. Franceschiello, I. Iacopini, S. Kéfi, V. Latora, Y. Moreno, et al.  The physics of higher-order interactions in complex systems.  _Nature Physics_ , 17(10):1093–1098, 2021.
  * Benson et al. [2018] A. R. Benson, R. Abebe, M. T. Schaub, A. Jadbabaie, and J. Kleinberg.  Simplicial closure and higher-order link prediction.  _Proceedings of the National Academy of Sciences_ , 115(48):E11221–E11230, Nov. 2018.  ISSN 0027-8424, 1091-6490.  10.1073/pnas.1800683115.
  * Bezanson et al. [2017] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah.  Julia: A fresh approach to numerical computing.  _SIAM Review_ , 59(1):65–98, 2017.
  * Bick et al. [2021] C. Bick, E. Gross, H. A. Harrington, and M. T. Schaub.  What are higher-order networks?  _arXiv:2104.11329 [nlin, stat]_ , Apr. 2021.
  * Bishop [2006] C. M. Bishop.  _Pattern Recognition and Machine Learning_.  Information Science and Statistics. Springer, New York, 2006.  ISBN 978-0-387-31073-2.
  * Blondel et al. [2008] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre.  Fast unfolding of communities in large networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2008(10):P10008, 2008.
  * Bojanek et al. [2020] K. Bojanek, Y. Zhu, and J. MacLean.  Cyclic transitions between higher order motifs underlie sustained asynchronous spiking in sparse recurrent networks.  _PLoS Computational Biology_ , 16(9):e1007409, 2020.
  * Bordenave et al. [2015] C. Bordenave, M. Lelarge, and L. Massoulié.  Non-backtracking spectrum of random graphs: Community detection and non-regular ramanujan graphs.  In _2015 IEEE 56th Annual Symposium on Foundations of Computer Science_ , pages 1347–1357. IEEE, 2015.
  * Caplar et al. [2017] N. Caplar, S. Tacchella, and S. Birrer.  Quantitative evaluation of gender bias in astronomical publications from citation counts.  _Nature Astronomy_ , 1(6):1–5, 2017.
  * Chang et al. [2020] J. Chang, Y. Chen, L. Qi, and H. Yan.  Hypergraph clustering using a new laplacian tensor with applications in image processing.  _SIAM Journal on Imaging Sciences_ , 13(3):1157–1178, 2020.
  * Chodrow [2020] P. S. Chodrow.  Configuration models of random hypergraphs.  _Journal of Complex Networks_ , 8(3):cnaa018, 2020\.
  * Chodrow et al. [2021] P. S. Chodrow, N. Veldt, and A. R. Benson.  Generative hypergraph clustering: From blockmodels to modularity.  _Science Advances_ , 7:eabh1303, 2021.
  * Contisciani et al. [2022] M. Contisciani, F. Battiston, and C. De Bacco.  Principled inference of hyperedges and overlapping communities in hypergraphs.  _arXiv:2204.05646 [physics, stat]_ , Apr. 2022.
  * Decelle et al. [2011] A. Decelle, F. Krzakala, C. Moore, and L. Zdeborová.  Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications.  _Physical Review E_ , 84(6):066106, 2011.
  * Dempster et al. [1977] A. P. Dempster, N. M. Laird, and D. B. Rubin.  Maximum likelihood from incomplete data via the EM algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Dion et al. [2018] M. L. Dion, J. L. Sumner, and S. M. Mitchell.  Gendered citation patterns across political science and social science methodology fields.  _Political Analysis_ , 26(3):312–327, 2018.
  * Dumitriu et al. [2021] I. Dumitriu, H. Wang, and Y. Zhu.  Partial recovery and weak consistency in the non-uniform hypergraph Stochastic Block Model.  _arXiv:2112.11671 [math, stat]_ , Dec. 2021.
  * Dworkin et al. [2020a] J. Dworkin, P. Zurn, and D. S. Bassett.  (in) citing action to realize an equitable future.  _Neuron_ , 106(6):890–894, 2020a.
  * Dworkin et al. [2020b] J. D. Dworkin, K. A. Linn, E. G. Teich, P. Zurn, R. T. Shinohara, and D. S. Bassett.  The extent and drivers of gender imbalance in neuroscience reference lists.  _Nature Neuroscience_ , 23(8):918–926, 2020b.
  * Fowler [2006a] J. H. Fowler.  Connecting the Congress: A study of cosponsorship networks.  _Political Analysis_ , 14(4):456–487, 2006a.
  * Fowler [2006b] J. H. Fowler.  Legislative cosponsorship networks in the US House and Senate.  _Social Networks_ , 28(4):454–465, 2006b.
  * Galuppi et al. [2021] F. Galuppi, R. Mulas, and L. Venturello.  Spectral theory of weighted hypergraphs via tensors.  _arXiv:2106.00277 [math]_ , June 2021.
  * Ghoshdastidar and Dukkipati [2017] D. Ghoshdastidar and A. Dukkipati.  Consistency of spectral hypergraph partitioning under planted partition model.  _The Annals of Statistics_ , 45(1), Feb. 2017.  ISSN 0090-5364.  10.1214/16-AOS1453.
  * Hu and Wang [2022] J. Hu and M. Wang.  Multiway Spherical Clustering via Degree-Corrected Tensor Block Models.  _arXiv:2201.07401 [math, stat]_ , Jan. 2022.
  * Ihara [1966] Y. Ihara.  On discrete subgroups of the two by two projective linear group over p-adic fields.  _Journal of the Mathematical Society of Japan_ , 18(3):219–235, 1966.
  * Jost and Mulas [2021] J. Jost and R. Mulas.  Normalized Laplace Operators for Hypergraphs with Real Coefficients.  _Journal of Complex Networks_ , 9(1):cnab009, Apr. 2021.  ISSN 2051-1310, 2051-1329.  10.1093/comnet/cnab009.
  * Kawamoto [2018] T. Kawamoto.  Algorithmic detectability threshold of the stochastic block model.  _Physical Review E_ , 97(3):032301, Mar. 2018\.  ISSN 2470-0045, 2470-0053.  10.1103/PhysRevE.97.032301.
  * Ke et al. [2019] Z. T. Ke, F. Shi, and D. Xia.  Community detection for hypergraph networks via regularized tensor power iteration.  _arXiv preprint arXiv:1909.06503_ , 2019.
  * Kempton [2016] M. C. Kempton.  Non-Backtracking Random Walks and a Weighted Ihara’s Theorem.  _Open Journal of Discrete Mathematics_ , 2016.
  * Kirkley et al. [2021] A. Kirkley, G. T. Cantwell, and M. E. J. Newman.  Belief propagation for networks with loops.  _Science Advances_ , 7(17):eabf1211, Apr. 2021\.  ISSN 2375-2548.  10.1126/sciadv.abf1211.
  * Krzakala et al. [2013] F. Krzakala, C. Moore, E. Mossel, J. Neeman, A. Sly, L. Zdeborová, and P. Zhang.  Spectral redemption in clustering sparse networks.  _Proceedings of the National Academy of Sciences_ , 110(52):20935–20940, 2013.
  * Lei and Rinaldo [2015] J. Lei and A. Rinaldo.  Consistency of spectral clustering in stochastic block models.  _The Annals of Statistics_ , 43(1):215–237, 2015\.
  * Lim [2005] L.-H. Lim.  Singular values and eigenvalues of tensors: A variational approach.  In _1st IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2005._ , pages 129–132. IEEE, 2005\.
  * Llorens et al. [2021] A. Llorens, A. Tzovara, L. Bellier, I. Bhaya-Grossman, A. Bidet-Caulet, W. K. Chang, Z. R. Cross, R. Dominguez-Faus, A. Flinker, Y. Fonken, et al.  Gender bias in academia: A lifetime problem that needs solutions.  _Neuron_ , 109(13):2047–2074, 2021.
  * MacQueen et al. [1967] J. MacQueen et al.  Some methods for classification and analysis of multivariate observations.  In _Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability_ , volume 1, pages 281–297. Oakland, CA, USA, 1967.
  * Maliniak et al. [2013] D. Maliniak, R. Powers, and B. F. Walter.  The gender citation gap in international relations.  _International Organization_ , 67(4):889–922, 2013.
  * Martin et al. [2014] T. Martin, X. Zhang, and M. E. Newman.  Localization and centrality in networks.  _Physical Review E_ , 90(5):052808, 2014.
  * Massoulié [2014] L. Massoulié.  Community detection thresholds and the weak Ramanujan property.  In _Proceedings of the Forty-Sixth Annual ACM Symposium on Theory of Computing_ , pages 694–703, 2014.
  * Mastrandrea et al. [2015] R. Mastrandrea, J. Fournet, and A. Barrat.  Contact Patterns in a High School: A Comparison between Data Collected Using Wearable Sensors, Contact Diaries and Friendship Surveys.  _PLoS ONE_ , 10(9):e0136497, Sept. 2015.  ISSN 1932-6203.  10.1371/journal.pone.0136497.
  * Mellor and Grusovin [2019] A. Mellor and A. Grusovin.  Graph comparison via the nonbacktracking spectrum.  _Physical Review E_ , 99(5):052309, 2019.
  * Mossel et al. [2015] E. Mossel, J. Neeman, and A. Sly.  Reconstruction and estimation in the planted partition model.  _Probability Theory and Related Fields_ , 162(3):431–461, 2015.
  * Mossel et al. [2018] E. Mossel, J. Neeman, and A. Sly.  A proof of the block model threshold conjecture.  _Combinatorica. An International Journal on Combinatorics and the Theory of Computing_ , 38(3):665–708, 2018.
  * Mulas and Zhang [2021] R. Mulas and D. Zhang.  Spectral Theory of Laplace Operators on Oriented Hypergraphs.  _Discrete Mathematics_ , 344(6):112372, June 2021\.  ISSN 0012365X.  10.1016/j.disc.2021.112372.
  * Mulas et al. [2021] R. Mulas, C. Kuehn, T. Böhle, and J. Jost.  Random walks and Laplacians on hypergraphs: When do they match?  _arXiv:2106.11663 [math]_ , June 2021.
  * Nadakuditi and Newman [2012] R. R. Nadakuditi and M. E. Newman.  Graph spectra and the detectability of community structure in networks.  _Physical Review Letters_ , 108(18):188701, 2012\.
  * Newman [2006] M. E. Newman.  Modularity and community structure in networks.  _Proceedings of the National Academy of Sciences_ , 103(23):8577–8582, 2006.
  * Peel et al. [2017] L. Peel, D. B. Larremore, and A. Clauset.  The ground truth about metadata and community detection in networks.  _Science Advances_ , 3(5):e1602548, May 2017.  ISSN 2375-2548.  10.1126/sciadv.1602548.
  * Porter et al. [2009] M. A. Porter, J.-P. Onnela, and P. J. Mucha.  Communities in networks.  _Notices of the AMS_ , 56(9):1082–1097, 2009\.
  * R Core Team [2022] R Core Team.  _R: A Language and Environment for Statistical Computing_.  R Foundation for Statistical Computing, Vienna, Austria, 2022.
  * Ramani et al. [2019] A. S. Ramani, N. Eikmeier, and D. F. Gleich.  Coin-flipping, ball-dropping, and grass-hopping for generating random graphs from matrices of edge probabilities.  _SIAM Review_ , 61(3):549–595, 2019.
  * Saad [2011] Y. Saad.  _Numerical methods for large eigenvalue problems: revised edition_.  SIAM, 2011.
  * Sebestyen [1962] G. S. Sebestyen.  _Decision-Making Processes in Pattern Recognition (ACM Monograph Series)_.  Macmillan Publishing Co., Inc., 1962.
  * Shi and Malik [2000] J. Shi and J. Malik.  Normalized cuts and image segmentation.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 22(8):888–905, 2000.
  * Stehlé et al. [2011] J. Stehlé, N. Voirin, A. Barrat, C. Cattuto, L. Isella, J.-F. Pinton, M. Quaggiotto, W. Van den Broeck, C. Régis, B. Lina, and P. Vanhems.  High-Resolution Measurements of Face-to-Face Contact Patterns in a Primary School.  _PLoS ONE_ , 6(8):e23176, Aug. 2011.  ISSN 1932-6203.  10.1371/journal.pone.0023176.
  * Stephan and Zhu [2022] L. Stephan and Y. Zhu.  Sparse random hypergraphs: Non-backtracking spectra and community detection.  _arXiv:2203.07346 [math, stat]_ , Mar. 2022.
  * Storm [2006] C. K. Storm.  The Zeta Function of a Hypergraph.  _The Electronic Journal of Combinatorics_ , 13(1):R84, Oct. 2006.  ISSN 1077-8926.  10.37236/1110.
  * Teich et al. [2021] E. G. Teich, J. Z. Kim, C. W. Lynn, S. C. Simon, A. A. Klishin, K. P. Szymula, P. Srivastava, L. C. Bassett, P. Zurn, J. D. Dworkin, et al.  Citation inequity and gendered citation practices in contemporary physics.  _arXiv preprint arXiv:2112.09047_ , 2021.
  * Torres et al. [2019] L. Torres, P. Suárez-Serrato, and T. Eliassi-Rad.  Non-backtracking cycles: Length spectrum theory and graph mining applications.  _Applied Network Science_ , 4(1):41, Dec. 2019\.  ISSN 2364-8228.  10.1007/s41109-019-0147-y.
  * Torres et al. [2021] L. Torres, K. S. Chan, H. Tong, and T. Eliassi-Rad.  Nonbacktracking Eigenvalues under Node Removal: X-Centrality and Targeted Immunization.  _SIAM Journal on Mathematics of Data Science_ , 3(2):656–675, Jan. 2021.  ISSN 2577-0187.  10.1137/20M1352132.
  * Von Luxburg [2007] U. Von Luxburg.  A tutorial on spectral clustering.  _Statistics and Computing_ , 17(4):395–416, 2007\.
  * Von Luxburg et al. [2008] U. Von Luxburg, M. Belkin, and O. Bousquet.  Consistency of spectral clustering.  _The Annals of Statistics_ , pages 555–586, 2008.
  * Wang et al. [2021] X. Wang, J. D. Dworkin, D. Zhou, J. Stiso, E. B. Falk, D. S. Bassett, P. Zurn, and D. M. Lydon-Staley.  Gendered citation practices in the field of communication.  _Annals of the International Communication Association_ , 45(2):134–153, 2021.
  * Wickham [2016] H. Wickham.  _ggplot2: Elegant Graphics for Data Analysis_.  Springer-Verlag New York, 2016.  ISBN 978-3-319-24277-4.
  * Zhou et al. [2006] D. Zhou, J. Huang, and B. Schölkopf.  Learning with hypergraphs: Clustering, classification, and embedding.  _Advances in Neural Information Processing Systems_ , 19:1601–1608, 2006.
  * Zurn et al. [2020] P. Zurn, D. S. Bassett, and N. C. Rust.  The citation diversity statement: a practice of transparency, a way of life.  _Trends in Cognitive Sciences_ , 24(9):669–672, 2020.

##  Appendix A Proof of Theorem 3.3

Our proof approach extends Kempton’s proof of the Ihara-Bass formula for
graphs [35]. The same approach was used by Stephan and Zhu [61] for the case
of uniform hypergraphs.

For each $k$, define operators $\mathbf{S}_{k}\in
L(\vec{\mathcal{E}},\mathcal{N})$, $\mathbf{T}_{k}\in
L(\mathcal{N},\vec{\mathcal{E}})$, and $\mathbf{W}_{k}\in
L(\vec{\mathcal{E}})$ with entries:

$\displaystyle s_{k;jR,i}$ $\displaystyle\coloneqq\begin{cases}1&\quad i\in
R\setminus j\;,\;\lvert R\rvert=k\\\ 0&\quad\text{otherwise.}\end{cases}$
$\displaystyle t_{k;i,jR}$ $\displaystyle\coloneqq\begin{cases}1&\quad
i=j\;,\;\lvert R\rvert=k\\\ 0&\quad\text{otherwise.}\end{cases}$
$\displaystyle w_{k;iQ,jR}$ $\displaystyle\coloneqq\begin{cases}1&\quad iQ\sim
jR\;,\;i\neq j\;,\;\lvert Q\rvert=\lvert R\rvert=k\\\
0&\quad\text{otherwise.}\end{cases}$

These operators satisfy several important relations. We begin with entrywise
calculations:

(27) $\displaystyle[\mathbf{T}_{k}\mathbf{S}_{k^{\prime}}]_{i,j}$
$\displaystyle=\sum_{\ell R\in\vec{\mathcal{E}}}t_{k;i,\ell
R}s_{k^{\prime};\ell
R,j}=\delta_{k,k^{\prime}}\lvert\left\\{R\in\mathcal{E}_{k}\;:\;i,j\in
R\right\\}\rvert\coloneqq\delta_{k,k^{\prime}}a_{k;i,j}$
$\displaystyle[\mathbf{T}_{k}\mathbf{W}_{\bar{k}}\mathbf{S}_{k^{\prime}}]_{i,j}$
$\displaystyle=\sum_{iQ,jR\in\vec{\mathcal{E}}}t_{k;i,iQ}w_{k;iQ,jR}s_{k^{\prime};jR,j}$
$\displaystyle=\delta_{k,k^{\prime}}\begin{cases}(k-1)d_{k;i,i}&\quad i=j\\\
(k-2)a_{k;i,j}&\quad\text{otherwise.}\end{cases}$ (28)
$\displaystyle=\delta_{k,k^{\prime}}\left[(k-1)\mathbf{D}_{k}+(k-2)\mathbf{A}_{k}\right]_{i,j}\;.$
$\displaystyle[\mathbf{S}_{k}\mathbf{T}_{k^{\prime}}-\delta_{k,k^{\prime}}\mathbf{W}_{k}]_{iQ,jR}$
$\displaystyle=\sum_{i\in\mathcal{N}}s_{k;iQ,i}t_{k^{\prime};i,jR}-\delta_{k,k^{\prime}}w_{k;iQ,jR}$
$\displaystyle=\begin{cases}1&\quad iQ\rightarrow jR\;,\;\lvert
iQ\rvert=k\;,\;\lvert jR\rvert=k^{\prime}\\\
0&\quad\text{otherwise}\end{cases}$ (29) $\displaystyle=b_{k\rightarrow
k^{\prime};iQ,jR}\;.$

Define block matrices

$\displaystyle\mathbf{S}\coloneqq\left[\begin{matrix}\mathbf{S}_{2}&&\\\
&\ddots&\\\
&&\mathbf{S}_{\bar{k}}\end{matrix}\right]\quad,\quad\mathbf{T}\coloneqq\left[\begin{matrix}\mathbf{T}_{2}&\cdots&\mathbf{T}_{\bar{k}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{T}_{2}&\cdots&\mathbf{T}_{\bar{k}}\end{matrix}\right]\quad\text{,
and}\quad\mathbf{W}\coloneqq\left[\begin{matrix}\mathbf{W}_{2}&&\\\
&\ddots&\\\ &&\mathbf{W}_{\bar{k}}\end{matrix}\right]\;.$

Direct multiplication and use of Eqs. 27, 28, and 29 gives the relations

(30) $\displaystyle\mathbf{T}\mathbf{S}$ $\displaystyle=\mathbf{A}$ (31)
$\displaystyle\mathbf{T}\mathbf{W}\mathbf{S}$
$\displaystyle=((\mathbf{K}-\mathbf{I}_{\kappa})\otimes\mathbf{I}_{n})\mathbf{D}+((\mathbf{K}-2\mathbf{I}_{\kappa})\otimes\mathbf{I}_{n})\mathbf{A}$
(32) $\displaystyle\mathbf{S}\mathbf{T}-\mathbf{W}$ $\displaystyle=\mathbf{B}$

We are now prepared for the main computation. The _push-through identity_
states that

(33)
$\displaystyle\det(\mathbf{X}+\mathbf{Y}\mathbf{Z})=\det(\mathbf{X})\det(\mathbf{I}+\mathbf{Z}\mathbf{X}^{-1}\mathbf{Y})\;,$

provided that $\mathbf{X}$ is invertible and all matrix products are well-
defined. Kempton [35] provides an elementary proof. Using Eq. 32 to write
$\mathbf{I}-\mu\mathbf{B}=\mathbf{I}-\mu\mathbf{S}\mathbf{T}+\mu\mathbf{W}$
and applying Eq. 33 to the righthand side gives

$\displaystyle\det(\mathbf{I}-\mu\mathbf{B})=\det(\mathbf{I}+\mu\mathbf{W})\det(\mathbf{I}-\mu\mathbf{T}(\mathbf{I}+\mu\mathbf{W})^{-1}\mathbf{S})\;.$

Focusing on the second factor, we compute

(34)
$\displaystyle(\mathbf{I}+\mu\mathbf{W})^{-1}=\left[\begin{matrix}\mathbf{P}_{2}(\mu)&&\\\
&\ddots&\\\ &&\mathbf{P}_{\bar{k}}(\mu)\;,\end{matrix}\right]$

where we have defined

$\displaystyle\mathbf{P}_{k}(\mu)\coloneqq
p_{k}(\mu)\mathbf{I}+q_{k}(\mu)\mathbf{W}_{k}$

with coefficients

$\displaystyle p_{k}(\mu)\coloneqq\frac{1+\mu(k-2)}{h_{k}(\mu)}\;,$
$\displaystyle\quad q_{k}(\mu)\coloneqq\frac{-\mu}{h_{k}(\mu)}\;,$
$\displaystyle\text{and}\;\quad h_{k}(\mu)\coloneqq(1-\mu)(1+\mu(k-1))\;.$

The derivation of this inverse uses the fact that $\mathbf{I}+\mu\mathbf{W}$
is a block-diagonal matrix with one block for each edge. Each block for an
edge of size $k$ has the form $\mathbf{I}+(\mu-1)\mathbf{E}$, where
$\mathbf{E}$ is a $k\times k$ matrix of ones.

Using Eqs. 30, 31, and 34, we now compute

$\displaystyle\mathbf{T}(\mathbf{I}+\mu\mathbf{W})^{-1}\mathbf{S}$
$\displaystyle=\left[\begin{matrix}p_{2}(\mu)\mathbf{A}_{2}&\cdots&p_{\bar{k}}(\mu)\mathbf{A}_{\bar{k}}\\\
\vdots&\ddots&\vdots\\\
p_{2}(\mu)\mathbf{A}_{2}&\cdots&p_{\bar{k}}(\mu)\mathbf{A}_{\bar{k}}\end{matrix}\right]$
$\displaystyle\quad+\left[\begin{matrix}q_{2}(\mu)\mathbf{D}_{2}&\cdots&q_{\bar{k}}(\mu)\left((K-1)\mathbf{D}_{K}+(K-2)\mathbf{A}_{K}\right)\\\
\vdots&\ddots&\vdots\\\
q_{2}(\mu)\mathbf{D}_{2}&\cdots&q_{\bar{k}}(\mu)\left((K-1)\mathbf{D}_{K}+(K-2)\mathbf{A}_{K}\right)\end{matrix}\right]$

Performing $n$ row multiplications by $h_{k}(\mu)$ for each $k$ yields

$\displaystyle\det(\mathbf{I}-\mu\mathbf{T}(\mathbf{I}+\mu\mathbf{W})^{-1}\mathbf{S})=\left(\prod_{k}h_{k}(\mu)^{-n}\right)\det\mathbf{M}(\mu)\;,$

where $\mathbf{M}(\mu)$ is the matrix

$\displaystyle\mathbf{M}(\mu)$
$\displaystyle\coloneqq\left[\begin{matrix}h_{2}(\mu)\mathbf{I}_{n}&&\\\
&\ddots&\\\
&&h_{\bar{k}}(\mu)\mathbf{I}_{n}\end{matrix}\right]-\mu\left[\begin{matrix}(1-\mu)\mathbf{A}_{2}&\cdots&(1+\mu(K-2))\mathbf{A}_{\bar{k}}\\\
\vdots&\ddots&\vdots\\\
(1-\mu)\mathbf{A}_{2}&\cdots&(1+\mu(K-2))\mathbf{A}_{\bar{k}}\end{matrix}\right]$
$\displaystyle\quad+\mu^{2}\left[\begin{matrix}\mathbf{D}_{2}&\cdots&(K-1)\mathbf{D}_{\bar{k}}+(K-2)\mathbf{A}_{\bar{k}}\\\
\vdots&\ddots&\vdots\\\
\mathbf{D}_{2}&\cdots&(K-1)\mathbf{D}_{\bar{k}}+(K-2)\mathbf{A}_{\bar{k}}\end{matrix}\right]$
$\displaystyle=\mathbf{I}_{\kappa}+\mu((\mathbf{K}-2\mathbf{I}_{\kappa})\otimes\mathbf{I}_{n}-\mathbf{A})+\mu^{2}((\mathbf{K}-\mathbf{I}_{\kappa})\otimes\mathbf{I}_{n})(\mathbf{D}-\mathbf{I}_{\kappa
n})\;.$

This gives second factor in the statement of Theorem 3.3, so our final step is
to address the factor $\det(\mathbf{I}+\mu\mathbf{W})$. We have

$\displaystyle\det(\mathbf{I}+\mu\mathbf{W})$
$\displaystyle=\prod_{k}(1-\mu)^{m_{k}(k-1)}(1+\mu(k-1))^{m_{k}}\;.$

We find in turn

$\displaystyle\frac{\det(\mathbf{I}+\mu\mathbf{W})}{\prod_{k}h_{k}(\mu)^{n}}=\prod_{k}(1-\mu)^{m_{k}(k-1)-n}(1+\mu(k-1))^{m_{k}-n}=f_{\mathcal{H}}(\mu)\;.$

This completes the computation and the proof.

###  A.1 Proof of 1

Recall that $\vec{m}$ is the total number of pointed edges. We make the
substitution $\mu=\frac{1}{\beta}$ and write

(35) $\displaystyle\beta^{-\vec{m}}p_{\mathbf{B}}(\beta)=\beta^{-2\kappa
n}f_{\mathcal{H}}(\beta^{-1})p_{\mathbf{B}^{\prime}}(\beta)\;.$

We can have $p_{\mathbf{B}}(\beta)=0$ only if either
$f_{\mathcal{H}}(\beta^{-1})=0$ or $p_{\mathbf{B}^{\prime}}(\beta)=0$. For
each $k$, if $m_{k}>n$, then $f_{\mathcal{H}}$ has $m_{k}-n$ roots of the form
$\beta=1-k$. Similarly, if $\sum_{k\in K}m_{k}(k-1)>\kappa n$, then
$f_{\mathcal{H}}$ has $\sum_{k\in K}m_{k}(k-1)-\kappa n$ roots of the form
$\beta=1$. These are the only roots of $f_{\mathcal{H}}$, and any remaining
roots of $p_{\mathbf{B}}$ must therefore be roots of
$p_{\mathbf{B}^{\prime}}$. If on the other hand
$p_{\mathbf{B}^{\prime}}(\beta^{\prime})=0$, then either
$p_{\mathbf{B}}(\beta)=0$ or $\beta^{-1}$ is a pole of $f_{\mathcal{H}}$. By
our factorization of $f_{\mathcal{H}}$, this can occur only if $\beta=1-k$ for
some $k$ or $\beta=1$. These cases can occur only if $m_{k}<n$ and $\sum_{k\in
K}m_{k}(k-1)<\kappa n$, respectively.

###  A.2 Proof of Lemma 3.4

We will establish two relationships between the vectors $\mathbf{x}_{1}$ and
$\mathbf{x}_{2}$. To establish the first relationship, begin by writing, for
any $i\in\mathcal{N}$ and $k\in K$,

$\displaystyle\beta
x_{1;k,i}=\beta\sum_{Q\in\mathcal{E}_{k}(i)}u_{iQ}=\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{R\in\mathcal{E}(i)\setminus
Q}\sum_{j\in R\setminus i}u_{jR}\;.$

The first equality is the definition of $\mathbf{x}_{1}$, while the second is
the eigenvector relation $\beta\mathbf{u}=\mathbf{B}\mathbf{u}$. Manipulating
the two inner sums gives

$\displaystyle\beta x_{1;k,i}$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\left[\sum_{R\in\mathcal{E}(i)}\sum_{j\in
R\setminus i}u_{jR}-\sum_{j\in Q\setminus i}u_{jQ}\right]$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\left[\sum_{k^{\prime}\in
K}\sum_{R\in\mathcal{E}_{k^{\prime}}(i)}\sum_{j\in R\setminus
i}u_{jR}-\sum_{j\in Q\setminus i}u_{jQ}\right]$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\left[\sum_{k^{\prime}\in
K}x_{2;k^{\prime},i}-\sum_{j\in Q\setminus i}u_{jQ}\right]$ (36)
$\displaystyle=d_{ik}\sum_{k^{\prime}\in K}x_{2;k^{\prime},i}-x_{2;k,i}\;,$

where $d_{ik}$ is the number of $k$-edges incident to node $i$ and $x_{2;k,i}$
is as defined in Eq. 2.

To establish the second relationship, we similarly write

$\displaystyle\beta x_{2;k,i}=\beta\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in
Q\setminus i}u_{jQ}=\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in Q\setminus
i}\sum_{R\in\vec{\mathcal{E}}(i)\setminus iQ}\sum_{h\in R\setminus
j}u_{hR}\;,$

which again follows from the definition of $\mathbf{x}_{2}$ and the
eigenvector relation $\beta\mathbf{u}=\mathbf{B}\mathbf{u}$. We further
calculate

$\displaystyle\beta x_{2;k,i}$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in Q\setminus
i}\left[\sum_{R\in\vec{\mathcal{E}}(i)}\sum_{h\in R\setminus
j}u_{hR}-\sum_{h\in Q\setminus j}u_{hQ}\right]$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in Q\setminus
i}\left[\sum_{k^{\prime}\in
K}\sum_{R\in\vec{\mathcal{E}}_{k^{\prime}}(i)}\sum_{h\in R\setminus
j}u_{hR}-\sum_{h\in Q\setminus j}u_{hQ}\right]$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in Q\setminus
i}\left[\sum_{k^{\prime}\in K}x_{2;k^{\prime},j}-\sum_{h\in Q\setminus
j}u_{hQ}\right]\;.$

Focusing now on the second term in brackets, we compute

$\displaystyle\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in Q\setminus i}\sum_{h\in
Q\setminus j}u_{hQ}$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\left[\sum_{j\in Q}\sum_{h\in
Q\setminus j}u_{hQ}-\sum_{h\in Q\setminus i}u_{hQ}\right]$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\left[(k-1)\sum_{h\in
Q}u_{hQ}-\sum_{h\in Q\setminus i}u_{hQ}\right]$
$\displaystyle=\sum_{Q\in\mathcal{E}_{k}(i)}\left[(k-1)\left(\sum_{h\in
Q\setminus i}u_{hQ}+u_{iQ}\right)-\sum_{h\in Q\setminus i}u_{hQ}\right]$
$\displaystyle=(k-2)x_{2;k,i}+(k-1)x_{1;k,i}\;.$

This gives the relation

(37) $\displaystyle\beta x_{2;k,i}=\sum_{Q\in\mathcal{E}_{k}(i)}\sum_{j\in
Q\setminus i}\sum_{k^{\prime}\in
K}x_{2;k^{\prime},j}-(k-2)x_{2;k,i}-(k-1)x_{1;k,i}\;.$

When expressed in matrix form, Eqs. 36 and 37 read

(38) $\displaystyle\beta\left(\begin{matrix}\mathbf{x}_{1}\\\
\mathbf{x}_{2}\end{matrix}\right)=\mathbf{B}^{\prime}\left(\begin{matrix}\mathbf{x}_{1}\\\
\mathbf{x}_{2}\end{matrix}\right)\;.$

Thus, $\mathbf{v}=(\mathbf{x}_{1},\mathbf{x}_{2})^{T}$ is an eigenvector of
$\mathbf{B}^{\prime}$ unless $\mathbf{v}=\mathbf{0}$, as was to be shown.

##  Appendix B Precise Statement and Proof of 2

Let $\mathcal{Z}$ be the alphabet of cluster labels, with
$\lvert\mathcal{Z}\rvert=\ell$. Let $\mathcal{M}$ be the space of possible
messages $\boldsymbol{\mu}$; we can identify $\mathcal{M}$ with a product of
probability $(\ell-1)$-simplices containing one factor for each node-tuple
pair. Let $\bar{\boldsymbol{\mu}}$ be the vector of messages with entries
$\bar{\mu}_{iR}^{(s)}=q^{(s)}$ for all nodes $i$, subsets
$R\in\mathcal{R}(i)$, and labels $s\in\mathcal{Z}$.

We will consider perturbations to the belief-propagation dynamics. The set of
relevant perturbations has some structure. The normalization condition
$\sum_{s\in\mathcal{Z}}\mu_{iR}^{(s)}=1$ on elements of $\mathcal{M}$ requires
that perturbations $\boldsymbol{\epsilon}$ to a message vector
$\boldsymbol{\mu}$ must satisfy must satisfy $\sum_{s}\epsilon_{iR}^{(s)}=0$.
Letting $\boldsymbol{\Pi}$ denote the projection operator onto the subspace
defined by this relation, we have

(39)
$\displaystyle[\boldsymbol{\Pi}\boldsymbol{\epsilon}]_{iR}^{(s)}=\epsilon_{iR}^{(s)}-\frac{1}{\ell}\sum_{t\in\mathcal{Z}}\epsilon_{iR}^{(t)}\;.$

For a given hypergraph realization, we can separate $\boldsymbol{\mu}$ into
components $(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})$, where entries of
$\boldsymbol{\mu}_{0}$ correspond to unrealized edges and entries of
$\boldsymbol{\mu}_{1}$ correspond to realized edges. We can similarly separate
the components of the function $\mathbf{F}$. This allows us to write the BP
update dynamics in the form

$\displaystyle\boldsymbol{\mu}_{0}$
$\displaystyle\leftarrow\mathbf{F}_{0}(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})$
$\displaystyle\boldsymbol{\mu}_{1}$
$\displaystyle\leftarrow\mathbf{F}_{1}(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})\;.$

We are now prepared to state a precise analog to to the heuristic 2.

######  Theorem B.1.

Let $\mathcal{H}$ be sampled from the sparse Bernoulli blockmodel. Then, as
$n$ grows large,
$\mathbf{F}(\bar{\boldsymbol{\mu}})\doteq\bar{\boldsymbol{\mu}}$. Furthermore,

(40)
$\displaystyle\mathbf{J}_{10}\coloneqq\boldsymbol{\Pi}\frac{\partial\mathbf{F}_{1}(\bar{\boldsymbol{\mu}})}{\partial\boldsymbol{\mu}_{0}}\doteq\mathbf{0}\quad\text{and}\quad\mathbf{J}_{11}\coloneqq\boldsymbol{\Pi}\frac{\partial\mathbf{F}_{1}(\bar{\boldsymbol{\mu}})}{\partial\boldsymbol{\mu}_{1}}\doteq\sum_{k\in
K}(\mathbf{G}_{k}\otimes\mathbf{B}_{k}).$

######  Proof B.2.

Throughout this proof, sums indexed by $\mathbf{z}_{R}$ are assumed to run
through $\mathcal{Z}^{\lvert R\rvert}$.

We’ll first compute several approximations describing how messages propagate
along unrealized edges, i.e. subsets $R$ such that $a_{R}=0$. Since
$\eta(a_{R}=0|\mathbf{z}_{R})=1-\omega(\mathbf{z}_{R})n^{1-\lvert R\rvert}$,
Eq. 13 becomes

$\displaystyle\nu_{Ri}^{(s)}$
$\displaystyle\leftarrow\frac{1}{Z_{Ri}}\sum_{\mathbf{z}_{R}:z_{i}=s}(1-\omega(\mathbf{z}_{R})n^{1-\lvert
R\rvert})\prod_{j\in R\setminus i}\mu_{jR}^{(z_{j})}$
$\displaystyle=\left(1-O(n^{1-\lvert
R\rvert})\right)\frac{1}{Z_{Ri}}\sum_{\mathbf{z}_{R}:z_{i}=s}\prod_{j\in
R\setminus i}\mu_{jR}^{(z_{j})}$ (41) $\displaystyle=\left(1-O(n^{1-\lvert
R\rvert})\right)\frac{1}{Z_{Ri}}\;.$

The last line follows from the normalization of the messages $\mu_{jR}^{(s)}$,
since the sum ranges over all possible labelings of the nodes in $R\setminus
i$. Since the smallest possible edge size is $k=2$, we have shown that
$\nu_{Ri}^{(s)}\doteq Z_{Ri}^{-1}$. In particular, $\nu_{Ri}^{(s)}$ is
approximately constant with respect to $s$.

We can also approximate $\mu_{iR}^{(s)}$ in the case $a_{R}=0$, using Eqs. 41,
12, and 14 to obtain

$\displaystyle\mu_{iR}^{(s)}$
$\displaystyle\leftarrow\frac{1}{Z_{iR}}q^{(s)}\prod_{Q\in\mathcal{R}(i)\setminus
R}\nu_{Qi}^{(s)}$ $\displaystyle=\left(1+O(n^{1-\lvert
R\rvert})\right)\frac{Z_{Ri}}{Z_{iR}}q^{(s)}\prod_{Q\in\mathcal{R}(i)}\nu_{Qi}^{(s)}$
(42) $\displaystyle\doteq\mu_{i}^{(s)}\;.$

Here, we are able to identify the normalizing constant
$Z_{i}=\frac{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 $\mathbf{F}_{1}$.
Substituting Eq. 13 into Eq. 12 and absorbing normalizing constants allows us
to eliminate the messages $\nu_{Ri}^{(s)}$ entirely, obtaining an explicit
form for $\mathbf{F}_{1}$:

(43)
$\displaystyle\mathbf{F}_{1}(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})_{iR}^{(s)}=\frac{1}{Z_{iR}}q^{(s)}\prod_{Q\in\mathcal{R}(i)\setminus
R}\sum_{\mathbf{z}_{Q}:z_{i}=s}\eta(a_{Q}|\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{jQ}^{(z_{j})}\;.$

The updates in $\boldsymbol{\mu}_{1}$ are now
$\boldsymbol{\mu}_{1}\leftarrow\mathbf{F}_{1}(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})$.

Let us write $\mathbf{F}_{1}$ in the form

$\displaystyle\mathbf{F}_{1}(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})_{iR}^{(s)}=\frac{1}{Z_{iR}}q^{(s)}M_{iR}^{(s)}N_{iR}^{(s)}\;,$

where $M_{iR}^{(s)}$ contains factors corresponding to sets $Q$ such that
$a_{Q}=1$, while $N_{iR}^{(s)}$ contains factors for sets $Q$ such that
$a_{Q}=0$. We can expand $N_{iR}^{(s)}$:

$\displaystyle N_{iR}^{(s)}$
$\displaystyle=\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus R\\\
a_{Q}=0\end{subarray}}\sum_{\mathbf{z}_{Q}:z_{i}=s}\eta(a_{Q}=0|\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{jQ}^{(z_{j})}$
$\displaystyle=\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus R\\\
a_{Q}=0\end{subarray}}\sum_{\mathbf{z}_{Q}:z_{i}=s}(1-\omega(\mathbf{z}_{Q})n^{1-\lvert
Q\rvert})\prod_{j\in Q\setminus i}\mu_{jQ}^{(z_{j})}$
$\displaystyle=\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus R\\\
a_{Q}=0\end{subarray}}\left(1-n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{jQ}^{(z_{j})}\right)\;,$

where we have used the normalization of the messages $\mu_{jQ}^{(z_{j})}$ in
the last line.

We will now approximate $N_{iR}^{(s)}$ 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_{iR}^{(s)}$, we approximate
$\mu_{jQ}^{(z_{j})}=(1+O(n^{1-\lvert Q\rvert}))\mu_{j}^{(z_{j})}$. Next,

$\displaystyle N_{iR}^{(s)}$
$\displaystyle=\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus R\\\
a_{Q}=0\end{subarray}}\left(1-n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}(1+O(n^{1-\lvert Q\rvert}))\mu_{j}^{(z_{j})}\right)$
$\displaystyle\doteq(1+O(n^{-1}))\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus
R\\\ a_{Q}=0\end{subarray}}\left(1-n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{j}^{(z_{j})}\right)$
$\displaystyle=(1+O(n^{-1}))\frac{\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\end{subarray}}\left(1-n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus
i}\mu_{j}^{(z_{j})}\right)}{\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\\\
a_{Q}=1\end{subarray}}\left(1-n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{j}^{(z_{j})}\right)}\;.$

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_{k}^{(s)}$. 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,

$\displaystyle N_{iR}^{(s)}$
$\displaystyle\doteq\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\end{subarray}}\left(1-n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{j}^{(z_{j})}\right)$ $\displaystyle\coloneqq N^{(s)}\;,$

where we have defined

(44) $\displaystyle
N^{(s)}\coloneqq\prod_{k}\prod_{\begin{subarray}{c}Q\in\mathcal{R}_{k}\end{subarray}}\left(1-n^{1-k}\sum_{\mathbf{z}_{Q}:z_{q_{1}}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus q_{1}}\mu_{j}^{(z_{j})}\right)$

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 Eq. 43
satisfies

$\displaystyle\mathbf{F}_{1}(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})_{iR}^{(s)}$
$\displaystyle\doteq\frac{N^{(s)}q^{(s)}}{Z_{iR}}\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus
R\\\ a_{Q}=1\end{subarray}}n^{1-\lvert
Q\rvert}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{jQ}^{(z_{j})}\;.$ (45)
$\displaystyle=\frac{N^{(s)}q^{(s)}}{Z_{iR}}\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus
R\\\
a_{Q}=1\end{subarray}}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\mu_{jQ}^{(z_{j})}\;.$

Here we have absorbed factors that do not depend on $s$ into $Z_{iR}$.
Importantly, $\mathbf{F}_{1}$ depends on $\boldsymbol{\mu}_{0}$ only through
the field term $N^{(s)}$.

We next consider the behavior of Eq. 45 near the point
$\bar{\boldsymbol{\mu}}$ with entries $\bar{\mu}_{iR}^{(s)}=q^{(s)}$. Let
$\boldsymbol{\epsilon}$ be a perturbation vector assumed small. We’ll first
consider the field term. Let $\Delta$ refer to terms of order
$O(n^{-1}+\lVert\boldsymbol{\epsilon}\rVert)$. The perturbed field term
$N^{(s)}$ now reads

$\displaystyle
N^{(s)}(\bar{\boldsymbol{\mu}}+\boldsymbol{\epsilon})_{iR}^{(s)}$
$\displaystyle=\prod_{k}\prod_{\begin{subarray}{c}Q\in\mathcal{R}_{k}\end{subarray}}\left(1-n^{1-k}\sum_{\mathbf{z}_{Q}:z_{q_{1}}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus q_{1}}\left(q^{(z_{j})}+\epsilon_{jQ}^{(z_{j})}\right)\right)$
$\displaystyle\doteq\prod_{k}\prod_{\begin{subarray}{c}Q\in\mathcal{R}_{k}\end{subarray}}\left(1-n^{1-k}\sum_{\mathbf{z}_{Q}:z_{q_{1}}=s}\omega(\mathbf{z}_{Q})\left(\prod_{j\in
Q\setminus q_{1}}q^{(z_{j})}+\sum_{j\in Q\setminus
q_{1}}\epsilon_{jQ}^{(z_{j})}\prod_{\ell\in Q\setminus
q_{1},j}q^{(z_{\ell})}\right)\right)$
$\displaystyle=\prod_{k}\prod_{Q\in\mathcal{R}_{k}}\left(1-(k-1)!n^{1-k}c_{k}-n^{1-k}\sum_{\mathbf{z}_{Q}:z_{q_{1}}=s}\omega(\mathbf{z}_{Q})\sum_{j\in
Q\setminus q_{1}}\epsilon_{jQ}^{(z_{j})}\prod_{\ell\in Q\setminus
q_{1},j}q^{(z_{\ell})}\right)$
$\displaystyle=(1+\Delta)\prod_{k}\prod_{Q\in\mathcal{R}_{k}}\left(1-(k-1)!n^{1-k}c_{k}\right)$
$\displaystyle\coloneqq(1+\Delta)N\;,$

where we have defined $N$ to absorb the products. We have shown that, near
$\bar{\boldsymbol{\mu}}$, the field term $N^{(s)}$ approximately does not
depend on $s$.

Paralleling the partition
$\boldsymbol{\mu}=(\boldsymbol{\mu}_{0},\boldsymbol{\mu}_{1})$, we can
partition the entries of $\boldsymbol{\epsilon}$ as
$\boldsymbol{\epsilon}=(\boldsymbol{\epsilon}_{0},\boldsymbol{\epsilon}_{1})$,
again corresponding to unrealized and realized edges. From Eq. 45,

$\displaystyle\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}+\boldsymbol{\epsilon})_{iR}^{(s)}$
$\displaystyle=\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0}+\boldsymbol{\epsilon}_{0},\bar{\boldsymbol{\mu}}_{1}+\boldsymbol{\epsilon}_{1})_{iR}^{(s)}$
$\displaystyle=(1+\Delta)\frac{Nq^{(s)}}{Z_{iR}}\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus
R\\\
a_{Q}=1\end{subarray}}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\prod_{j\in
Q\setminus i}\left(q^{(z_{j})}+\epsilon_{jQ}^{(z_{j})}\right)$
$\displaystyle\doteq(1+\Delta)\frac{Nq^{(s)}}{Z_{iR}}\prod_{\begin{subarray}{c}Q\in\mathcal{R}(i)\setminus
R\\\
a_{Q}=1\end{subarray}}\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\left[\prod_{j\in
Q\setminus i}q^{(z_{j})}+\sum_{j\in Q\setminus
i}\epsilon_{jQ}^{(z_{j})}\prod_{h\in Q\setminus i,j}q^{(z_{h})}\right]$
$\displaystyle=(1+\Delta)\frac{Nq^{(s)}}{Z_{iR}}\prod_{k}\prod_{\begin{subarray}{c}Q\in\mathcal{R}_{k}(i)\setminus
R\\\
a_{Q}=1\end{subarray}}\left((k-1)!c_{k}+\sum_{\mathbf{z}_{Q}:z_{i}=s}\omega(\mathbf{z}_{Q})\sum_{j\in
Q\setminus i}\epsilon_{jQ}^{(z_{j})}\prod_{h\in Q\setminus
i,j}q^{(z_{h})}\right)\;.$

By conditioning on the label of $j$, we can simplify the second term in the
factor:

$\displaystyle\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}+\boldsymbol{\epsilon})_{iR}^{(s)}$
$\displaystyle=(1+\Delta)\frac{Nq^{(s)}}{Z_{iR}}\prod_{k}\prod_{\begin{subarray}{c}Q\in\mathcal{R}_{k}(i)\setminus
R\\\
a_{Q}=1\end{subarray}}\left((k-1)!c_{k}+(k-2)!\sum_{t\in\mathcal{Z}}c_{k}^{(s,t)}\sum_{j\in
Q\setminus i}\epsilon_{jQ}^{(t)}\right)$ (46)
$\displaystyle=(1+\Delta)\frac{Nq^{(s)}}{Z_{iR}}\prod_{k}\prod_{\begin{subarray}{c}Q\in\mathcal{R}_{k}(i)\setminus
R\\\
a_{Q}=1\end{subarray}}\left(1+\sum_{t\in\mathcal{Z}}\frac{c_{k}^{(s,t)}}{(k-1)c_{k}}\sum_{j\in
Q\setminus i}\epsilon_{jQ}^{(t)}\right)\;.$

When $\boldsymbol{\epsilon}=\mathbf{0}$, we have that
$\mathbf{F}(\bar{\boldsymbol{\mu}})_{iR}^{(s)}\doteq\frac{N}{Z_{iR}}q^{(s)}$.
Since the messages must normalize in $s$, we have that, up to errors that can
be absorbed into the first factor, $\frac{N}{Z_{iR}}=1$. This shows that

(47)
$\displaystyle\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1})$
$\displaystyle\doteq\bar{\boldsymbol{\mu}}_{1}\;.$

A similar calculation shows that

$\displaystyle\mathbf{F}_{0}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1})$
$\displaystyle\doteq\bar{\boldsymbol{\mu}}_{0}\;,$

which relations jointly give
$\mathbf{F}(\bar{\boldsymbol{\mu}})\doteq\bar{\boldsymbol{\mu}}$. This proves
the first clause of the theorem.

Furthermore, since $\mathbf{F}_{1}$ depends on $\bar{\boldsymbol{\mu}}_{0}$
only through the factor
$N/Z_{iR}=1+\Delta=1+O(n^{-1}+\lVert\boldsymbol{\epsilon}\rVert)$, we have
that any derivative of $\mathbf{F}_{1}$ in a direction corresponding to
$\boldsymbol{\mu}_{0}$ is of order $O(n^{-1})$. Thus, the Jacobian
$\frac{\partial\mathbf{F}_{1}}{\partial\boldsymbol{\mu}_{0}}$, evaluated at
$\bar{\boldsymbol{\mu}}$, has entries of order $O(n^{-1})$. The projected
matrix
$\mathbf{J}_{10}=\boldsymbol{\Pi}\frac{\partial\mathbf{F}_{1}}{\partial\boldsymbol{\mu}_{0}}$
is also of order $O(n^{-1})$, proving the first equation in Eq. 40.

It remains to compute $\mathbf{J}_{11}$ at $\bar{\boldsymbol{\mu}}$. Expanding
the product in Eq. 46 to first order in $\boldsymbol{\epsilon}$ and separating
the arguments of $\mathbf{F}_{1}$ gives

$\displaystyle(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1}+\boldsymbol{\epsilon}_{1})_{iR}^{(s)}=(1+\Delta)q^{(s)}\left(1+\sum_{k\in
K}\sum_{t\in\mathcal{Z}}\frac{c_{k}^{(s,t)}}{(k-1)c_{k}}\sum_{\begin{subarray}{c}Q\in\mathcal{R}_{k}(i)\setminus
R\\\ a_{Q}=1\end{subarray}}\sum_{j\in Q\setminus
i}\epsilon_{jQ}^{(t)}\right)\;.$

Using Eq. 47 gives

$\displaystyle\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1}+\boldsymbol{\epsilon}_{1})-\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1})_{iR}^{(s)}=(1+\Delta)q^{(s)}\sum_{k\in
K}\sum_{t\in\mathcal{Z}}\frac{c_{k}^{(s,t)}}{(k-1)c_{k}}\sum_{\begin{subarray}{c}Q\in\mathcal{R}_{k}(i)\setminus
R\\\ a_{Q}=1\end{subarray}}\sum_{j\in Q\setminus i}\epsilon_{jQ}^{(t)}\;.$

We now apply the projection $\boldsymbol{\Pi}$ onto the subspace of admissible
perturbations, yielding

$\displaystyle[\boldsymbol{\Pi}\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1}+\boldsymbol{\epsilon}_{1})-\boldsymbol{\Pi}\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1})]_{iR}^{(s)}=(1+\Delta)q^{(s)}\sum_{k\in
K}\sum_{t\in\mathcal{Z}}\left(\frac{c_{k}^{(s,t)}}{(k-1)c_{k}}-1\right)\sum_{\begin{subarray}{c}Q\in\mathcal{R}_{k}(i)\setminus
R\\\ a_{Q}=1\end{subarray}}\sum_{j\in Q\setminus i}\epsilon_{jQ}^{(t)}\;.$

We can identify $R$ with an edge $e$, and the pair $(i,R)$ as a pointed edge
$\vec{e}$ with $p(\vec{e})=i$. Doing the same for $Q$ and $j$, we can
recognize the two rightmost sums as the action of $\mathbf{B}_{k}$ on the
perturbation vector $\boldsymbol{\epsilon}$. Using the definition of
$\mathbf{G}_{k}$, we can write this relation as

$\displaystyle\boldsymbol{\Pi}\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1}+\boldsymbol{\epsilon}_{1})-\boldsymbol{\Pi}\mathbf{F}_{1}(\bar{\boldsymbol{\mu}}_{0},\bar{\boldsymbol{\mu}}_{1})=(1+\Delta)\sum_{k\in
K}(\mathbf{G}_{k}\otimes\mathbf{B}_{k})\boldsymbol{\epsilon}\;.$

Ignoring the error term, this relation would define the Jacobian
$\mathbf{J}_{11}\coloneqq\boldsymbol{\Pi}\frac{\partial\mathbf{F}_{1}}{\partial\boldsymbol{\mu}_{1}}$
as equal to the righthand side. Allowing
$\lVert\boldsymbol{\epsilon}\rVert\rightarrow 0$, we conclude that
$\mathbf{J}_{11}$ satisfies

$\displaystyle\mathbf{J}_{11}\doteq(1+O(n^{-1}))\sum_{k\in
K}(\mathbf{G}_{k}\otimes\mathbf{B}_{k})\;,$

which establishes the second clause of Eq. 40 and completes the proof.

##  Appendix C Proof of Theorem 4.1

######  Proof C.1.

We will prove Eq. 9; the proof of Eq. 10 is similar but somewhat shorter. Let
$\tilde{\mathbf{v}}$ be a vector indexed by tuples $Q\in\binom{[n]}{k}$ and
nodes $i\in Q$ with entries

$\displaystyle\tilde{v}_{iQ}$ $\displaystyle\coloneqq\begin{cases}v_{iQ}&\quad
Q\in\mathcal{E}\\\ 0&\quad\text{otherwise.}\end{cases}$

Let $\tilde{\mathbf{B}}_{k}$ be the matrix with entries

$\displaystyle\tilde{b}_{k;iQ,jR}=\begin{cases}b_{k;iQ,jR}&\quad
Q,R\in\mathcal{E}\\\ 0&\quad\text{otherwise.}\end{cases}$

We are going to show that
$\mathbb{E}\left[(\tilde{\mathbf{B}}_{k}\tilde{\mathbf{v}})_{iQ}-\beta_{k}\tilde{v}_{iQ}|Q\in\mathcal{E}\right]\doteq
0$; since $\tilde{\mathbf{B}}_{k}$ and $\tilde{\mathbf{v}}$ agree with
$\mathbf{B}_{k}$ and $\mathbf{v}$ conditioned on the event $Q\in\mathcal{E}$,
this will imply Eq. 10.

We proceed via direct computation. Expanding the expectation, we can write

$\displaystyle\mathbb{E}[(\tilde{\mathbf{B}}_{k}\tilde{\mathbf{v}})_{iQ}|Q\in\mathcal{E}]$
$\displaystyle=\sum_{jR\in\mathcal{E}_{k}}\mathbb{E}[\tilde{b}_{iQ,jR}\tilde{v}_{jR}|Q\in\mathcal{E}]$
$\displaystyle=\sum_{jR\in\mathcal{E}_{k}}\mathbb{E}[\tilde{b}_{iQ,jR}\tilde{v}_{jR}|Q,R\in\mathcal{E}]\eta(R\in\mathcal{E})$
$\displaystyle=\sum_{jR\in\mathcal{E}_{k}}\eta(R\in\mathcal{E})b_{iQ,jR}v_{jR}$
$\displaystyle=\sum_{jR\in\mathcal{E}_{k}}\eta(R\in\mathcal{E})b_{iQ,jR}\sum_{\ell\in
R\setminus j}\sigma_{\ell}\;.$

The third line follows from the fact that, conditioned on the event
$Q,R\in\mathcal{E}$, $\tilde{b}_{iQ,jR}=b_{iQ,jR}$ and
$\tilde{v}_{jR}=v_{jR}$. Proceeding from the fourth line, we can evaluate
$b_{iQ,jR}$ and rearrange the sums:

$\displaystyle\mathbb{E}[(\tilde{\mathbf{B}}_{k}\tilde{\mathbf{v}})_{iQ}|Q\in\mathcal{E}]$
$\displaystyle=\sum_{j\in Q\setminus
i}\sum_{\begin{subarray}{c}R\in\mathcal{R}_{k}(j)\setminus
Q\end{subarray}}\eta(R\in\mathcal{E})\sum_{\ell\in R\setminus j}\sigma_{\ell}$
$\displaystyle=\sum_{j\in Q\setminus
i}\sum_{\begin{subarray}{c}R^{\prime}\in\binom{[n]\setminus j}{k-1}\\\
R^{\prime}\neq Q\setminus j\end{subarray}}\eta(R^{\prime}\cup
j\in\mathcal{E})\sum_{\ell\in R^{\prime}}\sigma_{\ell}$
$\displaystyle=\sum_{j\in Q\setminus i}\sum_{\ell\neq
j}\sigma_{\ell}\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus\\{j,\ell\\}}{k-2}\\\
R^{\prime\prime}\neq
Q\setminus\\{\ell,j\\}\end{subarray}}\eta(R^{\prime\prime}\cup\\{j,\ell\\}\in\mathcal{E})\;.$

The inner sum satisfies

$\displaystyle\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus\\{j,\ell\\}}{k-2}\\\
R^{\prime\prime}\neq
Q\setminus\\{\ell,j\\}\end{subarray}}\eta(R^{\prime\prime}\cup\\{j,\ell\\}\in\mathcal{E})\doteq\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus\\{j,\ell\\}}{k-2}\end{subarray}}\eta(R^{\prime\prime}\cup\\{j,\ell\\}\in\mathcal{E})\;,$

The asymptotic equality holds because there are $\binom{n-2}{k-2}$ terms, of
which the condition $R^{{}^{\prime\prime}}\neq Q\setminus\\{\ell,j\\}$
excludes only one.111In the edge case $k=2$, the two sides are in fact exactly
equal.  We proceed to compute the sum on the righthand side.

$\displaystyle\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus\\{j,\ell\\}}{k-2}\end{subarray}}\eta(R^{\prime\prime}\cup\\{j,\ell\\}\in\mathcal{E})$
$\displaystyle=n^{1-k}\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus\\{j,\ell\\}}{k-2}\end{subarray}}\omega(\mathbf{z}_{R^{\prime\prime}},z_{j},z_{\ell})$
$\displaystyle=n^{1-k}\sum_{\mathbf{z}\in\mathcal{Z}^{k-2}}\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus\\{j,\ell\\}}{k-2}\\\
\mathbf{z}_{R^{\prime\prime}}=\mathbf{z}\end{subarray}}\omega(\mathbf{z},z_{j},z_{\ell})\;.$

We can make progress by counting the number of subsets $R^{\prime\prime}$ that
realize each specified label vector $\mathbf{z}$. There are $\binom{n-2}{k-2}$
possible choices, and the proportion of these choices satisfying
$\mathbf{z}_{R}=\mathbf{z}$ is asymptotically $\prod_{s\in\mathbf{z}}q^{(s)}$.
This gives

$\displaystyle\sum_{\begin{subarray}{c}R^{\prime\prime}\in\binom{[n]\setminus
j,\ell}{k-2}\end{subarray}}\eta(R^{\prime\prime}\cup\\{j,\ell\\}\in\mathcal{E})$
$\displaystyle\doteq
n^{1-k}\binom{n-2}{k-2}\sum_{\mathbf{z}\in\mathcal{Z}^{k-2}}\omega(\mathbf{z},z_{j},z_{\ell})\prod_{s\in\mathbf{z}}q^{(s)}$
$\displaystyle\doteq\frac{1}{n}\frac{1}{(k-2)!}\sum_{\mathbf{z}\in\mathcal{Z}^{k-2}}\omega(\mathbf{z},z_{j},z_{\ell})\prod_{s\in\mathbf{z}}q^{(s)}$
$\displaystyle=\frac{1}{n}c_{k}^{(z_{j},z_{\ell})}\;,$

where we have used Eq. 7 in the final line. We therefore have

(48)
$\displaystyle\mathbb{E}[(\tilde{\mathbf{B}}_{k}\tilde{\mathbf{v}})_{iQ}|Q\in\mathcal{E}]\doteq\frac{1}{n}\sum_{j\in
Q\setminus i}\sum_{\ell\neq j}\sigma_{\ell}c_{k}^{(z_{j},z_{\ell})}\;.$

Let us split this sum according to whether $z_{\ell}=z_{j}$:

$\displaystyle\mathbb{E}[(\tilde{\mathbf{B}}_{k}\tilde{\mathbf{v}})_{iQ}|Q\in\mathcal{E}]$
$\displaystyle\doteq\frac{1}{n}\sum_{j\in Q\setminus
i}\left(\sum_{\begin{subarray}{c}\ell\neq j\\\
z_{\ell}=z_{j}\end{subarray}}\sigma_{\ell}c_{k}^{(z_{j},z_{\ell})}+\sum_{\begin{subarray}{c}\ell\neq
j\\\ z_{\ell}\neq
z_{j}\end{subarray}}\sigma_{\ell}c_{k}^{(z_{j},z_{\ell})}\right)$
$\displaystyle=\frac{1}{n}\sum_{j\in Q\setminus
i}\sigma_{j}\left(\sum_{\begin{subarray}{c}\ell\neq j\\\
z_{\ell}=z_{j}\end{subarray}}c_{k}^{\mathrm{in}}-\sum_{\begin{subarray}{c}\ell\neq
j\\\ z_{\ell}\neq z_{j}\end{subarray}}c_{k}^{\mathrm{out}}\right)$
$\displaystyle\doteq\frac{1}{2}\sum_{j\in Q\setminus
i}\sigma_{j}\left(c_{k}^{\mathrm{in}}-c_{k}^{\mathrm{out}}\right)$
$\displaystyle=\frac{1}{2}\left(c_{k}^{\mathrm{in}}-c_{k}^{\mathrm{out}}\right)\tilde{v}_{iQ}\;.$

For the third line, we have used the fact that there are approximately
$\frac{n}{2}$ terms in each sum. This completes the proof.

##  Appendix D Proof of Theorem 5.1

######  Theorem D.1.

There exists a matrix $\mathbf{J}^{\prime}\in L(2\times\mathcal{Z}\times
K\times\mathcal{N})$ such that any eigenpair $(\lambda,\boldsymbol{\epsilon})$
of $\mathbf{J}$ satisfies exactly one of the following two alternatives:

  1. 1.

$(\lambda,\mathbf{L}\boldsymbol{\epsilon})$ is an eigenpair of
$\mathbf{J}^{\prime}$.

  2. 2.

$(\lambda,\mathbf{M}\boldsymbol{\epsilon})$ is an eigenpair of $\mathbf{N}$.

The operators $\mathbf{L}$, $\mathbf{M}$, and $\mathbf{N}$ are defined in
Section 5.2.

######  Proof D.2.

Let $(\lambda,\boldsymbol{\epsilon})$ be an eigenpair of $\mathbf{J}$. Then,
writing out the eigenvector-eigenvalue relation, we have

$\displaystyle\lambda\epsilon_{iR}^{(s)}$ $\displaystyle=\sum_{k\in
K}\sum_{t\in\mathcal{Z}}g_{k}^{(s,t)}\sum_{Q\in\mathcal{E}_{k}(i)\setminus
R}\sum_{\begin{subarray}{c}j\in Q\setminus
i\end{subarray}}\epsilon_{jQ}^{(t)}$ (49) $\displaystyle=\sum_{k\in
K}\sum_{t\in\mathcal{Z}}g_{k}^{(s,t)}\left(\beta_{ik}^{(s)}-\delta_{\lvert
R\rvert,k}\sum_{j\in R\setminus i}\epsilon_{jR}^{(t)}\right)\;,$

where we have inserted Eq. 17. Suppose first that
$\boldsymbol{\epsilon}\in\ker\mathbf{M}$. Then,
$\boldsymbol{\alpha}=\boldsymbol{\beta}=\mathbf{0}$. Simplifying, Eq. 49
becomes

(50) $\displaystyle\lambda\epsilon_{iR}^{(s)}$
$\displaystyle=-\sum_{t\in\mathcal{Z}}g_{\lvert
R\rvert}^{(s,t)}\sum_{\begin{subarray}{c}j\in R\setminus
i\end{subarray}}\epsilon_{jR}^{(t)}\;.$

Summing over $i$, we obtain

(51) $\displaystyle\lambda u_{R}^{(s)}=(1-\lvert
R\rvert)\sum_{t\in\mathcal{Z}}g_{\lvert R\rvert}^{(s,t)}u_{R}^{(t)}\;,$

where $\mathbf{u}=\mathbf{M}\boldsymbol{\epsilon}$. This states that
$\mathbf{u}$ is an eigenvector of $\mathbf{N}$ with eigenvalue $\lambda$,
which establishes Case (2) of Theorem 5.1.

Now suppose that $\boldsymbol{\epsilon}\notin\ker\mathbf{M}$. Proceeding from
Eq. 49 and summing over edges of size $k$ incident to $i$ yields

$\displaystyle\lambda\alpha_{ik}^{(s)}$
$\displaystyle=\sum_{R\in\mathcal{E}_{k}(i)}\sum_{k^{\prime}\in
K}\sum_{t\in\mathcal{Z}}g_{k^{\prime}}^{(s,t)}\left(\beta_{i,k^{\prime}}^{(t)}-\delta_{\lvert
R\rvert,k^{\prime}}\sum_{j\in R\setminus i}\epsilon_{jR}^{(t)}\right)$
$\displaystyle=\sum_{k^{\prime}\in
K}\sum_{t\in\mathcal{Z}}g_{k^{\prime}}^{(s,t)}\sum_{R\in\mathcal{E}_{k}(i)}\left(\beta_{ik^{\prime}}^{(t)}-\delta_{\lvert
R\rvert,k^{\prime}}\sum_{j\in R\setminus i}\epsilon_{jR}^{(t)}\right)$
$\displaystyle=\sum_{k^{\prime}\in
K}\sum_{t\in\mathcal{Z}}g_{k^{\prime}}^{(s,t)}\left(d_{i,k}\beta_{i,k^{\prime}}^{(t)}-\delta_{kk^{\prime}}\sum_{R\in\mathcal{E}_{k}(i)}\sum_{j\in
R\setminus i}\epsilon_{jR}^{(t)}\right)$ (52)
$\displaystyle=d_{ik}\sum_{k^{\prime}\in
K}\sum_{t\in\mathcal{Z}}g_{k^{\prime}}^{(s,t)}\beta_{ik^{\prime}}^{(t)}-\sum_{t\in\mathcal{Z}}g_{k}^{(s,t)}\beta_{ik}^{(t)}\;,$

which a closed linear system in the variables $\boldsymbol{\alpha}$ and
$\boldsymbol{\beta}$. To obtain a second system, we sum over nodes in size-$k$
edges incident to $i$, obtaining

$\displaystyle\lambda\beta_{ik}^{(s)}$
$\displaystyle=\sum_{R\in\mathcal{E}_{k}(i)}\sum_{j\in R\setminus
i}\sum_{k^{\prime}\in