# On the efficacy of higher-order spectral clustering under weighted
stochastic block models

Xiao Guo†, Hai Zhang†, Xiangyu Chang‡
† School of Mathematics, Northwest University, China
‡ School of Management, Xi’an Jiaotong University, China Xiangyu Chang is the
corresponding author.

###### Abstract

Higher-order structures of networks, namely, small subgraphs of networks (also
called network motifs), are widely known to be crucial and essential to the
organization of networks. Several works have studied the community detection
problem–a fundamental problem in network analysis at the level of motifs. In
particular, the higher-order spectral clustering has been developed, where the
notion of _motif adjacency matrix_ is introduced as the algorithm’s input.
However, how the higher-order spectral clustering works and when it performs
better than its edge-based counterpart remain largely unknown. To elucidate
these problems, the higher-order spectral clustering is investigated from a
statistical perspective. The clustering performance of the higher-order
spectral clustering is theoretically studied under a weighted stochastic block
model, and the resulting bounds are compared with the corresponding results of
the edge-based spectral clustering. The upper bounds and simulations show that
when the network is dense with the weak signal of weights, higher-order
spectral clustering can lead to a performance gain in clustering. Real data
experiments also corroborate the merits of higher-order spectral clustering.

Keywords: Higher-order structures, Community detection, Weighted networks,
Network motifs

##  1 Introduction

The network is a standard representation of relationships among units of
complex systems in many scientific domains, including biology, ecology,
sociology, and information science, among many others (Newman, 2018;
Goldenberg et al., 2010; Kolaczyk, 2009). One of the most common features of
networks is that they have communities, clusters, or modules–groups of nodes
that are, in some sense, more similar to nodes within the same community than
to other nodes. Community detection, namely, detecting such communities based
on given network structures, has been one of the fundamental problems in
network analysis which helps us gain knowledge about the behavior and
functionality of networks.

Past decades have seen various community detection methods, including spectral
clustering, modularity maximization, likelihood methods, and semidefinite
programming. See Abbe (2018) for a survey. Most of the current community
detection procedures are edge-based, that is, those procedures essentially
make use of the network adjacency matrix, which stores the similarity of each
pair of nodes to detect communities. However, there has been increasing
evidence that _higher-order_ connectivity patterns are crucial to the
organization of real networks. They can help understand and gain more insights
into the behavior of networks (Holland and Leinhardt, 1977; Milo et al., 2002;
Mangan et al., 2003; Yang and Leskovec, 2014; Rosvall et al., 2014). The
higher-order structures of networks often refer to small network subgraphs, or
the so-called network _motifs_ , in contrast to the _lower-order_ network
edges. In particular, triangular motifs are common in many kinds of networks.
For example, in social networks, individuals who share a common friend are
more likely to become friends, leading to a large proportion of triangles to
the two-edge wedges (Rohe and Qin, 2013).

Very recently, several works have studied the problem of network clustering at
the level of motifs, which can provide new insights into network organization
beyond the edge-based clustering if the network exhibits rich higher-order
structures; see Benson et al. (2016); Tsourakakis et al. (2017); Laenen and
Sun (2020); Underwood et al. (2020), among others. Generally, to capture the
higher-order connective patterns, a _motif adjacency matrix_ is often designed
as a surrogate for the original edge adjacency matrix; see Benson et al.
(2016); Paul et al. (2018) for example. To be specific, the $(i,j)$’s entry of
the motif adjacency matrix is the number of certain motifs that nodes $i$ and
$j$ participate in commonly in binary networks. It is suggested that one could
use the motif adjacency matrix in the subsequent clustering analysis to obtain
a good empirical result. In particular, the _higher-order spectral clustering_
, namely, the spectral clustering with the motif adjacency matrix as its
input, is well-suited to the problem and is the focus of this work. However,
it remains largely unknown how higher-order spectral clustering actually
works. Moreover, when it is better than its edge-based counterpart remains to
be seen. From a statistical perspective, the aim should be to show the
efficacy of higher-order spectral clustering on a particular network dataset
and to understand its merits under some underlying mechanisms.

In this paper, we attempt to illustrate the efficacy of higher-order spectral
clustering in the context of weighted networks. To this end, we develop the
_weighted stochastic block models_ (WSBMs) to mimic the weighted networks with
communities, which is a generalization of the stochastic block models (SBMs)
(Holland et al., 1983) to weighted networks. In particular, we assume the
edges are first generated according to the SBM mechanism, where nodes are
partitioned into several distinct communities and conditioned on the
underlying community assignments. The edges are generated independently
according to the community membership of their end nodes. After that, positive
weights independently generated from a common distribution are assigned to all
the existing edges. We specify the expectation of weights without requiring
their distributions. As a generalization of the aforementioned motif adjacency
matrix, we define the _weighted motif adjacency matrix_ to capture the network
weights, whose $(i,j)$’s entry is the weights summation of a certain motif
that nodes $i$ and $j$ both involve. In particular, we focus on the undirected
triangular motif. With these at hand, we study the approximation error and
misclustering error of the higher-order spectral clustering under the WSBM
mechanism and compare the upper bounds with the corresponding results of edge-
based spectral clustering (i.e., spectral clustering with the weighted edge
adjacency matrix as the input). It turns out that when the network is dense
_but_ the conditional expectation of edge weight is small, and the higher-
order spectral clustering leads to a lower misclustering error than the edge-
based one does. The rationality is as follows. The signal of communities could
be weak under small edge weight, which may lead to unsatisfactory clustering
performance of the edge-based method. However, when the network is dense, the
triangles are rich, which enables each node pair to borrow the strength and
information of incident edges. Thus, the clustering quality could be improved
by utilizing the triangular motif. Note that the unequal weights are crucial
to ensure the advantage of the higher-order spectral clustering over the edge-
based one. In addition, it is worth noting that, theoretically, the weighted
motif adjacency matrix leads to a larger eigen-gap of the corresponding
population matrix than the edge adjacency matrix does. In this sense, the
signals of the weighted networks are thus enhanced. We conduct simulations to
verify the theoretical findings. In particular, the higher-order spectral
clustering has performance gain over the edge-based counterpart in various
regimes where the requirements for theory are not necessarily met. The real
data experiments on two statisticians networks and one ${\rm PM}_{2.5}$
network provide insightful and interesting results, showing the efficacy of
higher-order spectral clustering on weighted networks.

The reminder of this paper is organized as follows. Subsection 1.1 reviews the
related works. Section 2 provides the general framework including the
notation, weighted stochastic block models and higher-order spectral
clustering. Section 3 presents the theoretical results. Section 4 discusses
the theoretical and methodological aspects of our work and also poses possible
extensions. Sections 5 and 6 include the simulations and real data analysis,
respectively. Section 7 concludes the paper. The proofs can be found in the
Supplementary Material.

###  1.1 Related works

Analyzing networks at the level of motifs has received increasing attention in
recent years. In particular, Benson et al. (2016) proposed a conductance-based
method, which generalized the original edge-based conductance to the motif-
based one. The motif adjacency matrix was then introduced to simplify the
minimization of the motif-based conductance; more specifically, minimizing the
motif-based conductance can be reduced to minimizing the edge-based
conductance but using the motif adjacency matrix. Serrour et al. (2011)
studied the higher-order network clustering in the modularity maximization
framework, where they generalized the notion of modularity by transforming its
building block from edges to motifs. They employed the tensor to capture the
triangular patterns within the network but then transformed it to the
aforementioned motif adjacency matrix to accelerate the algorithm. There are
also works focusing on directed networks (Laenen and Sun, 2020; Cucuringu et
al., 2020; Underwood et al., 2020) and local spectral clustering (Yin et al.,
2017).

The statistical aspects of the motif-based methods are also studied by several
authors. Rohe and Qin (2013) show the blessings of transitivity in sparse
networks. They developed an algorithm that exploits the triangles built by
network transitivity and shows its clustering ability statistically under the
newly developed local stochastic block models. Tsourakakis et al. (2017) also
studied the motif-based conductance and provided a random walk explanation of
it. They showed that the random walk on the network corresponding to the motif
adjacency matrix is more likely to stay in the same true community of the
stochastic block model (Holland et al., 1983) than the random walk on the
original network. Paul et al. (2018) proposed a superimposed stochastic block
model, which is a superimposition of a classical dyadic (edge-based) random
graph and a triadic (triangle-based) random graph. They rigorously analyzed
the misclustering error bound of the higher-order spectral clustering under
such models. However, the analysis fails to disclose the advantage of higher-
order spectral clustering over the edge-based counterpart. In addition, the
model is somewhat restrictive since the networks generated by the superimposed
stochastic block models have edges lying in $\\{0,1,2\\}$, which infrequently
happens in real networks. Cucuringu et al. (2020) studies the spectral
clustering on complex-valued Hermitian matrix representations, which
implicitly use the higher-order structure of directed networks. They studied
the algorithm theoretically under stochastic block models.

There exists another line of works focusing on the hyper-graph clustering
problem; see Ghoshdastidar et al. (2017); Ghoshdastidar and Dukkipati (2017,
2014), among many others. In hyper-graphs, the hyper-edges are directly known
as a prior rather than constructed using the motifs. Theoretically, the motif
adjacency matrix brings extra dependence between entries. Nevertheless, in the
context of hyper-graph clustering, the entries of the adjacency tensor are
independent.

There are also some authors making use of the motifs to study the testing
problem of whether it has only one community or multiple communities under
SBMs or its variants; notably by Jin et al. (2018); Gao et al. (2017); Jin et
al. (2021); Cammarata and Ke (2022).

On the other hand, there exist a few works on community detection for weighted
networks; see Aicher et al. (2015); Xu et al. (2020); Cerqueira and Levina
(2023) for example and references therein. Compared with these works, we are
motivated to enhance the signal strength in weighted networks by using the
weighted motif adjacency matrix and studying its effect theoretically.
Generally, using the weighted motif adjacency matrix can be regarded as a
network preprocessing step, and the weighted motif adjacency matrix can
essentially be the input of other community detection algorithms.

##  2 Framework

In this section, we introduce the general framework of analysis. In
particular, the WSBMs and the higher-order spectral clustering are introduced.

We first introduce some notes and notation. The readers could also refer to
Table 1 for a brief summary. $\mathbb{M}_{n,K}$ denotes the set of all
$n\times K$ matrices which have exactly one 1 and $K-1$ 0’s in each row. Any
$\Theta\in\mathbb{M}_{n,K}$ is called a _membership matrix_ , where each row
corresponds to the community membership of the corresponding node. For
example, node $i$ belongs to community $g_{i}\in\\{1,...,K\\}$ if and only if
$\Theta_{ig_{i}}=1$. For $1\leq k\leq K$, denote $G_{k}=G_{k}(\Theta)=\\{1\leq
i\leq n:g_{i}=k\\}$ and $n_{k}=|G_{k}|$, where $G_{k}$ consists of nodes with
their community membership being $k$. For any matrix $A_{n\times n}$ and
$I,J\subseteq\\{1,...,n\\}$, $A_{I\ast}$ and $A_{\ast J}$ denote the submatrix
of $A$ consisting of the corresponding rows and columns of $A$, respectively.
$\|A\|_{\tiny{\rm F}}$, $\|A\|_{2}$, and $\|A\|_{\infty}$ denote the Frobenius
norm, spectral norm, and the element-wise maximum absolute value of $A$,
respectively. ${\rm diag}(A)$ denotes the diagonal matrix with its diagonal
elements being the same as those of $A$. We also use $\|a\|_{2}$ to denote the
Euclidean norm of any vector $a$.

Now we introduce the WSBMs which generalizes the SBMs (Holland et al., 1983)
such that the potential network can have weighted edges rather than the binary
edges in the SBMs. For an underlying network with $n$ nodes and $K$
communities, the two main parameter matrices of the WSBMs are the membership
matrix $\Theta\in\mathbb{M}_{n,K}$, and the connectivity matrix
$B\in[0,1]^{K\times K}$ where $B$ is of full rank, symmetric, and the entry
$B_{kl}$ of $B$ represents the edge probability between any node in community
$l$ and any node in community $k$. For simplicity, throughout the paper we
assume the underlying network has balanced community size $n/K$ with the
within-community probability being $p_{n}$ and the between-community
probability being $p_{n}(1-\lambda)$ for any pair of nodes, where $\lambda$ is
a constant in $(0,1)$. Whereas note that the theoretical results, say the
approximation error bound to be established in Theorem 1, also holds for
general $B$’s; see Remark 3. Under the simplified structural assumption, the
connectivity matrix $B$ takes the simple form

$B=p_{n}\lambda I_{K}+p_{n}(1-\lambda)1_{K}1_{K}^{\intercal},$ (2.1)

where $1_{K}$ represents a $K$-dimensional vector of 1’s. Further, we
introduce a probability distribution $f$ supported on the positive line which
is actually used to generate the edge weights. We do not specify the
distribution of $f$. We assume that the expectation of its corresponding
random variable is $\alpha_{n}$, and the variance is larger than 0. Although
we assume an community-independent weight here, we note that the results can
be extended to community-independent case; see Remarks 1 and 3.

Given $\Theta$, $B$, and $f$, the weighted network adjacency matrix
$W=(W_{ij})_{1\leq i,j\leq n}$ is generated as

$W_{ij}=\begin{cases}{\rm Bernoulli}(B_{g_{i}g_{j}})\cdot V_{ij}&\mbox{if
}\;i<j,\\\ 0,&\mbox{if }\;i=j,\\\ W_{ji},&\mbox{if }\;i>j.\end{cases}$

where ${\rm Bernoulli}(p)$ denotes the bernoulli random variable which is 1
with probability $p$ and is 0 with probability $1-p$, $V_{ij}$’s denote the
random variables generated from distribution $f$, and $V_{ij}$’s are mutually
independent and they are also independent of $B_{g_{i}g_{j}}$’s. Define
$\mathcal{W}=\alpha_{n}\Theta B\Theta^{\intercal}$, and it is then easy to see
that $\mathcal{W}$ is the population of $W$ in the sense that
$\mathbb{E}(W)=\mathcal{W}-{\rm diag}(\mathcal{W})$. In order to capture the
higher-order connectivity patterns of the weighted network, we define the
following weighted motif (triangular) adjacency matrix $W^{M}=(W^{M}_{ij})$ to
be

$W^{M}_{ij}=\sum_{k=1}^{n}\mathbf{1}(W_{ij}\cdot W_{ik}\cdot
W_{jk}>0)(W_{ij}+W_{ik}+W_{jk}),$ (2.3)

where $\mathbf{1}(\cdot)$ stands for the indicator function. Recall that in
Benson et al. (2016); Tsourakakis et al. (2017), the motif (triangular)
adjacency matrix was defined to be a matrix with its $(i,j)$’s entry being the
number of triangles that nodes $i$ and $j$ participate in commonly. Thus (2.3)
is a simple generalization of the motif adjacency matrix proposed by Benson et
al. (2016); Tsourakakis et al. (2017) in order to handle the weighted
networks. Now, to facilitate further analysis, let us have a closer look at
the expectation of $W^{M}$. Recall (2.1) and (2), and then we can easily
obtain the following observations. When $g_{i}=g_{j}$,

$\displaystyle\mathbb{E}(W_{ij}^{M})$ $\displaystyle=3\alpha_{n}\cdot
p_{n}\\{(\frac{n}{K}-2)p_{n}^{2}+(K-1)\frac{n}{K}p_{n}^{2}(1-\lambda)^{2}\\}:={\equiv
h_{1}},$ (2.4)

and when $g_{i}\neq g_{j}$,

$\displaystyle\mathbb{E}(W_{ij}^{M})=3\alpha_{n}\cdot p_{n}(1-\lambda)$
$\displaystyle\\{2(\frac{n}{K}-1)p_{n}^{2}(1-\lambda)+(K-2)\frac{n}{K}p_{n}^{2}(1-\lambda)^{2}\\}{:=h_{2}}.$
(2.5)

Define

$\mathcal{W}^{M}=\Theta((h_{1}-h_{2})I_{K}+h_{2}1_{K}1_{K}^{\intercal})\Theta^{\intercal},$
(2.6)

then $\mathcal{W}^{M}$ is the population of $W^{M}$ by noting

$\mathbb{E}(W^{M})=\mathcal{W}^{M}-\mbox{diag}(\mathcal{W}^{M}).$

###### Remark 1

(2.4) and (2.5) can be similarly computed under more general $B$’s and
community-dependent weights. For example, assume that the within- and between-
community probabilities are $a_{n}$ and $b_{n}$, respectively, and the within-
and between-community weights are $\alpha_{n}$ and $\beta_{n}$, respectively,
then within- and between-community expectations of $W^{M}$ are

$h_{1}:=a_{n}(\frac{n}{K}-2)a_{n}^{2}\cdot
3\alpha_{n}+a_{n}(K-1)\frac{n}{K}b_{n}^{2}(\alpha_{n}+2\beta_{n}),$
$h_{2}:=2b_{n}(\frac{n}{K}-1)a_{n}b_{n}(\alpha_{n}+2\beta_{n})+b_{n}(K-2)\frac{n}{K}b_{n}^{2}\cdot(3\beta_{n}),$

which are crucial for bounding the eigen-gap and thus the misclustering rate;
see (3.4).

With these formulations at hand, we now introduce the algorithm for community
detection where the goal is to recover the membership matrix $\Theta$ up to
some column permutations. In this paper, we study the spectral clustering (Von
Luxburg, 2007) which generally consists of two steps. The first step is to
perform the eigenvalue decomposition of a suitable matrix representing the
network. The next step is to run $k$-means on the resulting $K$ leading
eigenvectors. We consider two kinds of spectral clustering. We will mainly
deal with the so-called higher-order spectral clustering or the motif-based
spectral clustering; see Algorithm 1 for details. And we would compare it with
the edge-based spectral clustering, namely, Algorithm 1 without _step 1_ and
with _step 2_ replaced by the eigenvalue decomposition of the weighted
adjacency matrix $W$. The goal is to study how the higher-order spectral
clustering performs under the WSBMs and to understand how it can enhance the
clustering performance compared with that of the edge-based spectral
clustering.

###### Remark 2

Similar to Algorithm 1, Benson et al. (2016) and Serrour et al. (2011) also
proposed spectral-based algorithms by utilizing higher-order motifs. The
sweep-cut algorithm in Benson et al. (2016) essentially minimizes the motif-
based conductance, which is the normalized number of motifs that are cut of
two detected communities. The spectral-based algorithm in Serrour et al.
(2011) approximately maximizes the motif-based modularity, which further
reduces to the modularity with respect to motif adjacency matrix. In this
sense, we expect that Algorithm 1 would lead to larger modularity (w.r.t.
weighted motif adjacency matrix) than the edge-based counterpart does. This
actually will be verified by simulated experiments in Section 5.

The next lemma shows that $\mathcal{W}^{M}$, namely, the population version of
$W^{M}$, has eigenvectors that reveal the true communities.

###### Lemma 1

For a WSBM with $K$ communities and $B$ structured as in (2.1), suppose
$\mathcal{W}^{M}=U\Sigma U^{\intercal}$ is the eigenvalue decomposition of
$\mathcal{W}^{M}$ (see 2.6), then $U=\Theta X$. Specifically, for
$\Theta_{i\ast}=\Theta_{j\ast}$, we have ${U}_{i\ast}={U}_{j\ast}$; while for
$\Theta_{i\ast}\neq\Theta_{j\ast}$, we have
$\|{U}_{i\ast}-{U}_{j\ast}\|_{2}=\sqrt{(n_{g_{i}})^{-1}+(n_{g_{j}})^{-1}}$.

It is easy to see from Lemma 1 that the population eigenvector has $K$
distinct rows ($n$ rows in total) and two rows are identical if and only if
the corresponding nodes are in the same underlying community. Therefore, the
higher-order spectral clustering would cluster well if the sample version
eigenvectors are concentrated around its expectation. We will discuss its
theoretical properties in the next section.

Algorithm 1 Higher-order spectral clustering for $K$ communities

0: Community number $K$, a weighted adjacency matrix $W\in\mathbb{R}^{n\times
n}$ representing the network;

0: Community membership;

1: Motif matrix construction: Construct the weighted motif adjacency matrix
$W^{M}$ according to (2.3).

2: Eigen-decomposition: Find the $K$ leading eigenvectors of $W^{M}$ and form
the $n\times K$ eigen-matrix.

3: $k$-means: Treat each row of the eigen-matrix as a point in
$\mathbb{R}^{K}$ and run the Lloyd’s algorithm on these points to solve
$k$-means.

Table 1: A summary of the notes and notations  Notation | Definition
---|---
$n$ | Number of nodes
$K$ | Number of communities
$\Theta$ | Membership matrix
$B$ | Connectivity matrix
$\mathcal{W}$ | Population edge-based weighted adjacency matrix
$\mathcal{W}^{M}$ | Population motif-based weighted adjacency matrix
$W$ | Edge-based weighted adjacency matrix
$W^{M}$ | Motif-based weighted adjacency matrix
$g_{i}$ |  The community which node $i$ belongs to
$G_{k}$ |  Set of nodes from community $k$
$n_{k}$ |  The cardinality of $G_{k}$
$p_{n}$ |  The maximum link probability in $\mathcal{W}$
$\alpha_{n}$ | The expectation of edge weights when edge exists

##  3 Theoretical analysis

In this section, we theoretically justify the clustering performance of the
higher-order spectral clustering under the model set-up of the WSBMs and then
we compare the theoretical bounds with those of the edge-based spectral
clustering.

###  3.1 Higher-order spectral clustering

Recall that $W$ is the weighted adjacency matrix of a weighted graph generated
from the WSBMs with $n$ nodes and $K$ communities (see (2)), and $\mathcal{W}$
denotes the population of $W$. $p_{n}$ is the maximum linking probability
between two nodes in the WSBMs (see (2.1)) and $\alpha_{n}$ is the conditional
expectation of an edge given that the edge exists. For simplicity, we assume
that $W_{ij}\leq 1$ for any $1\leq i,j\leq n$. This is reasonable since we
will see in (3.2) that we require a vanishing $\alpha_{n}$. On the other hand,
one could specify the distribution of weights and using concentration
inequalities to derive a ${\rm log}n$-type bound. $W^{M}$ is the corresponding
weighted motif adjacency matrix (see (2.3)), where each entry is defined to be
the weighted sum of the triangles that the corresponding nodes both join, and
$\mathcal{W}^{M}$ denotes the population of $W^{M}$ except the diagonal
elements. The next theorem provides the concentration bound of the weighted
motif adjacency matrix around its expectation.

###### Theorem 1

Assume $p_{n}\geq c^{\prime}\sqrt{\frac{{\rm log}n}{n}}$ for some constant
$c^{\prime}>0$, and $\alpha_{n}p_{n}\geq c^{\prime\prime}\frac{({\rm
log}n)^{\beta}}{n}$ for some constants $c^{\prime\prime}>0$ and any
$0<\beta<1$. Let ${\tau_{\rm max}=np_{n}^{2}}$ and

${D=\tau_{\rm max}^{2}n\alpha_{n}p_{n}=n^{3}\alpha_{n}p_{n}^{5}}.$

Then for some constant $r>0$, there exists a constant $c$ such that with
probability at least $1-n^{-r}-{\rm exp}(-r({\rm log}n)^{\beta})$,

${\|W^{M}-\mathbb{E}(W^{M})\|_{2}\leq c\sqrt{D}},$ (3.1)

where $c$ depends on $r$ and $\beta$.

###### Remark 3

Note that we make use of neither the structure of the connectivity matrix $B$
nor the weights distribution to derive the results. Indeed, (3.1) holds for
more general WSBMs, as long as the weights’ expectations and connectivity
matrix’s entries are upper bounded.

Theorem 1 requires that the network is not too sparse in the sense that
$p_{n}\geq c^{\prime}\sqrt{\frac{{\rm log}n}{n}}$ for some constant
$c^{\prime}>0$, where actually ensures that the network motifs are not rare.
We also require the overall signal strength is lower bounded, namely,
$\alpha_{n}p_{n}\geq c^{\prime\prime}\frac{({\rm log}n)^{\beta}}{n}$ for some
constants $c^{\prime\prime}>0$ and $0<\beta<1$. Theorem 1 indicates that the
weighted motif adjacency matrix concentrates around its population version at
the rate of $\sqrt{n^{3}\alpha_{n}p_{n}^{5}}$. The resulting bound is similar
to the Theorem 2 of Paul et al. (2018) except that $\alpha_{n}$ arises in our
result. Recall that for the edge adjacency matrix, its spectral bound is the
square root of the maximum expected degree (Lei and Rinaldo, 2015). Applying
this law to the weighted motif adjacency matrix, its spectral bound would read
as the square root of the maximum expected “motif degree”, which is actually
$\sqrt{n^{2}\alpha_{n}p_{n}^{3}}$ as can be seen in the proof. Comparing with
$\sqrt{n^{2}\alpha_{n}p_{n}^{3}}$, the bound $\sqrt{n^{3}\alpha_{n}p_{n}^{5}}$
in (3.1) is loose in this sense if we note that $p_{n}\geq
c^{\prime}\sqrt{\frac{{\rm log}n}{n}}$ for some constant $c^{\prime}>0$.

Now we are ready to study the clustering performance of the higher-order
spectral clustering. Specifically, we use the following metric to evaluate the
quality of clustering,

$\displaystyle L(\hat{{\Theta}},\Theta)=\underset{J\in E_{K}}{\rm
min}\,\underset{1\leq k\leq
K}{\sum}\;(2n_{k})^{-1}\|(\hat{{\Theta}}J)_{G_{k\ast}}-\Theta_{G_{k\ast}}\|_{0},$
(3.2)

where $\hat{{\Theta}}$ denotes the estimated membership matrix by the higher-
order spectral clustering, and $E_{K}$ denotes the set of all $K\times K$
permutation matrices. Obviously, $L$ measures the sum of the fractions of the
misclustered nodes within each community. The following theorem provides an
upper bound on $L$.

###### Theorem 2

Let $\hat{{\Theta}}$ be the estimated membership matrix by the higher-order
spectral clustering, and denote the minimum non-zero eigenvalue of
$\mathcal{W}^{M}$ as $\lambda_{K}(\mathcal{W}^{M})$. Suppose the assumptions
in Theorem 1 hold, and there exists a constant $c>0$ such that, if

$\frac{K^{3}}{n\alpha_{n}p_{n}}\leq 1/c,$ (3.3)

then with probability at least $1-n^{-r}-{\rm exp}(-r({\rm log}n)^{\beta})$
for some constant $r>0$ and any $0<\beta<1$, there exist subsets
$S_{k}\subseteq G_{k}$ for $k=1,...,K$ such that

$L(\hat{{\Theta}},\Theta)\leq\sum_{k=1}^{K}\frac{|S_{k}|}{n_{k}}\leq
c\frac{K\|W^{M}-\mathcal{W}^{M}\|_{2}^{2}}{\lambda_{K}^{2}(\mathcal{W}^{M})}\leq
c\frac{K^{3}}{n\alpha_{n}p_{n}}.$ (3.4)

Moreover, for $G=\cup_{k=1}^{K}(G_{k}\backslash S_{k})$, there exists a
$K\times K$ permutation matrix $J$ such that

$\hat{\Theta}_{G_{\ast}}J=\Theta_{G_{\ast}},$ (3.5)

where we recall that $\hat{\Theta}_{G_{\ast}}$ and ${\Theta}_{G_{\ast}}$
denote the submatrix of $\hat{\Theta}$ and $\Theta$ consisting of the rows
indexed by $G$.

Theorem 2 indicates that the misclustering rate of the motif-based spectral
clustering is bounded by $O(\frac{K^{3}}{n\alpha_{n}p_{n}})$ with high
probability. (3.3) is a technical condition that provides the range of
parameters $(K,n,\alpha_{n},p_{n})$ under which the conclusions hold. $S_{k}$
in (3.4) is actually the set of misclustered nodes in the true community
$G_{k}$. The result in (3.4) is high-dimensional in that each parameter can
vary with the number of nodes $n$. The bound vanishes under several parameter
settings. For example, when $\alpha_{n}\asymp{\rm log}n/{n}$ and $p_{n}\asymp
1/({\rm log}n)^{{1-\beta}}$ where $\beta\in(0,1)$ and $f(n)\asymp g(n)$ if
$cg(n)\leq f(n)\leq Cg(n)$ for some constants $0<c<C<\infty$, then $K=o(({\rm
log}n)^{\tiny\beta/3})$ ensures a vanishing misclustering bound.

###  3.2 Comparison with edge-based spectral clustering

To understand how the higher-order spectral clustering enhances the
misclustering performance, we compare its misclustering bounds with those from
edge-based spectral clustering. Note that edge-based spectral clustering has
been studied by Lei and Rinaldo (2015); Rohe et al. (2011) under the
_unweighted_ SBMs. One can easily modify the proofs in Lei and Rinaldo (2015)
to obtain the corresponding results under the WSBMs, hence we omit the
details. It is easy to learn from Theorem 2 that the spectral bound
($\|A-\mathcal{A}\|_{2}$) and the minimum non-zero eigenvalue of the
population input matrix ($\lambda_{K}(\mathcal{A})$) have the following
relationship with the misclustering rate $L$ defined in (3.2),

$L\leq c\frac{K\|A-\mathcal{A}\|_{2}^{2}}{\lambda_{K}(\mathcal{A})^{2}},$
(3.6)

where $c$ is some constant, $A$ can be $W$ or $W^{M}$, and $\mathcal{A}$ can
be $\mathcal{W}$ or $\mathcal{W}^{M}$ correspondingly. Hence, we list these
three metrics of the two spectral clustering methods in Table 2. The bounds in
Table 2 show that when

$\frac{cK^{3}}{n\alpha_{n}p_{n}}\leq\frac{c^{\prime\prime}K^{3}\mbox{log}n}{n^{2}\alpha_{n}^{2}p_{n}^{2}},$
(3.7)

the higher-order spectral clustering has a lower misclustering rate than does
the edge-based counterpart. Combining (3.7) with the parameter assumption in
Theorem 1 leads to the following parameter ranges

$\begin{cases}\alpha_{n}\leq c\frac{{\rm log}n}{n}/p_{n},\\\ p_{n}\geq
c^{\prime}\sqrt{\frac{{\rm log}n}{n}},\\\ \alpha_{n}p_{n}\geq
c^{\prime\prime}\frac{({\rm log}n)^{\beta}}{n},\end{cases}$

which can be all met when, for example, $\alpha_{n}\asymp{\rm log}n/{n}$ and
$p_{n}\asymp 1/({\rm log}n)^{{1-\beta}}$, where $\beta\in(0,1)$ is any
constant that arises in Theorem 1. The first inequality in (3.2) requires that
the signal strength is weak if we notice the second inequality. The second
inequality says that the network density should be large at a certain level.
And the third inequality could be thought of as a requirement of the overall
signal strength of weight and network density. As a result, when the network
is dense with a weak signal of weights, the higher-order spectral clustering
can be better than its edge-based counterpart in terms of the clustering
error. In addition, although we only specify the expectation of conditional
weights, we implicitly assume that the variance of weight is larger than 0. In
fact, the unweighted SBMs can not lead to the advantage of the higher-order
spectral clustering over the edge-based counterpart, which is consistent with
the results in Paul et al. (2018) (see Section 3.5 in Paul et al. (2018)).

###### Remark 4

It is worth noting that there exists statistical minimax misclustering error
rate under SBMs and a general weighted SBMs (Gao et al., 2017; Xu et al.,
2020). However, we here focus on enhancing the clustering performance of the
spectral clustering by using the weighted motif adjacency matrix. We show that
the higher-order spectral clustering leads to a lower misclustering upper
bound if the network is dense with a weak expected signal of weights, which is
partially because the motif-based adjacency matrix enlarges the eigen-gap
between the smallest non-zero eigenvalue and 0 of the population matrix (see
Table 2). Similar observations are obtained in Wang et al. (2018), and our
empirical results in Section 5 also validate the theoretical results. It is
left as future work to compare these two spectral clustering methods more
delicately, for example, deriving tighter bounds or studying the strong
consistency of both methods.

Table 2: A summary of the main results of two spectral clustering methods under the WSBMs.  Bounds | Edge-based | Motif-based
---|---|---
Spectral bound | $O_{p}(\sqrt{\mbox{max}\\{n\alpha_{n}p_{n},c\mbox{log}n\\}})$ | $O_{p}(\sqrt{n^{3}\alpha_{n}p_{n}^{5}})$
$\lambda_{K}$ | $O(\frac{n\alpha_{n}p_{n}}{K})$ | $O(\frac{n^{2}\alpha_{n}p_{n}^{3}}{K})$
misclustering rate | $O_{p}(\mbox{max}\\{\frac{K^{3}}{n\alpha_{n}p_{n}},\frac{K^{3}{\rm log}n}{n^{2}\alpha_{n}^{2}p_{n}^{2}}\\})$ | $O_{p}(\frac{K^{3}}{n\alpha_{n}p_{n}})$

##  4 Discussions

In this section, we discuss our work from methodology to theory. Possible
extensions are also posed.

#### Methodology.

Conceptually, although we focused on the spectral clustering, the notion of
motif adjacency matrix could be used as the input of many community detection
methods. From this point of view, it can be regarded as a data preprocessing
technique. Qin and Rohe (2013) showed that the regularized graph Laplacian in
which the degree matrix is regularized with a small constant can lead to
better clustering results than the original spectral clustering. Similarly,
using the motif adjacency matrix instead of the adjacency matrix can be
thought of as one network denoising technique.

On the other hand, driven by the evidence that network edges often contain
sensitive information, there is a growing body of works on the privacy-
preserving analysis of networks; see Karwa and Slavković (2016); Karwa et al.
(2017), among others. From this perspective, the entries of the motif
adjacency matrix can be regarded as the summary statistics, and thus it can
help protect the original edges to some extent.

Computationally, there is no doubt that counting the number of triangles for
each pair of nodes is costly on large-scale networks. One could use sampling
techniques to improve the computational performance (Seshadhri et al., 2014;
Benson et al., 2016; Chen and Chen, 2018; Guo et al., 2020; Zhang et al.,
2022).

#### Theory.

Theorem 1 and Theorem 2 are our main results. Theorem 1 investigates the
deviation of the weighted motif adjacency matrix from its population version
in the sense of the spectral norm. We generally use the $\epsilon$-net
technique to make a discretization, and then use a combinatorial method to
bound the approximation error. This technique was first developed by Feige and
Ofek (2005). Since that time it has been widely used by statistics and machine
learning communities; see Lei and Rinaldo (2015); Paul et al. (2018); Gao et
al. (2017); Chin et al. (2015), among many others. Note that the entries of
$W^{M}$ are dependent which creates difficulty in deriving the bound. To
tackle this issue, we use similar arguments developed in Paul et al. (2018).
In particular, we make use of the _typical bounded difference inequality_
established in Warnke (2016) to handle the dependency. Note that as the
network edges in our model are not Bernoulli distributed and the resulting
weighted motif adjacency matrix is not identical to theirs, the details of our
proof are slightly different from Paul et al. (2018). It remains unclear
whether the bound in (3.1) could be improved by using other techniques.
Theorem 2 bounds the misclustering rate of the motif-based spectral
clustering. The general idea is to use the Davis-Kahan theorem (Davis and
Kahan, 1970) to bound the perturbation of eigenvectors from the approximation
error bound (3.1). Such framework is widely used in Lei and Rinaldo (2015);
Rohe et al. (2011), among others.

#### Extensions.

We mainly study the full rank WSBMs, that is, the rank of $\mathcal{W}$ equals
the target community number $K$. Actually, all the results could be
generalized to rank-deficient WSBMs (i.e., ${\rm rank}(\mathcal{W})\leq K$) by
investigating the population eigen-structure of such WSBMs and adding extra
conditions on $B$, just as the argument in Zhang et al. (2022). In addition,
though we mainly deal with the clustering performance, we can estimate the
connectivity matrix $B=(B_{ql})_{q,l=1,...,K}$ via the following simple plug-
in estimator,

$\hat{B}_{ql}:=\frac{\sum_{1\leq i,j\leq
n}W_{ij}^{M}\hat{\Theta}_{iq}\hat{\Theta}_{jl}}{\sum_{1\leq i,j\leq
n}\hat{\Theta}_{iq}\hat{\Theta}_{jl}},\quad 1\leq q,l\leq K.$

Moreover, we could evaluate its theoretical performance by using Theorem 1 and
2. Following the advantage of motif-based spectral clustering in terms of
clustering, we can imagine that the motif-based method would lead to better
estimates of $B$ than the edge-based method does. For the sake of space, we
leave all these extensions as our future work.

##  5 Simulation studies

In this section, we empirically compare the finite sample performance of the
higher-order spectral clustering, namely, the motif-based spectral clustering
(denoted by _motif-based_) with that of the edge-based spectral clustering
(denoted by _edge-based_) under the WSBMs. Note that we use higher-order
spectral clustering and motif-based spectral clustering interchangeably in
this section. To that end, we first carry out a sensitivity analysis to
evaluate our theoretical findings. Then, we provide several extended
experiments to show that high-order spectral clustering has an advantage under
a broad regime even when some theoretical assumptions are violated.

First, we evaluate the performance of two methods and validate the theoretical
results. To be consistent with Section 3, we use the following three metrics
to assess the performance of each method. The first is the spectral deviation
of the weighted motif adjacency matrix $W^{M}$ from its population version
$\mathcal{W}^{M}$, denoted by _spectral bound_. The second is the minimum non-
zero eigenvalue of the population $\mathcal{W}^{M}$, denoted by _eigen gap_.
The third is the summation of the ratio of misclustered nodes within each true
community (see (3.2)), denoted by _misclustering rate_. We study the effect of
sample size $n$, the effect of maximum link probability $p$, the effect of
out-in-ratio (the ratio of the between-community probability $q$ over the
within-community probability $p$), and the effect of the number of communities
$K$ respectively with the following experimental set-ups,

  * •

_Effect of $n$:_ $K=2,p=0.5,q=p(1-0.4)$, and $n$ varies;

  * •

_Effect of $p$:_ $K=2,n=60,q=p(1-0.4)$, and $p$ varies;

  * •

_Effect of out-in-ratio:_ $K=2,n=60,p=0.5$, and $\frac{q}{p}$ varies;

  * •

_Effect of $K$:_ $n=120,p=0.5,q=p(1-0.4)$, and $K$ varies;

where the weights are all _i.i.d._ generated from uniform distribution ${\rm
Uniform}(0.01,1)$ provided that there is an edge.

Figures 1-4 show the results for the above four experimental set-ups. As
indicated in Section 3, the spectral bound and the eigen-gap with respect to
the two spectral clustering methods have different scales. Hence we use two
axes to show their tendency. Specifically, the left and right axes correspond
to motif-based and edge-based spectral clustering, respectively. The numerical
results show satisfactory consistency with the theoretical results. For the
spectral bound, the motif-based spectral clustering grows super-linearly with
$n$ and $p$, while the edge-based spectral clustering grows sub-linearly with
$n$ and $p$ (see Figures 1(a) and 2(a)). Note that we have not taken the
community structure into consideration when bounding the spectral error. While
as indicated in Figures 3(a) and 4(a), the spectral bounds grow with the
increasing of out-in-ratio and drop with the increasing number of communities.
Hence, it would be beneficial to incorporate this information when bounding
the spectral error. We leave this for our future work. For the eigen-gap, as
$n$ and $p$ increase, the motif-based spectral clustering grows faster than
the edge-based spectral clustering does, where the latter has linear growth
with $n$ and $p$ (see Figures 1(b) and 2(b)). In addition, the eigen-gap
decreases with the growth of $K$ for both methods. As for the misclustering
rate, we can see from Figures 1(c), 3(c), and 4(c) that when $n$, out-in-
ratio, and $K$ are intermediate, the motif-based spectral clustering has great
advantage over the edge-based method in clustering. When these terms are too
small or too large, the signal for the communities is weak such that neither
method can recover the communities, or the signal is strong such that the
edge-based clustering can do well. In particular, we see from Figure 2(c) that
when the network is sparse, the edge-based spectral clustering is better than
the motif-based method. As the maximum linking probability $p$ increases, the
motif-based spectral clustering performs gradually better than the edge-based
method does. This is consistent with our theoretical findings.

(a) Spectral bound (b) Eigen gap (c) Misclustering rate

Figure 1: The average effect of $n$ on the three metrics over 100
replications. Other parameters $K=2,p=0.5,q=p(1-0.4)$, and the edge weights
are _i.i.d._ generated from ${\rm Uniform}(0.01,1)$ provided that the edges
exist. In (a) and (b), the left and right axes correspond to motif-based and
edge-based spectral clustering, respectively.

(a) Spectral bound (b) Eigen gap (c) Misclustering rate

Figure 2: The average effect of $p$ on the three metrics over 100
replications. Other parameters $K=2,n=60,q=p(1-0.4)$, and the edge weights are
_i.i.d._ generated from ${\rm Uniform}(0.01,1)$ provided that the edges exist.
In (a) and (b), the left and right axes correspond to motif-based and edge-
based spectral clustering, respectively.

(a) Spectral bound (b) Eigen gap (c) Misclustering rate

Figure 3: The average effect of out-in-ratio on the three metrics over 100
replications. Other parameters $K=2,n=60,p=0.5$, and the edge weights are
_i.i.d._ generated from ${\rm Uniform}(0.01,1)$ provided that the edges exist.
In (a) and (b), the left and right axes correspond to motif-based and edge-
based spectral clustering, respectively.

(a) Spectral bound (b) Eigen gap (c) Misclustering rate

Figure 4: The average effect of $K$ on the three metrics over 100
replications. Other parameters $n=120,p=0.5,q=p(1-0.4)$, and the edge weights
are _i.i.d._ generated from ${\rm Uniform}(0.01,1)$ provided that the edges
exist. In (a) and (b), the left and right axes correspond to motif-based and
edge-based spectral clustering, respectively.

Next, we provide additional experiments to show that motif-based spectral
clustering has the advantage in wide scenarios where the theoretical
requirements are not necessarily met. We use Adjusted Rand Index (ARI)
(Manning et al., 2010), Normalized Mutual Information (NMI) (Manning et al.,
2010), and Modularity (Newman, 2006) to justify the clustering performance of
the two methods, where larger scores indicate better clustering performance.
In particular, ARI and NMI are computed by comparing the estimated communities
of each of the methods with the true underlying communities in SBMs,
respectively. Modularity is computed based on the _weighted motif adjacency
matrix_ , which measures the difference between the strength of edges (i.e.,
the weighted edges in the weighted motif adjacency matrix) between any two
nodes in the same community and the expected strength of edges between them.
As explained in Remark 2, we expect that motif-based spectral clustering would
lead to larger modularity. In particular, we conduct the following
experiments.

#### Core-periphery structure.

The previous experiments focused on recovering the _affinity_ structure in
networks. That is, the two communities have relatively high within-community
link probability compared to the between-community link probability. We here
assume that the network exhibit _core-periphery_ structure, that is, one of
the two communities has a relatively high within-community link probability
compared to both the other cluster’s within-community link probability and the
between-community link probability (Priebe et al., 2019). Specifically, we
assume the connectivity matrix $B:=[0.5,0.3;0.3;0.05]$ and let the sample size
$n$ vary. Figure 5 shows that in our set-up, the motif-based method has the
advantage over the edge-based counterpart in terms of all three metrics, which
corroborates the ability of the motif-based method for finding core-periphery
structure (Benson et al., 2016).

(a) ARI (b) NMI (c) Modularity

Figure 5: The average performance of two methods on networks with core-
periphery structure over 100 replications. The parameters are
$B:=[0.5,0.3;0.3;0.05]$, $K=2,n=60$. The weights are _i.i.d._ generated from
${\rm Uniform}(0.5,0.6)$ provided that the edges exist.

#### Disassortative weights.

In theoretical analysis, we assumed that all the weights are _i.i.d._
regardless of the communities. We assume that the weights are _disassortative_
in that the within-community edges have less weight than the between-community
edges. In particular, the between-community weights are _i.i.d._ ${\rm
Uniform}(0.5,1)$, and the within-community weights are _i.i.d._ ${\rm
Uniform}(0.01,0.5)$. While the communities are _assortative_ in that the
within-community edges have a larger link probability than between-community
edges. It is obvious that this structure brings difficulty to the estimation
of communities. Figure 6 shows that the motif-based method remains to be
better than the edge-based method in this regime. It is worth mentioning that
when the weights and communities are both assortative, the edge-based method
is good enough to recover the communities since the signal is strong enough
then.

(a) ARI (b) NMI (c) Modularity

Figure 6: The average performance of two methods on networks with
disassortative weights over 100 replications. The parameters are
$K=2,n=60,p=0.5$. The between-community weights are _i.i.d._ ${\rm
Uniform}(0.5,1)$, and the within-community weights are _i.i.d._ ${\rm
Uniform}(0.01,0.5)$.

#### Other weight distributions.

The aforementioned experiments involved uniformly distributed network weights.
Here, we test the efficacy of the motif-based method on two other
distributions. The first is the chi-squared distribution with degree of
freedom 1, denoted by $\chi^{2}(1)$. The second is the exponential
distribution with mean 1, denoted by ${\rm exp}(1)$. The corresponding results
are shown in Figure 7. It can be seen that the motif-based spectral clustering
performs better than the edge-based counterpart, which shows that the efficacy
of the motif-based method is not restrictive in terms of the distribution of
weights.

(I) Weights from $\chi^{2}(1)$

(II) Weights from ${\rm exp}(1)$

(a) ARI (b) NMI

(c) Modularity (d) ARI

(e) NMI (f) Modularity

Figure 7: The average performance of two methods on networks with two kinds of
weight distributions. (I) and (II) correspond to the results with the network
weights _i.i.d._ generated from $\chi^{2}(1)$ and ${\rm exp}(1)$,
respectively. Other parameters are $n=60,p=0.3,K=2$.

#### Other motifs.

From both the theoretical and algorithmic sides, we have mainly focused on the
triangular motif in motif-based spectral clustering. It is of interest to see
whether motif-based spectral clustering remains powerful on other motifs. We
here test the algorithm on two other motifs. The first is the wedge motif,
which is the subgraph with three nodes and two edges. The second is the four-
nodes clique motif, which is the totally connected subgraph with four nodes
and six edges. The weighted motif adjacency matrix is then formulated by
adding the weights in the motif that the node pairs participate in. Figure 8
shows the results. It turns out that motif-based spectral clustering is also
effective on these two motifs. Hence, it would be interesting to theoretically
study the consistency.

(I) Wedge motif

(II) Four-nodes clique

(a) ARI (b) NMI (c) Modularity (d) ARI (e) NMI (f) Modularity

Figure 8: The performance of two methods, where the motif-based spectral
clustering is based on (I) the wedge motif and (II) the four-nodes clique,
respectively. The parameters are $K=2,n=60,p=0.5$, and the edge weights are
_i.i.d._ $\chi^{2}(1)$.

To sum up, the above experiments show that the motif-based spectral clustering
has the advantage over the edge-based counterpart when the clustering problem
is relatively hard, namely, the signal of communities is weak. And the
applicable range of motif-based spectral clustering is far beyond the limits
of theory. It is also worth noting that by our simulated examples, the
modularity is a good measure to judge between the two methods when the true
clustering labels are not known.

##  6 Real data analysis

In this section, we show the merits of higher-order spectral clustering using
three real datasets, including the statisticians’ citation network,
statisticians’ coauthor network, and the ${\rm PM}_{2.5}$ data set. In the
sequel, we introduce the datasets and their corresponding clustering results,
respectively. To save space, all the clustering results are included in the
Supplementary Material.

Statisticians citation network. This dataset was initially collected by Ji and
Jin (2016) based on all published papers from 2003 to the first half of 2012
in four of the top statistical journal: Annals of Statistics, Biometrika,
Journal of American Statistical Association and Journal of Royal Statistical
Society (Series B), which results in 3248 papers and 3607 authors in total. We
here study a citation network, where each node represents an author and the
number of edges between any pair of nodes is equal to the number of one-way
citations between the two authors. We consider the largest connected
components of this network which includes 2654 nodes. We set $K=3$ as in Ji
and Jin (2016). Figure 10 and 11 shows the clustering results for the higher-
order and the edge-based spectral clustering, respectively. We can see from
Figure 10(a) and 11(a) that both of the two methods find two small communities
and one large community, which can be regarded as the background community.
Hence, we take a closer look at the two small communities. As shown in Figure
10(b) and 10(c), the motif-based spectral clustering detects two communities
with many triangles. One community (see Figure 10(b)) consists of
statisticians in high-dimensional statistics, including but without limiting
to the authors of the lasso, group lasso, adaptive lasso, SCAD, graphical
models, which are pioneering works of high-dimensional statistics in the past
20 years. The other community (10(c)) consists of 5 statisticians engaged in a
functional analysis or non-parametric statistics. Turning to the edge-based
spectral clustering, we find that one community (see Figure 11(c)) includes 7
statisticians engaged in functional analysis and non-parametric statistics,
which is very similar to the non-parametric statistics community found by the
motif-based clustering because only two more statisticians are included. After
examining the other community (see Figure 11(b)) carefully, we find that this
community is comprised of two statistician communities, namely, high-
dimensional statisticians and Bayesian statisticians, and statistician Michael
Jordan bridges these two communities. The results are quite interesting, but
from the clustering point of view, motif-based spectral clustering leads to
better results since the resulting communities are purer.

Statisticians coauthor network. This network was also generated based on the
aforementioned dataset (Ji and Jin, 2016). In particular, each node represents
an author and the number of edges between any pair of nodes equals to the
number of papers they coauthored. We consider the largest connected components
of the network which result in 2263 nodes. We also set $K=3$ as in Ji and Jin
(2016). Figure 12 and 13 display the clustering results for the motif-based
and edge-based spectral clustering, respectively. We can see clearly that both
methods detect two main communities, and the remaining communities are so
large that we regard them as the background. In particular, the motif-based
spectral clustering detected two communities (see Figure 12(b) and 12(c)),
which include statisticians engaged in biostatistics/medical statistics and
Bayesian statistics, respectively. By contrast, the edge-based spectral
clustering detects one community (see Figure 13(b)) with 3 biostatisticians
who have worked or studied at Harvard. And the other community (see Figure
13(c)) includes two disconnected components, where one consists of
statisticians in biostatistics/medical statistics, and the other includes
statisticians in non-parametric statistics. Hence, motif-based spectral
clustering detects more reasonable communities than edge-based ones.

${{\rm\bf PM}_{\bf 2.5}}$ pollution data. We collected the ${\rm PM}_{2.5}$
data which consists of daily averaged ${\rm PM}_{2.5}$ concentration in the
year 2015 for each of the 31 Chinese capital cities. In more detail, for each
city, we have its daily averaged ${\rm PM}_{2.5}$ concentrations for 354 days
except 11 days containing missing data. We aim to study the communities of
cities in terms of the ${\rm PM}_{2.5}$ pollution. For this purpose, we should
construct the ${\rm PM}_{2.5}$ pollution network first. Graphical models
provide a useful tool for constructing the network from such data.
Specifically, we treat the ${\rm PM}_{2.5}$ concentration of each city as a
random variable (i.e., the node in the resulting network) and employ the
graphical lasso (Yuan and Lin, 2007; Friedman et al., 2008) with the tuning
parameter selected by eBIC (Foygel and Drton, 2010) to obtain a weighted
network, whereby the rationality of graphical models, the absolute weight
between any pair of nodes (cities) is proportional to the conditional
correlation of the two corresponding random variables, given the remaining
variables. We set $K=5$, which is consistent with the output by the fast
greedy modularity optimization algorithm (Clauset et al., 2004) that decides
the number of communities automatically. Figure 14 and 15 show the clustering
results corresponding to the motif-based and the edge-based spectral
clustering, where for ease of interpretation, we show the cities on the
Chinese map. It can be seen from Figure 14 and 15 that, both methods detect
communities of cities such that cities within each community are closely
located. It could be explained by the fact that the cities located closely
share similar meteorological, economic, and industrial patterns, thus leading
to the conditional dependence of their ${\rm PM}_{2.5}$ pollution. However,
the communities detected by the two methods are different in some sense. There
are two main differences. The first difference is that two communities
detected by the motif-based method in Figure 14(c) and 14(f) are generally
merged into one single community (see Figure 15(e)) by the edge-based method.
The second difference is that the community found by the motif-based method in
Figure 14(d) is generally divided into two communities by the edge-based
method (see Figure 15(d) and 15(f)). One who has common sense about Chinese
geography will feel that the communities detected by the motif-based method
are more reasonable. As an explanation, we conduct a two-sample t-test to test
whether the average ${\rm PM}_{2.5}$ value of cities within the two
communities in Figure 14(c) and 14(f) (motif-based communities) have a
significant difference. The answer is yes, with the p-value smaller than
$2.2\times 10^{-16}$. Similarly, we also test the difference between the two
communities in Figure 15(d) and 15(f) (edge-based communities) in terms of the
mean for ${\rm PM}_{2.5}$. Whereas the answer is no, with the p-value being
0.819. The above analysis indicates that the motif-based method yields more
reasonable communities than the edge-based method does.

#### Comparison of modularity.

The aforementioned analysis shows that the motif-based spectral clustering can
yield more interpretable communities than the edge-based counterpart, which is
consistent with the results in Benson et al. (2016). As a complement, we also
compute the modularity of the communities obtained by both methods; see
Section 5 for details of modularity. Figure 9 displays the box plot of the
modularity of two methods on three datasets over 100 replications. It turns
out that the motif-based method has slightly larger modularity than the edge-
based counterpart on the latter two datasets. While it is slightly worse than
the edge-based one on the first dataset. In addition, it is obvious that the
motif-based method is far more stable than the edge-based method, which makes
the motif-based method more desirable in real applications.

(a) Citation network (b) Coauthor network (c) Pollution network

Figure 9: The box plot of modularity of two methods on three datasets over 100
replications. The modularity is computed based on the weighted motif adjacency
matrix.

##  7 Conclusion

The higher-order structure of networks and the corresponding higher-order
spectral clustering have been proven to be insightful and helpful in real
network clustering problems. However, there are rare works disclosing the
merits of higher-order spectral clustering over the edge-based one
systematically and theoretically. Motivated by this, we theoretically studied
the clustering performance of the higher-order spectral clustering from a
statistical perspective, where we typically assumed the underlying networks
are generated from a weighted stochastic block model. The results showed that
when the network is dense and the signal of weights is weak, the higher-order
spectral clustering could produce better misclustering upper bound than the
edge-based counterpart. The methodologies and theories in this paper can be
generally used in other problems. For example, the weighted motif adjacency
matrix can be thought of as a network denoising trick and can be used for any
downstream network analysis. In addition, the theoretical tools in this paper
can be used to analyze networks with summary statistics but without edge
information, which is of independent interest in the context of privacy-
preserving network analysis.

There are many ways that the content of this work could be extended. First, we
only specify the expectation of weights in the WSBM. It would be beneficial to
derive more delicate results if the distribution of weights is incorporated.
Second, we found by experiments that higher-order spectral clustering also
works better than edge-based under the degree-corrected stochastic block
models. It would be interesting to understand this phenomenon theoretically,
which is currently under our consideration. Third, we considered triangles in
undirected networks. More types of higher-order motifs and the clustering in
directed networks deserve further research. Besides, it is of interest to see
whether the motif-based method can be generalized to the clustering of multi-
layer networks (Paul and Chen, 2020). Finally, we take a step forward in
understanding how higher-order spectral clustering works, and there is still
an urgent need to understand it thoroughly. In addition, it is crucial to
develop a network-generating model that really captures the higher-order
connectivity patterns in real networks.

## Acknowledgement

This work is partially done when Xiao Guo was a visiting student at Columbia
University. The authors would like to thank Professor Ming Yuan for
enlightening discussions in the early stage of this project. Our research is
partially supported by National Natural Science Foundation for Outstanding
Young Scholars (No.72122018), National Natural Science Foundation of China
(No.U1811461), and Natural Science Foundation of Shaanxi Province (No.
2021JQ-429 and No.2021JC-01).

SUPPLEMENTARY MATERIAL

This file includes the proofs for theorems and the figures on clustering
results of real data analysis.

## A. Proof of main results

### Proof of Lemma 1

Let $\Delta={\rm diag}(\sqrt{n_{1}},...,\sqrt{n_{K}})$ and recall (2.6), then
we have

$\displaystyle\mathcal{W}^{M}$
$\displaystyle=\Theta((h_{1}-h_{2})I_{K}+h_{2}1_{K}1_{K}^{\intercal})\Theta^{\intercal}$
$\displaystyle=\Theta\Delta^{-1}\Delta((h_{1}-h_{2})I_{K}+h_{2}1_{K}1_{K}^{\intercal})\Delta(\Theta\Delta^{-1})^{\intercal},$

where $h_{1}>h_{2}>0$. Recalling the structure of $\Theta$, one can easily
show $\Theta\Delta^{-1}$ is orthonormal. Denote the eigen-decomposition of
$\Delta(h_{1}-h_{2})I_{K}+h_{2}1_{K}1_{K}^{\intercal}\Delta$ as

$Z\Sigma
Z^{\intercal}=\Delta((h_{1}-h_{2})I_{K}+h_{2}1_{K}1_{K}^{\intercal})\Delta.$

Then we have

$\mathcal{W}^{M}=\Theta\Delta^{-1}Z\Sigma(\Theta\Delta^{-1}Z)^{\intercal}=U\Sigma
U^{\intercal}.$

Thus $U=\Theta\Delta^{-1}Z=\Theta X$ with $X:=\Delta^{-1}Z$ as
$\Theta\Delta^{-1}$ and $Z$ are orthonormal, respectively. Moreover, as $Z$ is
orthonormal and $\Delta^{-1}Z$ is a square matrix, we can verify that the rows
of $\Delta^{-1}Z$ are perpendicular to each other and the $k$th row has length
$\sqrt{1/n_{k}}$, which implies
$\|X_{k\ast}-X_{l\ast}\|_{2}=\sqrt{n_{k}^{-1}+n_{l}^{-1}}$. $\blacksquare$

### Proof of Theorem 1

The general framework of the following proof is adapted from Feige and Ofek
(2005), whose arguments are also used in Lei and Rinaldo (2015); Paul et al.
(2018); Chin et al. (2015); Gao et al. (2017), among others. The proof
includes three major steps.
Step 1: Discretization. We first reduce controlling
$\|W^{M}-\mathbb{E}(W^{M})\|_{2}$ to the problem of bounding the supremum of
$|x^{\intercal}(W^{M}-\mathbb{E}(W^{M}))y|$ over all pairs of vectors $x,y$ in
a finite set of grid points. For any given pair $(x,y)$ in the grid, the
quantity $x^{\intercal}(W^{M}-\mathbb{E}(W^{M}))y$ is decomposed into the sum
of two parts. The first part corresponds to the small entries of both $x$ and
$y$, called the light pairs, the other part corresponds to large entries of
$x$ or $y$, called the heavy pairs.
Step 2: Bounding the light pairs. Since the elements of $W^{M}$ are dependent,
we can not use the Bernstein’s inequality to control the contribution of the
light pairs. Instead, we make use of the typical bounded differences
inequality established in Warnke (2016) to bound the light pairs. Note that
Paul et al. (2018) also used similar arguments to do this in the unweighted
SBMs, but the details are different from those we do in the weighted case.
Step 3: Bounding the heavy pairs. The contribution from heavy pairs will be
bounded using a combinatorial argument on the event that the edge weights
(w.r.t. the weighted motif adjacency matrix $W^{M}$) in a collection of
subgraphs do not deviate much from their expectation. To that end, the
concentration inequality in Warnke (2017) will be used to bound the edge
weighted degree.

It should be noted that throughout the proof, we will use $c$, $c^{\prime}$,
$c^{\prime\prime}$ likewise to denote constants which may be different from
line to line. Now we proceed to prove Theorem 1. Let $S$ denote the unit ball
in the $n$-dimensional Euclidean space. Fix $\epsilon\in(0,1)$, for example
$\epsilon=\frac{1}{2}$, define an $\epsilon$-net of the ball as follows:

$\mathcal{N}=\\{x=(x_{1},...,x_{n})\in S:\forall
i,\sqrt{n}x_{i}/\epsilon\in\mathbb{Z}\\},$ (A.1)

where $\mathbb{Z}$ denotes the set of all integers. Thus $\mathcal{N}$
consists of all grid points of size $\epsilon\sqrt{n}$ within the unit ball.
The Lemma 2.1 in Lei and Rinaldo (2015) shows that for any
$A\in\mathbb{R}^{n\times n}$,

$\|A\|_{2}\leq(1-\epsilon)^{-2}\underset{x,y\in\mathcal{N}}{{\rm
sup}}|x^{\intercal}Ay|.$ (A.2)

Therefore, to bound $\|W^{M}-\mathbb{E}(W^{M})\|_{2}$, we only need to bound
$|x^{\intercal}(W^{M}-\mathbb{E}(W^{M}))y|$ over all possible pairs
$(x,y)\in\mathcal{N}$. For any pair of vector $(x,y)$, we have

$x^{\intercal}(W^{M}-\mathbb{E}(W^{M}))y=\sum_{i,j}x_{i}y_{j}(W^{M}-\mathbb{E}(W^{M}))_{ij}.$
(A.3)

We split the pairs $(x_{i},y_{j})$ into _light pairs_

${\mathcal{L}=\mathcal{L}(x,y)=\\{(i,j):|x_{i}y_{j}|\leq\frac{\sqrt{D}}{n\tau_{{\rm
max}}}\\},}$

and into _heavy pairs_

${\mathcal{H}=\mathcal{H}(x,y)=\\{(i,j):|x_{i}y_{j}|>\frac{\sqrt{D}}{n\tau_{{\rm
max}}}\\},}$

where $D$ and $\tau_{{\rm max}}$ are defined in Theorem 1. For the light
pairs, we have the following results.

######  Lemma 2 (_Light pairs_)

For some constant $r_{1}>0$, there exists a constant $c_{1}(r_{1})>0$ such
that with probability at least $1-n^{-r_{1}}$,

$\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}\underset{(i,j)\in\mathcal{L}}{\sum}x_{i}y_{j}(W^{M}-\mathbb{E}(W^{M}))_{ij}\big{|}\leq
c_{1}\sqrt{D}.$ (A.4)

To bound the contribution of heavy pairs, we first have the following
observation,

$\displaystyle\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}\underset{(i,j)\in\mathcal{H}}{\sum}x_{i}y_{j}(W^{M}-\mathbb{E}(W^{M}))_{ij}\big{|}\leq\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}\underset{(i,j)\in\mathcal{H}}{\sum}x_{i}y_{j}W^{M}_{ij}\big{|}+\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}\underset{(i,j)\in\mathcal{H}}{\sum}x_{i}y_{j}\mathbb{E}(W^{M}_{ij})\big{|}.$
(A.5)

Recall that $\mathbb{E}(W_{ij})\leq\alpha_{n}p_{n}$, then

$\displaystyle\mathbb{E}(W^{M}_{ij})$
$\displaystyle=\mathbb{E}\sum_{k}\mathbf{1}(W_{ij}\cdot W_{ik}\cdot
W_{jk}>0)(W_{ij}+W_{ik}+W_{jk})\leq 3n\alpha_{n}p_{n}^{3}.$ (A.6)

Then the second term in (A.5) can be bounded as follows,

$\displaystyle\big{|}\underset{(i,j)\in\mathcal{H}}{\sum}x_{i}y_{j}\mathbb{E}(W^{M}_{ij})\big{|}$
$\displaystyle\leq\underset{(i,j)\in\mathcal{H}}{\sum}\frac{x_{i}^{2}y_{j}^{2}}{|x_{i}y_{j}|}\mathbb{E}(W^{M}_{ij})$
$\displaystyle\leq\frac{n\tau_{{\rm max}}}{\sqrt{D}}{\rm
max}_{ij}\mathbb{E}(W^{M}_{ij})\sum_{(i,j)}x_{i}^{2}y_{j}^{2}$
$\displaystyle\leq{\frac{n\tau_{{\rm
max}}}{\sqrt{D}}3n\alpha_{n}p_{n}^{3}\leq\frac{cD}{\sqrt{D}}\leq c\sqrt{D}},$
(A.7)

where the second inequality follows from the definition of heavy pairs
$\mathcal{H}$, the third inequality follows from (A.6) and the vectors $x$ and
$y$ are within the unit ball, and the penultimate inequality follows since

$n\tau_{{\rm max}}\cdot 3n\alpha_{n}p_{n}^{3}\leq
cn^{3}\alpha_{n}p_{n}^{5}\leq cD,$ (A.8)

by recalling $\tau_{\rm max}=np_{n}^{2}$. For the first term in (A.5), we have
the following result.

######  Lemma 3 (_Heavy pairs_)

For some constants $r_{2}>0$, $r_{3}>0$, and any $0<\beta<1$, there exists a
constant $c_{2}(r_{2},r_{3},\beta)>0$ such that with probability at least
$1-n^{-r_{2}}-{\rm exp}(-r_{3}({\rm log}n)^{\beta})$,

$\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}\underset{(i,j)\in\mathcal{H}}{\sum}x_{i}y_{j}W^{M}_{ij}\big{|}\leq
c_{2}\sqrt{D}.$ (A.9)

Combining the results for the light and heavy pairs and recalling (A.2), we
obtain

$\displaystyle\|W^{M}-\mathbb{E}(W^{M})\|_{2}$ $\displaystyle\leq
4\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}{\sum}x_{i}y_{j}(W^{M}-\mathbb{E}(W^{M}))_{ij}\big{|}\leq
c\sqrt{D}.$ (A.10)

$\blacksquare$

### Proof of Theorem 2

We make use of the framework in Lei and Rinaldo (2015) to bound the
misclustering error rate. To fix ideas, we give some notation now. $U$ and
$\tilde{U}$ denote the $K$ leading eigenvectors of $\mathcal{W}^{M}$ and
$W^{M}$, respectively. $\hat{{U}}:=\hat{{\Theta}}\hat{{Y}}$ corresponds to the
optimal solution of the higher-order spectral clustering algorithm in that,

$(\hat{{\Theta}},\hat{{Y}})=\underset{{\Theta\in\mathbb{M}_{n,K},Y\in\mathbb{R}^{K\times
K}}}{{\rm arg\;min}}\;\|\Theta Y-\tilde{U}\|_{\rm\tiny F}^{2}.$

By Lemma 1, in the WSBMs, two nodes are in the same community if and only if
the corresponding rows of the population eigenvector matrix $U$ are the same.
Building on this result, in what follows, we first bound the deviation of
$\hat{U}$ from $U$. Then, within each true cluster, we bound the size of nodes
that correspond to a large deviation of $\hat{U}$ from $U$, we bound their
size. After that, we show for the remaining nodes that the estimated and true
communities are consistent.

First, we bound the deviation of $\hat{U}$ from $U$. Davis-Kahan $\rm
sin\Theta$ theorem (Theorem VII.3.1 of Bhatia (1997)) provides a useful
technical tool for bounding the perturbation of eigenvectors from the
perturbation of matrices. In particular, by Theorem 3.1 of Lei and Rinaldo
(2015), there exists a $K\times K$ orthogonal matrix $O$ such that,

$\|\tilde{U}-UO\|_{\tiny\rm
F}\leq\frac{2\sqrt{2K}}{\lambda_{K}(\mathcal{W}^{M})}\|W^{M}-\mathcal{W}^{M}\|_{2},$
(A.11)

where we recall that $\lambda_{K}(\mathcal{W}^{M})$ denotes the minimum non-
zero eigenvalue of the population matrix $\mathcal{W}^{M}$. Now we proceed to
bound the deviation of $\hat{{U}}$ from $U$. Note that

$\displaystyle\|\hat{U}-UO\|_{\tiny\rm F}^{2}$
$\displaystyle=\|\hat{{U}}-\tilde{U}+\tilde{U}-UO\|_{\tiny\rm F}^{2}$
$\displaystyle\leq\|UO-\tilde{U}\|_{\tiny\rm F}^{2}+\|\tilde{U}-UO\|_{\tiny\rm
F}^{2}$ $\displaystyle=2\|\tilde{U}-UO\|_{\tiny\rm F}^{2},$ (A.12)

where the first inequality follows from our assumption that $\hat{{U}}$ is the
global solution of the higher-order spectral clustering algorithm and $UO$ is
a feasible solution. Then combine (Proof of Theorem 2) with (A.11), we have

$\displaystyle\|\hat{{U}}-UO\|_{\tiny\rm
F}^{2}\leq\frac{16K}{\lambda_{K}^{2}(\mathcal{W}^{M})}\|W^{M}-\mathcal{W}^{M}\|_{2}^{2}.$
(A.13)

Now we calculate the terms on the RHS of (A.13). Recall the definition of
$\mathcal{W}^{M}$ in (2.6) and the definitions of $h_{1}$ and $h_{2}$ in (2.4)
and (2.5), we can easily obtain

$\lambda_{K}(\mathcal{W}^{M})=\frac{n}{K}(h_{1}-h_{2})\asymp\frac{n^{2}\alpha_{n}p_{n}^{3}}{K},$
(A.14)

and

$\displaystyle\|W^{M}-\mathcal{W}^{M}\|_{2}$
$\displaystyle=\|W^{M}-\mathbb{E}(W^{M})-{\rm
diag}(\mathcal{W}^{M})\|_{2}\asymp\|W^{M}-\mathbb{E}(W^{M})\|_{2}.$ (A.15)

Combining (A.14), (A.15), and the bound in Theorem 1 with (A.13), we have

$\displaystyle\|\hat{{U}}-UO\|_{\tiny\rm
F}^{2}\leq\frac{c16K^{3}}{n\alpha_{n}p_{n}}:={\rm
err}(K,n,c,\alpha_{n},p_{n}),$ (A.16)

where for notational convenience, we use ${\rm err}(K,n,c,\alpha_{n},p_{n})$
to denote ${c16K^{3}}/{n\alpha_{n}p_{n}}$ in what follows.

In the sequel, we proceed to bound the fraction of misclustered nodes within
each true cluster. By Lemma 1, we can write $U=\Theta X$, where
$\|X_{k\ast}-X_{l\ast}\|_{2}=\sqrt{n_{k}^{-1}+n_{l}^{-1}}$ for all $1\leq
k<l\leq K$. Hence $UO=\Theta XO=\Theta X^{\prime}$ with $X^{\prime}=XO$, and
$\|X^{\prime}_{k\ast}-X^{\prime}_{l\ast}\|_{2}=\sqrt{n_{k}^{-1}+n_{l}^{-1}}$
by the orthogonality of $O$. Define

$\delta_{k}={\rm min}_{l\neq
k}\;\|X^{\prime}_{k\ast}-X^{\prime}_{l\ast}\|_{2}=\sqrt{\frac{1}{n_{k}}+\frac{1}{{\rm
max}\\{n_{l}:l\neq k\\}}},$ (A.17)

and

$S_{k}=\\{i\in G_{k}(\Theta):\;\|(\hat{U})-(UO)_{i\ast}\|_{\tiny\rm
F}>\frac{\delta_{k}}{2}\\},$ (A.18)

where $S_{k}$ is essentially the number of misclustered nodes in the true
cluster $k$ (after some permutation) as we will see soon. By the definition of
$S_{k}$, it is obvious to see

$\sum_{k=1}^{K}|S_{k}|\delta_{k}^{2}/4\leq\|\hat{{U}}-UO\|_{\tiny\rm
F}^{2}={\rm err}(K,n,c,\alpha_{n},p_{n}).$ (A.19)

Recall that $\delta_{k}={\rm min}_{l\neq
k}\;\|X^{\prime}_{k\ast}-X^{\prime}_{l\ast}\|_{2}=\sqrt{\frac{1}{n_{k}}+\frac{1}{{\rm
max}\\{n_{l}:l\neq k\\}}}$, so $n_{k}\delta_{k}^{2}\geq 1$. Therefore, (A.19)
entails that

$\sum_{k=1}^{K}\frac{|S_{k}|}{n_{k}}\leq 4\cdot\|\hat{{U}}-UO\|_{\tiny\rm
F}^{2}=4\cdot{\rm err}(K,n,c,\alpha_{n},p_{n}).$ (A.20)

At last, we show that the nodes within in true community $k$ (after some
permutation) but outside $S_{k}$ are correctly clustered. Before moving on, we
first prove $|S_{k}|<n_{k}$. We have by (A.19) that

$|S_{k}|\leq\frac{4}{\delta_{k}^{2}}{\rm err}(K,n,c,\alpha_{n},p_{n}).$ (A.21)

As $n_{k}\delta_{k}^{2}\geq 1$, it suffices to prove

$4\cdot{\rm err}(K,n,c,\alpha_{n},p_{n})<1.$ (A.22)

which actually follows from the assumption (3.3) after some modification of
$c$ ($c$ can be different from line to line). As a result, we have
$|S_{k}|<n_{k}$ for every $1\leq k\leq K$. And thus, $T_{k}\equiv
G_{k}\backslash S_{k}\neq\emptyset$, where we recall that $G_{k}$ denotes the
nodes in the true community $k$. Let $T=\cup_{k=1}^{K}T_{k}$, we now show that
the rows in $(UO)_{T\ast}$ has a one to one correspondence with those in
$\hat{{U}}_{T\ast}$. On the one hand, for $i\in T_{k}$ and $j\in T_{l}$ with
$l\neq k$, $\hat{{U}}_{i\ast}\neq\hat{{U}}_{j\ast}$, otherwise one can have
the following contradiction

$\displaystyle{\rm max}\\{\delta_{k},\delta_{l}\\}$
$\displaystyle\leq\|(UO)_{i\ast}-(UO)_{j\ast}\|_{2}$
$\displaystyle\leq\|(UO)_{i\ast}-\hat{{U}}_{i\ast}\|_{2}+\|(UO)_{j\ast}-\hat{{U}}_{j\ast}\|_{2}$
$\displaystyle<\frac{\delta_{k}}{2}+\frac{\delta_{l}}{2},$ (A.23)

where the first and last inequality follows from (A.17) and (A.18),
respectively. On the other hand, for $i,j\in T_{k}$, we have
$\hat{{U}}_{i\ast}=\hat{{U}}_{j\ast}$, because otherwise $U_{T\ast}$ has more
than $K$ distinct rows contradicting the fact that the output community size
is $K$.

As a result, we have proved the membership is correctly recovered outside of
$\cup_{k=1}^{K}S_{k}$ and the rate of misclustered nodes in $S_{k}$ is bounded
as in (A.20), and therefore the conclusion of Theorem 2 is arrived.
$\blacksquare$

## B. Proof of auxiliary Lemmas

### Proof of Lemma 2

We will make use of the typical bounded differences inequality established in
Warnke (2016) to control the contribution of the light pairs. To fix ideas, we
reproduce the result that we will use in the following proposition.

######  Proposition 1 (Theorem 1.2 of Warnke (2016))

Let $X=(X_{1},...,X_{N})$ be a family of independent random variables with
$X_{k}$ taking values in a set $\Lambda_{k}$. Let
$\Gamma\subseteq\prod_{j\in[N]}\Lambda_{j}$ be an event and assume that
function $f:\prod_{j\in[N]}\Lambda_{j}\rightarrow\mathbb{R}$ satisfies the
following typical Lipschitz condition.
(TL) There are numbers $(c_{k})_{k\in[N]}$ and $(d_{k})_{k\in[N]}$ with
$c_{k}\leq d_{k}$ such that, whenever
$x,\tilde{x}\in\prod_{j\in[N]}\Lambda_{j}$ differ only in the $k$th
coordinate, we have

$\big{|}f(x)-f(\tilde{x})\big{|}\leq\begin{cases}c_{k},&\mbox{if
}\;\Gamma\;\mbox{holds},\\\ d_{k},&otherwise.\end{cases}$ (B.1)

Then for all numbers $(\gamma_{k})_{k\in[N]}$ with $\gamma_{k}\in(0,1]$ there
is an event $\mathcal{B}=\mathcal{B}(\Gamma,(\gamma_{k})_{k\in[N]})$
satisfying

$\mathbb{P}(\mathcal{B})\leq\mathbb{P}(X\notin\Gamma)\underset{k\in[N]}{\sum}{\gamma_{k}^{-1}}\quad{and}\quad\urcorner\mathcal{B}\subseteq\Gamma,$
(B.2)

such that for all $t\geq 0$ we have

$\mathbb{P}(f(X)\geq\mathbb{E}f(X)+t\;and\;\urcorner\mathcal{B})\leq{\rm
exp}\big{(}-\frac{t^{2}}{2\sum_{k\in[N]}(c_{k}+e_{k})^{2}}\big{)},$ (B.3)

where $e_{k}=\gamma_{k}(d_{k}-c_{k}).$ In many applications, the following
simple consequence of (B.2) and (B.3) suffices:

$\displaystyle\mathbb{P}(f(X)\geq\mathbb{E}f(X)+t)\leq{\rm
exp}\big{(}-\frac{t^{2}}{2\sum_{k\in[N]}(c_{k}+\gamma_{k}(d_{k}-c_{k}))^{2}}\big{)}+\mathbb{P}(X\notin\Gamma)\underset{k\in[N]}{\sum}\gamma_{k}^{-1}.$
(B.4)

Now we proceed to bound the light pairs. Define

$u_{ij}=x_{i}y_{j}\mathbf{1}((i,j)\in\mathcal{L})+x_{j}y_{i}\mathbf{1}((j,i)\in\mathcal{L})$

for all $i,j=1,...,n$. Then

$\displaystyle\underset{(i,j)\in\mathcal{L}}{\sum}x_{i}y_{j}(W^{M}-\mathbb{E}(W^{M}))_{ij}$
$\displaystyle=\underset{i<j}{\sum}\underset{k\neq
i,j}{\sum}\mathbf{1}(W_{ij}W_{ik}W_{jk}>0)(W_{ij}+W_{ik}+W_{jk})u_{ij}$
$\displaystyle-\underset{i<j}{\sum}\underset{k\neq
i,j}{\sum}\mathbb{E}\big{(}\mathbf{1}(W_{ij}W_{ik}W_{jk}>0)(W_{ij}+W_{ik}+W_{jk})\big{)}u_{ij}.$
(B.5)

To use the result in Proposition 1, we define

${f(W)=\underset{i<j}{\sum}\underset{k\neq
i,j}{\sum}\mathbf{1}(W_{ij}W_{ik}W_{jk}>0)(W_{ij}+W_{ik}+W_{jk})u_{ij}.}$
(B.6)

It is obvious that $f$ is a function of independent variables
$W_{ij}(i,j=1,...,n)$. Now we proceed to bound the effect on $f$ when only one
element of $W_{ij}$ is changed, namely, we specify (B.1) in our setting.
Suppose $W_{ij}$ changes, then the effect on $f$ may be “large” on the term
involving $u_{ij}$. In particular, the effect can be bounded as

$\underset{k\neq i,j}{\sum}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})u_{ij}.$ (B.7)

On the other hand, the effect on $f$ may be “small” on the term involving
$u_{ik}$ or $u_{jk}$. In particular, the effect can be bounded as $3u_{ik}$ or
$3u_{jk}$, respectively, provided that $k$ is the common neighborhood of $i$
and $j$. To further bound the “large” effect and “small” effect, we define

${\tau_{ij}=\underset{k\neq i,j}{\sum}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})},$ (B.8)

by recalling the formula in (B.7) and define the “good set” $\Gamma$

${\Gamma=\\{(W_{ij}):\underset{ij}{\rm max}\,\tau_{ij}\leq c\tau_{\rm
max}\\},}$ (B.9)

where $c$ is large enough. We want to prove that the probability of $W$
outside the good set vanishes when $n$ goes to infinity. To this end, we first
bound the expectation of

$\delta_{ij}^{k}:=\mathbf{1}(W_{ik}\cdot W_{jk}>0)(1+W_{ik}+W_{jk}).$

We have

$\displaystyle\mathbb{E}(\delta_{ij}^{k})$
$\displaystyle{=\mathbb{E}\big{(}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})\big{)}}\leq c^{\prime}p_{n}^{2}.$ (B.10)

Now we use the Bernstein’s inequality to bound the the probability of $W$
outside the good set $\Gamma$.

$\displaystyle\mathbb{P}(W\notin\Gamma)$ $\displaystyle\leq
n^{2}\mathbb{P}(\tau_{ij}>c\tau_{\rm max})$ $\displaystyle\leq
n^{2}\mathbb{P}\big{(}\underset{k\neq i,j}{\sum}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})>c\tau_{\rm max}\big{)}$ $\displaystyle\leq
n^{2}\mathbb{P}\big{(}\underset{k\neq i,j}{\sum}\big{[}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})-\mathbb{E}(\delta_{ij}^{k})\big{]}>c\tau_{\rm
max}-c^{\prime}np_{\rm max}\big{)}$ $\displaystyle\leq
n^{2}\mathbb{P}\big{(}\underset{k\neq i,j}{\sum}\big{[}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})-\mathbb{E}(\delta_{ij}^{k})\big{]}>c^{\prime\prime}\tau_{\rm
max}\big{)}$ $\displaystyle\leq n^{2}{\rm exp}\big{(}-\frac{c\tau_{\rm
max}^{2}}{c^{\prime}np_{n}^{2}+c^{\prime\prime}\tau_{\rm max}}\big{)}\leq
n^{2}{\rm exp}\big{(}-c^{\prime\prime\prime}\tau_{\rm max}\big{)}$
$\displaystyle\leq n^{2}{\rm exp}(-{\rm log}n^{c})\leq n^{-c+2}\leq
n^{-c^{\prime}},$ (B.11)

where we used the following facts, $\tau_{\rm max}\geq c{\rm log}n$ which is
implied by $\tau_{\rm max}=np_{n}^{2}$ and $p_{n}^{2}\geq c{\rm log}n/n$,
$\mathbb{E}((\delta_{ij}^{k})^{2})\leq c\mathbb{E}(\delta_{ij}^{k})$ and
(B.10).

Next, we continue to bound the effect of changing one element $W_{ij}$ of $W$
on $f(W)$. Before that, we first have the following observation. Recalling
(B.8), under the good event $\Gamma$, we have

$\displaystyle\underset{k\neq i,j}{\sum}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)\leq\tau_{ij}=\underset{k\neq i,j}{\sum}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})\leq c\tau_{\rm max},$ (B.12)

which means that the number of $k$ being the common neighborhood of $i$ and
$j$ is not larger than $c\tau_{\rm max}$ under the good event $\Gamma$. Recall
that, the “large” effect can be bounded as

$\underset{k\neq i,j}{\sum}\mathbf{1}(W_{ik}\cdot
W_{jk}>0)(1+W_{ik}+W_{jk})u_{ij},$

and the “small” effect can be bounded as $3u_{ik}$ or $3u_{jk}$, respectively,
provided that $k$ is the common neighborhood of $i$ and $j$. Therefore,
combining the “large” and “small” effects, under the good event $\Gamma$, the
total effect of changing one element $W_{ij}$ in $f(W)$ can be upper bounded
as

$c_{ij}\leq c\tau_{\rm max}u_{ij}+3\sum_{k:W_{ik}W_{jk}>0}u_{ik}\leq
c\tau_{\rm max}u_{ij}+c^{\prime}\tau_{\rm max}u_{ik}.$ (B.13)

On the other hand, when the bad event $\urcorner\Gamma$ occurs, the total
effect of changing one element $W_{ij}$ in $f(W)$ becomes

$d_{ij}\leq 3nu_{ij}+3nu_{ik}.$ (B.14)

We are now moving towards (B.3) in our setting. Define
$\gamma_{ij}=\frac{1}{n}$ for all $i,j$. Then
$e_{ij}=\gamma_{ij}(d_{ij}-c_{ij})=o(c_{ij})$ by recalling the definition of
$\tau_{\rm max}$. And we have

$\displaystyle\sum_{ij}c_{ij}^{2}\leq c\tau_{\rm max}^{2}\sum_{ij}u_{ij}^{2},$
(B.15)

and

$\displaystyle\mathbb{P}(\mathcal{B})$
$\displaystyle\leq\mathbb{P}(W\notin\Gamma)\underset{i,j}{\sum}\gamma_{ij}^{-1}=n^{3}\mathbb{P}(W\notin\Gamma)\leq{\rm
exp}\big{(}-(c^{\prime}-3){\rm log}n\big{)}.$ (B.16)

Therefore, combining (B.15) with (B.3), we have for large enough $c$ that

$\displaystyle{\mathbb{P}(f(W)-\mathbb{E}f(W)>c\sqrt{D}\cap\mathcal{\urcorner}\mathcal{B})}$
$\displaystyle\leq{\rm exp}\big{(}-\frac{c^{2}D}{c^{\prime}\tau_{\rm
max}^{2}\sum_{ij}u_{ij}^{2}}\big{)}$ $\displaystyle\leq{\rm
exp}\big{(}-\frac{c^{2}\tau_{\rm max}^{2}n\alpha_{n}p_{n}}{c^{\prime}\tau_{\rm
max}^{2}\sum_{ij}u_{ij}^{2}}\big{)}$ $\displaystyle\leq{\rm
exp}(-cn\sqrt{n\alpha_{n}p_{n}})\leq{\rm exp}(-c^{\prime\prime}n),$ (B.17)

where the second inequality follows from the fact that

$\displaystyle\sum_{ij}u_{ij}^{2}$ $\displaystyle\leq c\,{\rm
max}_{ij}|u_{ij}|\sum_{i,j}|x_{i}y_{j}|\leq c\frac{\sqrt{D}}{n\tau_{\rm
max}}\|x\|_{2}\|y\|_{2}\leq c^{\prime}\frac{\sqrt{D}}{n\tau_{\rm max}},$

and the last inequality is implied by $\alpha_{n}p_{n}>c^{\prime\prime}({\rm
log}n)^{\beta}/n$ for some small constant $c^{\prime\prime}$. Finally, it is
obvious that the event $\\{W\notin\Gamma\\}$ does not depend on the choice of
vectors $x$ and $y$, thus taking the supremum over all $x$ and $y$, and
combining (B.16) and (Proof of Lemma 2) with (B.4), we obtain

$\displaystyle\mathbb{P}(\underset{x,y\in\mathcal{N}}{{\rm
sup}}\big{|}\underset{(i,j)\in\mathcal{L}}{\sum}x_{i}y_{j}(W^{M}-\mathbb{E}(W^{M}))_{ij}\big{|}\geq
c_{1}\sqrt{D})\leq{\rm exp}\big{(}-(c^{\prime\prime}-{\rm log}5)n\big{)}+{\rm
exp}\big{(}-(c^{\prime}-3){\rm log}n\big{)},$ (B.18)

where we used the fact that $|\mathcal{N}|\leq{\rm exp}(n{\rm log}5)$. Since
$c^{\prime\prime}$ and $c^{\prime}$ are large enough, the probability in
(B.18) decays polynomially with $n$. $\blacksquare$

### Proof of Lemma 3

The proof follows a similar strategy as in Feige and Ofek (2005); Paul et al.
(2018); Lei and Rinaldo (2015), among others. We will focus on the heavy pairs
$(i,j)$ such that $x_{i}>0,y_{j}>0$ and the set

$\mathcal{H}_{1}=\\{(i,j)\in\mathcal{H}:x_{i}>0,y_{j}>0\\}.$

The other three cases can be treated similarly. We will need the following
notation:

  * $\bullet$

$I_{1}=\\{i:\frac{1}{2n^{1/2}}\leq x_{i}\leq\frac{1}{n^{1/2}}\\}$,
$I_{s}=\\{i:\frac{2^{s-1}}{2n^{1/2}}\leq x_{i}\leq\frac{2^{s}}{2n^{1/2}}\\}$
for $s=2,3,...,\lceil{\rm log}_{2}2n^{1/2}\rceil$.
$J_{1}=\\{j:\frac{1}{2n^{1/2}}\leq y_{j}\leq\frac{1}{n^{1/2}}\\}$,
$J_{t}=\\{j:\frac{2^{t-1}}{2n^{1/2}}\leq y_{j}\leq\frac{2^{t}}{2n^{1/2}}\\}$
for $t=2,3,...,\lceil{\rm log}_{2}2n^{1/2}\rceil$.

  * $\bullet$

$e(I,J)$ denotes the weights of distinct edges in the motif-based network
$W^{M}$ between nodes sets $I$ and $J$. $e(I,J)=\sum_{i\in I}\sum_{j\in
J}W^{M}_{ij},\,\mbox{if }\,I\cap J=\emptyset$; $\sum_{(i,j)\in I\times
J\backslash(I\cap J)^{2}}W^{M}_{ij}+\sum_{(i,j)\in(I\cap
J)^{2},i<j}W^{M}_{ij},\,\mbox{if}\;I\cap J\neq\emptyset$.

  * $\bullet$

$\mu(I,J)=\mathbb{E}(e(I,J))$, and $\bar{\mu}=|I||J|\frac{\Delta}{n}$, where
we define $\Delta=n^{2}\alpha_{n}p_{n}^{3}.$

  * $\bullet$

$\lambda_{st}=e(I_{s},J_{t})/\bar{\mu}_{st}$, where
$\bar{\mu}_{st}=\bar{\mu}(I_{s},J_{t})$.

  * $\bullet$

$\alpha_{s}=|I_{s}|2^{2s}/n,\beta_{t}=|J_{t}|2^{2t}/n$, and
$\sigma_{st}=\lambda_{st}\frac{\sqrt{D}}{\tau_{\rm max}}2^{-(s+t)}$.

We begin by providing the following two results which play a key role in the
proof.

######  Lemma 4 (Bounded degree)

Let

$d_{i}^{M}=\sum_{j}\sum_{k\neq i,j}(\mathbf{1}(W_{ij}\cdot W_{ik}\cdot
W_{jk}>0)(W_{ij}+W_{ik}+W_{jk}))$

denote the weighted triangle degree of node $i$. If
$\alpha_{n}p_{n}^{2}>c^{\prime\prime}({\rm log}n)^{\beta}/n$ for any
$0<\beta<1$, then for all $i$ and some constant $r_{4},r_{5}>0$, there exists
a constant $c_{3}(r_{4},r_{5})>0$ such that

$d_{i}^{M}\leq c_{3}\Delta,$

with probability larger than $1-n^{-r_{4}}-{\rm exp}(-r_{5}({\rm
log}n)^{\beta})$.

######  Lemma 5 (Bounded discrepancy)

For a constant $r_{6}>0$, there exists constants $c_{4}=c_{4}(r_{6})>0$,
$c_{5}=c_{5}(r_{6})>0$ such that with probability larger than $1-2n^{-r_{6}}$,
for any $I,J\subseteq[n]$ with $[I]\leq[J]$, at least one of the following
statements holds:

$(a)\;\frac{e(I,J)}{\bar{\mu}(I,J)}\leq ec_{4}$,

$(b)\;e(I,J){\rm log}\frac{e(I,J)}{\bar{\mu}(I,J)}\leq c_{5}\tau_{\rm
max}|J|{\rm log}\frac{n}{|J|}$.

Define $\mathcal{E}$ be the event that the results of Lemma 4 and 5 hold. In
the sequel, we proceed to bound the heavy pairs under the event $\mathcal{E}$.
We first have the following facts.

$\displaystyle\underset{(i,j)\in\mathcal{H}_{1}}{\sum}x_{i}y_{j}W^{M}_{ij}$
$\displaystyle\leq 2\underset{2^{(s+t)}\geq\frac{\sqrt{D}}{\tau_{\rm
max}}}{\sum}e(I_{s},J_{t})\frac{2^{s}}{2n^{1/2}}\frac{2^{t}}{2n^{1/2}}$
$\displaystyle\leq
2\frac{1}{4}\underset{2^{(s+t)}\geq\frac{\sqrt{D}}{\tau_{\rm
max}}}{\sum}\frac{e(I_{s},J_{t})}{|I_{s}||J_{t}|n\alpha_{n}p_{\rm
max}^{3}}\frac{D}{\tau_{\rm max}}\frac{2^{s}|I_{s}|2^{t}|J_{t}|}{n^{2}}$
$\displaystyle=\frac{1}{2}\sqrt{D}\underset{2^{(s+t)}\geq\frac{\sqrt{D}}{\tau_{\rm
max}}}{\sum}\frac{e(I_{s},J_{t})}{\bar{\mu}(I_{s},J_{t})}\frac{\sqrt{D}}{\tau_{\rm
max}}2^{-(s+t)}\frac{2^{2s}|I_{s}|2^{2t}{J_{t}}}{n^{2}}$
$\displaystyle=\frac{1}{2}\sqrt{D}\underset{2^{(s+t)}\geq\frac{\sqrt{D}}{\tau_{\rm
max}}}{\sum}\alpha_{s}\beta_{t}\sigma_{st}.$

We will prove that ${\sum}_{2^{(s+t)}\geq\frac{\sqrt{D}}{\tau_{\rm
max}}}\alpha_{s}\beta_{t}\sigma_{st}$ is bounded by some constant. To that
end, following Feige and Ofek (2005); Lei and Rinaldo (2015); Paul et al.
(2018), we split the set of pairs
$\mathcal{C}:\\{(s,t):2^{(s+t)}\geq\frac{\sqrt{D}}{\tau_{\rm
max}},|I_{s}|\leq|J_{t}|\\}$ into six sets and we will show that the
contribution of each part is bounded. We will use the following facts
repeatedly,

$\sum_{s}\alpha_{s}\leq 4(1/2)^{-2}=16,\quad\quad\sum_{t}\beta_{t}\leq
4(1/2)^{-2}=16.$

  * $\bullet$

$\mathcal{C}_{1}:\\{(s,t)\in\mathcal{C},\sigma_{st}\leq 1\\}$.

$\sum_{s,t}\alpha_{s}\beta_{t}\sigma_{st}\mathbf{1}\\{(s,t)\in\mathcal{C}_{1}\\}\leq\sum_{s,t}\alpha_{s}\beta_{t}\leq
256.$

  * $\bullet$

$\mathcal{C}_{2}:\\{(s,t)\in\mathcal{C}\backslash\mathcal{C}_{1},\lambda_{st}\leq
ce\\}$.
Note that

$\sigma_{st}=\lambda_{st}\frac{\sqrt{D}}{\tau_{\rm
max}}2^{-(s+t)}\leq\lambda_{st}\leq ce,$

so we have

$\sum_{s,t}\alpha_{s}\beta_{t}\sigma_{st}\mathbf{1}\\{(s,t)\in\mathcal{C}_{2}\\}\leq
ce.$

  * $\bullet$

$\mathcal{C}_{3}:\\{(s,t)\in\mathcal{C}\backslash(\mathcal{C}_{1}\cup\mathcal{C}_{2}),2^{s-t}\geq\frac{\sqrt{D}}{\tau_{\rm
max}}\\}$.
Note that by Lemma 4,

$\lambda_{st}=\frac{e(I_{s},J_{t})}{\bar{\mu}_{st}}\leq
c\frac{|I_{s}|\Delta}{|I_{s}||J_{t}|\Delta/n}\leq c\frac{n}{|J_{t}|}.$

Hence

$\sigma_{st}=\lambda_{st}\frac{\sqrt{D}}{\tau_{\rm max}}2^{-(s+t)}\leq
c\frac{n}{|J_{t}|}\frac{\sqrt{D}}{\tau_{\rm max}}2^{-(s+t)}.$

And consequently

$\displaystyle\sum_{t}\beta_{t}\sigma_{st}$ $\displaystyle\leq
c\sum_{t}2^{2t}\frac{|J_{t}|}{n}\cdot\frac{n}{|J_{t}|}\frac{\sqrt{D}}{\tau_{\rm
max}}2^{-(s+t)}=c\sum_{t}\frac{\sqrt{D}}{\tau_{\rm max}2^{s-t}}\leq c,$

where the last inequality follows from the fact that the non-zero summands
over $t$ are bounded by 1 for a geometric sequence. Finally, we have

$\sum_{s,t}\alpha_{s}\beta_{t}\sigma_{st}\mathbf{1}\\{(s,t)\in\mathcal{C}_{3}\\}\leq
c\sum_{s}\alpha_{s}\leq c^{\prime}.$

  * $\bullet$

$\mathcal{C}_{4}:\\{(s,t)\in\mathcal{C}\backslash(\mathcal{C}_{1}\cup\mathcal{C}_{2}\cup\mathcal{C}_{3}),{\rm
log}\lambda_{st}>\frac{1}{4}[2t\,{\rm log}2+{\rm log}(1/\beta_{t})]\\}$.
By Lemma 5, we have

$\displaystyle\lambda_{st}{\rm log}\lambda_{st}|I_{s}||J_{t}|\frac{\Delta}{n}$
$\displaystyle=\frac{e(I_{s},J_{t})}{\bar{\mu}(I_{s},J_{t})}{\rm
log}\frac{e(I_{s},J_{t})}{\bar{\mu}(I_{s},J_{t})}\bar{\mu}(I_{s},J_{t})\leq
c_{5}\tau_{\rm max}|J_{t}|{\rm log}\frac{n}{|J_{t}|}.$

Note that $\frac{\tau_{\rm max}}{\Delta}=\frac{\tau_{\rm max}^{2}}{D}$ by the
definition of $D,\Delta$ and $p_{n}\geq\frac{{\rm log}n}{n}$. Then we have

$\displaystyle\sigma_{st}\alpha_{s}$
$\displaystyle\leq\lambda_{st}\frac{\sqrt{D}}{\tau_{\rm
max}}2^{-(s+t)}\frac{|I_{s}|2^{2s}}{n}$ $\displaystyle\leq c_{5}\frac{1}{{\rm
log}\lambda_{st}}\frac{\tau_{\rm max}^{2}}{D}{\rm
log}\big{(}\frac{n}{|J_{t}|}\big{)}\frac{\sqrt{D}}{\tau_{\rm max}}2^{s-t}$
$\displaystyle\leq c_{5}\frac{1}{{\rm
log}\lambda_{st}}\frac{2^{s-t}}{\sqrt{D}/\tau_{\rm max}}\big{[}2t{\rm
log}2+{\rm log}(\frac{1}{\beta_{t}})\big{]}$ (B.19) $\displaystyle\leq
4c_{5}\frac{2^{s-t}}{\sqrt{D}/\tau_{\rm max}}.$

And then

$\displaystyle\underset{s,t}{\sum}\alpha_{s}\beta_{t}\sigma_{st}\mathbf{1}\\{(s,t)\in\mathcal{C}_{4}\\}=\underset{t}{\sum}\beta_{t}\underset{s}{\sum}\sigma_{st}\alpha_{s}\mathbf{1}\\{(s,t)\in\mathcal{C}_{4}\\}$
$\displaystyle\leq
4c_{5}\underset{t}{\sum}\beta_{t}\underset{s}{\sum}\frac{2^{s-t}}{\sqrt{D}/\tau_{\rm
max}}\mathbf{1}\\{(s,t)\in\mathcal{C}_{4}\\}\leq
4c_{5}\underset{t}{\sum}\beta_{t}\leq c^{\prime},$

where we used the fact that $(s,t)\notin\mathcal{C}_{3}$.

  * $\bullet$

$\mathcal{C}_{5}:\\{(s,t)\in\mathcal{C}\backslash(\mathcal{C}_{1}\cup\mathcal{C}_{2}\cup\mathcal{C}_{3}\cup\mathcal{C}_{4}),2t{\rm
log}2\geq{\rm log}(1/\beta_{t})\\}$.
Since $(s,t)\notin\mathcal{C}_{4}$, we have

${\rm log}\lambda_{st}\leq\frac{1}{4}[2t{\rm log}2+{\rm
log}\beta_{t}^{-1}]\leq ct{\rm log}2,$

and thus $\lambda_{st}\leq c2^{t}$. Because $(s,t)\notin\mathcal{C}_{1}$, we
have

$1\leq\sigma_{st}=\lambda_{st}\frac{\sqrt{D}}{\tau_{\rm max}}2^{-(s+t)}\leq
c\frac{\sqrt{D}}{\tau_{\rm max}}2^{-s}.$

And because $(s,t)\notin\mathcal{C}_{2}$, we have ${\rm log}\lambda_{st}\geq
c^{\prime}$. Thus by ($\bullet$ ‣ Proof of Lemma 3) and the definition of
$\mathcal{C}_{5}$, we have

$\sigma_{st}\alpha_{s}\leq c\frac{2^{s-t}}{\sqrt{D}/\tau_{\rm max}}\cdot
4t{\rm log}2.$

Consequently,

$\displaystyle\underset{s,t}{\sum}\alpha_{s}\beta_{t}\sigma_{st}\mathbf{1}\\{(s,t)\in\mathcal{C}_{5}\\}$
$\displaystyle=\underset{t}{\sum}\beta_{t}\underset{s}{\sum}\sigma_{st}\alpha_{s}\mathbf{1}\\{(s,t)\in\mathcal{C}_{5}\\}$
$\displaystyle\leq
c\underset{t}{\sum}\beta_{t}\underset{s}{\sum}\frac{2^{s-t}}{\sqrt{D}/\tau_{\rm
max}}4t$ $\displaystyle\leq
c\underset{t}{\sum}\beta_{t}t2^{-t}\underset{s}{\sum}\frac{2^{s}}{\sqrt{D}/\tau_{\rm
max}}$ $\displaystyle\leq c\underset{t}{\sum}\beta_{t}\leq c^{\prime},$

where the last inequality follows from the fact that the non-zero summands
over $s$ are bounded by 1 for a geometric sequence.

  * $\bullet$

$\mathcal{C}_{6}:\\{(s,t)\in\mathcal{C}\backslash(\mathcal{C}_{1}\cup\mathcal{C}_{2}\cup\mathcal{C}_{3}\cup\mathcal{C}_{4}\cup\mathcal{C}_{5})\\}$.
Since $(s,t)\notin\mathcal{C}_{4}\cup\mathcal{C}_{2}\cup\mathcal{C}_{5},$ , we
have

$0<{\rm log}\lambda_{st}\leq\frac{1}{4}[2t{\rm log}2+{\rm
log}\beta_{t}^{-1}]\leq\frac{1}{2}{\rm log}\beta_{t}^{-1}\leq{\rm
log}\beta_{t}^{-1}.$

Then

$\displaystyle\underset{s,t}{\sum}\alpha_{s}\beta_{t}\sigma_{st}\mathbf{1}\\{(s,t)\in\mathcal{C}_{6}\\}$
$\displaystyle=\underset{s}{\sum}\alpha_{s}\underset{t}{\sum}\beta_{t}\lambda_{st}\frac{\sqrt{D}}{\tau_{\rm
max}}2^{-(s+t)}\mathbf{1}((s,t)\notin\mathcal{C}_{6})$
$\displaystyle\leq\underset{s}{\sum}\alpha_{s}\underset{t}{\sum}\frac{\sqrt{D}}{\tau_{\rm
max}}2^{-(s+t)}\mathbf{1}((s,t)\notin\mathcal{C}_{6})$ $\displaystyle\leq
c\underset{s}{\sum}\alpha_{s}\leq c^{\prime}.$

Putting these pieces together, and combing (A.6), for any fixed
$x,y\in\mathcal{H}$, we have

$\underset{(i,j)\in\mathcal{H}}{\sum}x_{i}y_{j}W^{M}_{ij}\leq c\sqrt{D},$

under the event $\mathcal{E}$. Taking the supremum and combining the results
of Lemma 4 and 5, we arrive the conclusion of Lemma 3. $\blacksquare$

### Proof of Lemma 4

We will make use of Warnke (2017) to bound the degree. For reference, we
reproduce the result as the following proposition.

######  Proposition 2 (Theorem 9 of Warnke (2017))

Let $(Y_{i}),i\in\mathcal{I}$ be a collection of non-negative random variables
with $\sum_{i\in\mathcal{I}}\mathbb{E}(Y_{i})\leq\mu$. Assume that $\sim$ is a
symmetric relation on $\mathcal{I}$ such that each $Y_{i}$ with
$i\in\mathcal{I}$ is independent of $\\{Y_{j}:j\in\mathcal{I},j\not\sim i\\}$.
Let $Z_{C}={\rm max}\sum_{i\in\mathcal{J}}Y_{i}$, where the maximum is taken
over all sets $\mathcal{J}\subset\mathcal{I}$ such that ${\rm
max}_{j\in\mathcal{J}}\sum_{i\in\mathcal{J},i\sim j}Y_{i}\leq C$. Then for all
$C,t>0$, we have

$\mathbb{P}(Z_{C}\geq\mu+t)\leq{\rm min}\Big{\\{}{\rm
exp}\Big{(}-\frac{t^{2}}{2C(\mu+t/3)}\Big{)},\Big{(}1+\frac{t}{2\mu}\Big{)}^{-t/2C}\Big{\\}}.$

Now we proceed to bound the degree. First, we have the observation that

$\displaystyle\mathbb{E}d_{i}^{M}=\mathbb{E}[\sum_{j}\sum_{k\neq
i,j}(\mathbf{1}(W_{ij}W_{ik}W_{jk}>0)(W_{ij}+W_{ik}+W_{jk}))]$
$\displaystyle\leq 3n^{2}\alpha_{n}p_{\rm max}^{3}\leq c\Delta,$

where we recall that $\Delta=n^{2}\alpha_{n}p_{n}^{3}.$ Then let $I_{i}$ be
the set of all triangles that includes node $i$ and denote

$M_{\theta}:=\mathbf{1}(W_{ij}\cdot W_{ik}\cdot
W_{jk}>0)(W_{ij}+W_{ik}+W_{jk})$

with the index $\theta=\\{i,j,k\\}\in I_{i}$. Recall the good event $\Gamma$
we defined in (B.9), then under the good event, every pair of nodes has at
most $C=c\tau_{\rm max}$ common neighbors, and the good event holds with
probability larger than $1-n^{-c^{\prime\prime}}$. It is clear that two
triangles belonging to the set $I_{i}$ are independent if they do not share
any edges. Let $\sim$ denote a relation such that $\theta_{1}\sim\theta_{2}$
holds if $\theta_{1}$ and $\theta_{2}$ share an edge. For any $M_{ijk}$, the
good event restricts the number of triangles in the set $I_{i}$ that are
dependent with $M_{ijk}$ to $2C$.

Define $J_{i}\subseteq I_{i}$ as

$J_{i}=\\{\theta:\underset{\theta_{1}\in J_{i}}{{\rm max}}\;|\theta_{2}\in
J;\theta_{1}\;{\mbox{and}}\;\theta_{2}\;\mbox{ share at least one edge}|\leq
2C\\}.$

Then, we immediately have

$\underset{\theta_{2}\in J_{i}}{{\rm max}}\;\underset{\theta_{1}\in
J_{i}:\theta_{1}\sim\theta_{2}}{\sum}M_{\theta}\leq 3\cdot 2C\leq c\tau_{\rm
max},\quad\sum_{\theta\in\mathcal{I}_{i}}\mathbb{E}(M_{\theta})\leq\Delta.$

Take $t=\mu=\Delta$, then the results in Proposition 2 imply,

$\displaystyle\mathbb{P}(\underset{J_{i}}{{\rm max}}\;\underset{\theta\in
J_{i}}{\sum}{M_{\theta}}\geq 2\Delta)$ $\displaystyle\leq{\rm
min}\Big{\\{}{\rm
exp}\Big{(}-\frac{\Delta^{2}}{6C(\Delta+\Delta/3)}\Big{)},\Big{(}1+\frac{\Delta}{2\Delta}\Big{)}^{-\Delta/6C}\Big{\\}}$
$\displaystyle={\rm min}\Big{\\{}{\rm exp}\Big{(}-\frac{\Delta}{c\tau_{\rm
max}}\Big{)},(\frac{3}{2})^{-\Delta/c\tau_{\rm max}}\Big{\\}}$
$\displaystyle\leq{\rm exp}(-c^{\prime}({\rm log}n)^{\beta}),$

where the last inequality follows from the assumption that
$\frac{\Delta}{\tau_{\rm max}}\geq cn\alpha_{n}p_{n}\geq c^{\prime}{({\rm
log}n)^{\beta}}$ with $0<\beta<1$. Under the good event, we have
$J^{*}_{i}=I_{i}$, and thus $\underset{J_{i}}{{\rm max}}\;\underset{\theta\in
J_{i}}{\sum}{M_{\theta}}=d_{i}^{M}$. Consequently,

$\displaystyle\mathbb{P}(d_{i}^{M}\geq 2\Delta)$ $\displaystyle\leq{\rm
exp}(-c^{\prime}({\rm log}n)^{\beta})+P(W\notin\Gamma)$ $\displaystyle\leq
n^{-c^{\prime\prime}}+{\rm exp}(-c^{\prime}({\rm log}n)^{\beta}).$

$\blacksquare$

### Proof of Lemma 5

Recall that $\bar{\mu}=|I||J|\frac{\Delta}{n}$ with
$\Delta=n^{2}\alpha_{n}p_{n}^{3}.$ So if $|J|>\frac{n}{e}$, then

$\frac{e(I,J)}{\bar{\mu}(I,J)}=\frac{\sum_{i\in I}{\rm
max}_{i}d^{M}_{i}}{|I|\Delta/e}\leq\frac{c|I|\Delta}{|I|\Delta/e}\leq
c^{\prime}e,$

which is the case $(a)$ in the lemma.

Now we prove the case $(b)$ also holds. Denote $S$ as the set of all 3-tuples
such that each tuple has one vertex in each of the sets $I$ and $J$. Then

$e(I,J)=\sum_{\theta\in S(I,J)}W^{M}_{\theta},$

and

$\mathbb{E}(\sum_{\theta\in
S(I,J)}W^{M}_{\theta})\leq|I||J|\Delta/n=\bar{\mu}(I,J).$

Since $e(I,J)$ is the sum of dependent variables, next we use the
concentration inequality in Proposition 2 to bound $e(I,J)$. Recall the good
event $\Gamma$ we defined in (A.19), under the good event, every pair of nodes
has at most $C=c\tau_{\rm max}$ common neighbors, and the good event holds
with probability larger than $1-n^{-c^{\prime\prime}}$. Define $S^{*}\subset
S$ be the set such that any $W_{\theta}^{M}$ depends on $2C$ other
$W_{\theta}^{M}$. Then under the good event, we have

$\underset{\theta_{1}\in S^{*}}{{\rm max}}\underset{\theta_{2}\in
S^{*}(I,J),\theta_{2}\sim\theta_{1}}{\sum}W^{M}_{\theta}\leq c\tau_{\rm max}.$

Let $t=(l-1)\bar{\mu}(I,J)$, and apply Proposition 2, we can obtain

$\displaystyle\mathbb{P}(e(I,J)\geq l\bar{\mu}(I,J))$ $\displaystyle\leq{\rm
exp}\big{(}-\frac{t}{4\tau_{\rm max}}{\rm
log}(1+\frac{t}{2\bar{\mu}(I,J)})\big{)}$ $\displaystyle\leq{\rm
exp}\big{(}-\frac{(\bar{\mu}(I,J)+t){\rm
log}(1+\frac{t}{\bar{\mu}(I,J)})-t}{\tau_{\rm max}}\big{)}$
$\displaystyle={\rm exp}\big{(}\frac{l\;{\rm
log}\;l\bar{\mu}(I,J)-(l-1)\bar{\mu}(I,J)}{\tau_{\rm max}}\big{)},$
$\displaystyle\leq{\rm exp}\big{(}-\frac{1}{2}\frac{l\;{\rm
log}\;l\bar{\mu}(I,J)}{\tau_{\rm max}}\big{)},$

where the second inequality follows from the fact that

$\bar{\mu}(I,J){\rm log}\;l-(l-1)\bar{\mu}(I,J)\leq(l-1)\bar{\mu}(I,J)/2{\rm
log}(1+(l-1)/2),$

for large enough $l$, and the last inequality holds when $l\geq 8$. The
remaining part of the proof is almost the same as that in Paul et al. (2018);
Feige and Ofek (2005); Lei and Rinaldo (2015), hence we omit it.
$\blacksquare$

## C. Clustering results in real data analysis

(a) Detected communities

(b) High-dimensional statistics community (c) Non-parametric statistics and
Functional analysis community

Figure 10: The clustering results for the statisticians’ citation network by
the motif-based spectral clustering. (a) is the general clustering results,
and (b) and (c) are specific communities.

(a) Detected communities

(b) High-dimensional and Bayesian statistics _mixed_ community (c) Non-
parametric statistics and Functional analysis community

Figure 11: The clustering results for the statisticians’ citation network by
the edge-based spectral clustering. (a) is the general clustering results, and
(b) and (c) are specific communities.

(a) Detected communities

(b) Biostatistics and medical statistics community (c) Bayesian statistics
community

Figure 12: The clustering results for the statisticians’ coauthor network by
the motif-based spectral clustering. (a) is the general clustering results,
and (b) and (c) are specific communities.

(a) Detected communities

(b) Biostatistics and medical statistics (at Harvard) community (c)
Biostatistics/medical statistics and non-parametric statistics _mixed_
community

Figure 13: The clustering results for the statisticians’ coauthor network by
edge-based spectral clustering. (a) is the general clustering results, and (b)
and (c) are specific communities.

(a) Detected communities

(b) 0.6486 (c) 0.7895 (d) 0.8667 (e) 1 (f) 0.7612

Figure 14: The clustering results for the ${\rm PM}_{2.5}$ data by the motif-
based spectral clustering. (a) is the general clustering results, and (b)-(f)
are specific communities where their clustering coefficients are also shown.

(a) Detected communities

(b) 1 (c) 0.6400 (d) 1 (e) 0.6279 (f) 0.75

Figure 15: The clustering results for the ${\rm PM}_{2.5}$ data by the edge-
based spectral clustering. (a) is the general clustering results, and (b)-(f)
are specific communities where their clustering coefficients are also shown.

## References

  * Abbe (2018) Emmanuel Abbe.  Community detection and stochastic block models: recent developments.  _The Journal of Machine Learning Research_ , 18(1):6446–6531, 2018.
  * Aicher et al. (2015) Christopher Aicher, Abigail Z Jacobs, and Aaron Clauset.  Learning latent block structure in weighted networks.  _Journal of Complex Networks_ , 3(2):221–248, 2015.
  * Benson et al. (2016) Austin R Benson, David F Gleich, and Jure Leskovec.  Higher-order organization of complex networks.  _Science_ , 353(6295):163–166, 2016.
  * Bhatia (1997) Rajendra Bhatia.  Graduate texts in mathematics: Matrix analysis, 1997.
  * Cammarata and Ke (2022) Louis Cammarata and Zheng Tracy Ke.  Power enhancement and phase transitions for global testing of the mixed membership stochastic block model.  _arXiv preprint arXiv:2204.11109_ , 2022.
  * Cerqueira and Levina (2023) Andressa Cerqueira and Elizaveta Levina.  A pseudo-likelihood approach to community detection in weighted networks.  _arXiv preprint arXiv:2303.05909_ , 2023.
  * Chen and Chen (2018) Yinghan Chen and Yuguo Chen.  An efficient sampling algorithm for network motif detection.  _Journal of Computational and Graphical Statistics_ , 27(3):503–515, 2018.
  * Chin et al. (2015) Peter Chin, Anup Rao, and Van Vu.  Stochastic block model and community detection in sparse graphs: A spectral algorithm with optimal rate of recovery.  In _Conference on Learning Theory_ , pages 391–423, 2015.
  * Clauset et al. (2004) Aaron Clauset, Mark EJ Newman, and Cristopher Moore.  Finding community structure in very large networks.  _Physical review E_ , 70(6):066111, 2004.
  * Cucuringu et al. (2020) Mihai Cucuringu, Huan Li, He Sun, and Luca Zanetti.  Hermitian matrices for clustering directed graphs: insights and applications.  In _International Conference on Artificial Intelligence and Statistics_ , pages 983–992. PMLR, 2020.
  * Davis and Kahan (1970) Chandler Davis and William Morton Kahan.  The rotation of eigenvectors by a perturbation. iii.  _SIAM Journal on Numerical Analysis_ , 7(1):1–46, 1970.
  * Feige and Ofek (2005) Uriel Feige and Eran Ofek.  Spectral techniques applied to sparse random graphs.  _Random Structures & Algorithms_, 27(2):251–275, 2005.
  * Foygel and Drton (2010) Rina Foygel and Mathias Drton.  Extended bayesian information criteria for gaussian graphical models.  In _Advances in neural information processing systems_ , pages 604–612, 2010.
  * Friedman et al. (2008) Jerome Friedman, Trevor Hastie, and Robert Tibshirani.  Sparse inverse covariance estimation with the graphical lasso.  _Biostatistics_ , 9(3):432–441, 2008.
  * Gao et al. (2017) Chao Gao, Zongming Ma, Anderson Y Zhang, and Harrison H Zhou.  Achieving optimal misclassification proportion in stochastic block models.  _The Journal of Machine Learning Research_ , 18(1):1980–2024, 2017.
  * Ghoshdastidar and Dukkipati (2014) Debarghya Ghoshdastidar and Ambedkar Dukkipati.  Consistency of spectral partitioning of uniform hypergraphs under planted partition model.  In _Advances in Neural Information Processing Systems_ , pages 397–405, 2014.
  * Ghoshdastidar and Dukkipati (2017) Debarghya Ghoshdastidar and Ambedkar Dukkipati.  Uniform hypergraph partitioning: Provable tensor methods and sampling techniques.  _The Journal of Machine Learning Research_ , 18(1):1638–1678, 2017.
  * Ghoshdastidar et al. (2017) Debarghya Ghoshdastidar, Ambedkar Dukkipati, et al.  Consistency of spectral hypergraph partitioning under planted partition model.  _The Annals of Statistics_ , 45(1):289–315, 2017\.
  * Goldenberg et al. (2010) Anna Goldenberg, Alice X Zheng, Stephen E Fienberg, Edoardo M Airoldi, et al.  A survey of statistical network models.  _Foundations and Trends® in Machine Learning_ , 2(2):129–233, 2010.
  * Guo et al. (2020) Xiao Guo, Yixuan Qiu, Hai Zhang, and Xiangyu Chang.  Randomized spectral co-clustering for large-scale directed networks.  _arXiv preprint arXiv:2004.12164_ , 2020.
  * Holland and Leinhardt (1977) Paul W Holland and Samuel Leinhardt.  A method for detecting structure in sociometric data.  In _Social Networks_ , pages 411–432. Elsevier, 1977.
  * Holland et al. (1983) Paul W Holland, Kathryn Blackmond Laskey, and Samuel Leinhardt.  Stochastic blockmodels: First steps.  _Social networks_ , 5(2):109–137, 1983.
  * Ji and Jin (2016) Pengsheng Ji and Jiashun Jin.  Coauthorship and citation networks for statisticians.  _The Annals of Applied Statistics_ , 10(4):1779–1812, 2016.
  * Jin et al. (2018) Jiashun Jin, Zheng Ke, and Shengming Luo.  Network global testing by counting graphlets.  In _International conference on machine learning_ , pages 2333–2341. PMLR, 2018.
  * Jin et al. (2021) Jiashun Jin, Zheng Tracy Ke, and Shengming Luo.  Optimal adaptivity of signed-polygon statistics for network testing.  _The Annals of Statistics_ , 49(6):3408–3433, 2021.
  * Karwa and Slavković (2016) Vishesh Karwa and Aleksandra Slavković.  Inference using noisy degrees: Differentially private beta-model and synthetic graphs.  _The Annals of Statistics_ , 44(1):87–112, 2016\.
  * Karwa et al. (2017) Vishesh Karwa, Pavel N Krivitsky, and Aleksandra B Slavković.  Sharing social network data: differentially private estimation of exponential family random-graph models.  _Journal of the Royal Statistical Society: Series C (Applied Statistics)_ , 66(3):481–500, 2017.
  * Kolaczyk (2009) Eric D Kolaczyk.  In _Statistical Analysis of Network Data: Methods and Models_. Springer, 2009.
  * Laenen and Sun (2020) Steinar Laenen and He Sun.  Higher-order spectral clustering of directed graphs.  _Advances in Neural Information Processing Systems_ , 33:941–951, 2020.
  * Lei and Rinaldo (2015) Jing Lei and Alessandro Rinaldo.  Consistency of spectral clustering in stochastic block models.  _The Annals of Statistics_ , 43(1):215–237, 2015\.
  * Mangan et al. (2003) Shmoolik Mangan, Alon Zaslaver, and Uri Alon.  The coherent feedforward loop serves as a sign-sensitive delay element in transcription networks.  _Journal of molecular biology_ , 334(2):197–204, 2003.
  * Manning et al. (2010) Christopher Manning, Prabhakar Raghavan, and Hinrich Schütze.  Introduction to information retrieval.  _Natural Language Engineering_ , 16(1):100–103, 2010.
  * Milo et al. (2002) Ron Milo, Shai Shen-Orr, Shalev Itzkovitz, Nadav Kashtan, Dmitri Chklovskii, and Uri Alon.  Network motifs: simple building blocks of complex networks.  _Science_ , 298(5594):824–827, 2002.
  * Newman (2018) Mark Newman.  _Networks_.  Oxford university press, 2018.
  * Newman (2006) Mark EJ Newman.  Modularity and community structure in networks.  _Proceedings of the national academy of sciences_ , 103(23):8577–8582, 2006.
  * Paul and Chen (2020) Subhadeep Paul and Yuguo Chen.  Spectral and matrix factorization methods for consistent community detection in multi-layer networks.  _Annals of Statistics_ , 48(1):230–250, 2020\.
  * Paul et al. (2018) Subhadeep Paul, Olgica Milenkovic, and Yuguo Chen.  Higher-order spectral clustering under superimposed stochastic block model.  _arXiv preprint arXiv:1812.06515_ , 2018.
  * Priebe et al. (2019) Carey E Priebe, Youngser Park, Joshua T Vogelstein, John M Conroy, Vince Lyzinski, Minh Tang, Avanti Athreya, Joshua Cape, and Eric Bridgeford.  On a two-truths phenomenon in spectral graph clustering.  _Proceedings of the National Academy of Sciences_ , 116(13):5995–6000, 2019.
  * Qin and Rohe (2013) Tai Qin and Karl Rohe.  Regularized spectral clustering under the degree-corrected stochastic blockmodel.  In _Advances in Neural Information Processing Systems_ , pages 3120–3128, 2013.
  * Rohe and Qin (2013) K Rohe and T Qin.  The blessing of transitivity in sparse and stochastic networks.  _arXiv_ , pages 1–29, 2013.
  * Rohe et al. (2011) Karl Rohe, Sourav Chatterhee, and Bin Yu.  Spectral clustering and the high-dimensional stochastic block model.  _The Annals of Statistics_ , 39(4):1878–1915, 2011.
  * Rosvall et al. (2014) Martin Rosvall, Alcides V Esquivel, Andrea Lancichinetti, Jevin D West, and Renaud Lambiotte.  Memory in network flows and its effects on spreading dynamics and community detection.  _Nature communications_ , 5:4630, 2014.
  * Serrour et al. (2011) Belkacem Serrour, Alex Arenas, and Sergio Gómez.  Detecting communities of triangles in complex networks using spectral optimization.  _Computer Communications_ , 34(5):629–634, 2011\.
  * Seshadhri et al. (2014) C Seshadhri, Ali Pinar, and Tamara G Kolda.  Wedge sampling for computing clustering coefficients and triangle counts on large graphs.  _Statistical Analysis and Data Mining: The ASA Data Science Journal_ , 7(4):294–307, 2014.
  * Tsourakakis et al. (2017) Charalampos E Tsourakakis, Jakub Pachocki, and Michael Mitzenmacher.  Scalable motif-aware graph clustering.  In _Proceedings of the 26th International Conference on World Wide Web_ , pages 1451–1460. International World Wide Web Conferences Steering Committee, 2017.
  * Underwood et al. (2020) William G Underwood, Andrew Elliott, and Mihai Cucuringu.  Motif-based spectral clustering of weighted directed networks.  _Applied Network Science_ , 5(1):1–41, 2020.
  * Von Luxburg (2007) Ulrike Von Luxburg.  A tutorial on spectral clustering.  _Statistics and computing_ , 17(4):395–416, 2007\.
  * Wang et al. (2018) Bo Wang, Armin Pourshafeie, Marinka Zitnik, Junjie Zhu, Carlos D Bustamante, Serafim Batzoglou, and Jure Leskovec.  Network enhancement as a general method to denoise weighted biological networks.  _Nature communications_ , 9(1):1–8, 2018.
  * Warnke (2016) Lutz Warnke.  On the method of typical bounded differences.  _Combinatorics, Probability and Computing_ , 25(2):269–299, 2016.
  * Warnke (2017) Lutz Warnke.  Upper tails for arithmetic progressions in random subsets.  _Israel Journal of Mathematics_ , 221(1):317–365, 2017.
  * Xu et al. (2020) Min Xu, Varun Jog, and Po-Ling Loh.  Optimal rates for community estimation in the weighted stochastic block model.  _The Annals of Statistics_ , 48(1):183–204, 2020\.
  * Yang and Leskovec (2014) Jaewon Yang and Jure Leskovec.  Overlapping communities explain core–periphery organization of networks.  _Proceedings of the IEEE_ , 102(12):1892–1902, 2014.
  * Yin et al. (2017) Hao Yin, Austin R Benson, Jure Leskovec, and David F Gleich.  Local higher-order graph clustering.  In _Proceedings of the 23rd ACM SIGKDD international conference on knowledge discovery and data mining_ , pages 555–564, 2017.
  * Yuan and Lin (2007) Ming Yuan and Yi Lin.  Model selection and estimation in the gaussian graphical model.  _Biometrika_ , 94(1):19–35, 2007.
  * Zhang et al. (2022) Hai Zhang, Xiao Guo, and Xiangyu Chang.  Randomized spectral clustering in large-scale stochastic block models.  _Journal of Computational and Graphical Statistics_ , 31(3):887–906, 2022.