##footnotetext: Maoyu Zhang and Jingfei Zhang are joint first authors.

# Fast Community Detection in Dynamic and Heterogeneous Networks

Maoyu Zhang${}^{\text{1}}$, Jingfei Zhang${}^{\text{2}}$ and Wenlin
Dai${}^{\text{1}}$
${}^{\text{1}}$ Institute of Statistics and Big Data, Renmin University of
China
${}^{\text{2}}$ Department of Management Science, University of Miami

###### Abstract

Dynamic heterogeneous networks describe the temporal evolution of interactions
among nodes and edges of different types. While there is a rich literature on
finding communities in dynamic networks, the application of these methods to
dynamic heterogeneous networks can be inappropriate, due to the involvement of
different types of nodes and edges and the need to treat them differently.In
this paper, we propose a statistical framework for detecting common
communities in dynamic and heterogeneous networks. Under this framework, we
develop a fast community detection method called DHNet that can efficiently
estimate the community label as well as the number of communities. An
attractive feature of DHNet is that it does not require the number of
communities to be known a priori, a common assumption in community detection
methods. While DHNet does not require any parametric assumptions on the
underlying network model, we show that the identified label is consistent
under a time-varying heterogeneous stochastic block model with a temporal
correlation structure and edge sparsity. We further illustrate the utility of
DHNet through simulations and an application to review data from Yelp, where
DHNet shows improvements both in terms of accuracy and interpretability over
existing solutions.

Keywords: dynamic heterogeneous network, modularity, community detection, null
model, consistency, Yelp reviews.

##  1 Introduction

One of the fundamental problems in network data analysis is community
detection that aims to divide the network into non-overlapping groups of nodes
such that nodes within the same community are densely connected and nodes from
different communities are relatively sparsely connected. Community detection
can provide valuable insights on the organization of a network and greatly
facilitate the analysis of network characteristics. As such, community
detection methods have been applied to numerous scientific fields such as
social science (Moody and White, 2003), biology (Sørlie et al., 2001) and
business (Linden et al., 2003). Over the past few decades, the problem of
community detection has been approached from methodological, algorithmic and
theoretical perspectives with substantial developments. We refer to Fortunato
(2010) and Abbe (2017) for comprehensive reviews on this topic.

While the majority of existing community detection methods are developed for a
homogeneous network or a dynamic network, networks that are dynamic and
heterogeneous are fast emerging in recent years. For example, in a dynamic
healthcare network, nodes can be patients, diseases, doctors and hospitals and
edges can be in the type of patient-disease (patient treated for disease) and
patient-doctor (patient treated by doctor) and doctor-hospital (doctor works
at hospital). These edges are expected to evolve with time as patients may
develop new diseases that are treated by different doctors at possibly
different hospitals. Figure 1 provides an illustration of a dynamic
heterogeneous Yelp review network, which is analyzed in Section 6. In this
figure, there are three types of nodes including users, businesses and
categories and three types of edges including user-user (user is friend with
user), user-business (business is reviewed by user) and business-category
(business is labeled with category). As users review different businesses over
time, this network is both heterogeneous and dynamic.

Figure 1: The dynamic heterogeneous Yelp network with two communities,
including three types of nodes (user, business, category) and edges (user is
friend with user, business is reviewed by user, business is labeled with
category).

Due to the rich information embedded in a dynamic heterogeneous network, many
methods have been developed recently for its analysis, such as network
embedding (Wang et al., 2022; Zhang et al., 2022), representation learning
(Yin et al., 2019) and link prediction (Xue et al., 2020; Jiang et al., 2021).
However, community detection in dynamic heterogeneous networks is less
studied. One relevant work is Sun et al. (2010), which provides a mixture
model-based generative model for estimating the community structure, which is
assumed to be time-varying. Other works on this topic include Sengupta and
Chen (2015) and Zhang and Cao (2017), though they only focus on a single
heterogeneous network.

In our work, we focus on detecting common communities in a dynamic
heterogeneous network, that is, the community assignment does not vary with
time but the interactions within and between communities do. Finding common
communities are useful in many applications. For example, in genetic studies
and brain connectivity studies, the common communities represent functional
groups of genes or brain regions that are coordinated in biological processes,
and identifying them is of keen scientific interests (Zhang and Cao, 2017;
Zhang et al., 2020). In the Yelp review network, it is plausible that
businesses, categories and the majority of users have a common community
structure over time, as the service offered by a business and the interests of
users (e.g., pets, parks, fine dining) are often stable over a period of time.
One notable advantage of considering a common community structure is that the
networks observed at different time points are allowed to be highly sparse if
$S$, the number of time points, increases. For example, we show in Theorem 1
that consistent community detection is achievable as long as $\lambda
S\rightarrow\infty$, where $\lambda$ is the average degree, while the single
network case requires $\lambda\rightarrow\infty$ to achieve community
detection consistency. Moreover, our approach allows the community strength to
be highly variable over time. For example, a community needs to be active for
only a very short period of time for it to be consistently identified; see
more discussions after Theorem 1.

In this paper, we propose a statistical framework for modularity-based common
community detection in the dynamic heterogeneous network, where no parametric
assumptions are made on the model underlying the observed networks. Under this
framework, we develop a fast community detection method called DHNet that can
efficiently estimate the community label as well as the number of communities.
An attractive feature of DHNet is that it does not require the number of
communities to be known a priori, a common assumption in community detection
methods. Although DHNet does not rely on parametric assumptions on the
underlying network model, we propose a new dynamic heterogeneous stochastic
block model with a temporal correlation structure and edge sparsity, and show
that DHNet can consistently estimate the community label under this model.
This provides theoretical justifications of the proposed method and also sheds
lights on how different network properties (e.g., sparsity, size, community
strength) affect its performance. The consistency property of our method when
applied to dynamic bi-partite or multi-partite networks follows as special
cases.

The remainder of the article is organized as follows. Section 2 describes a
community detection framework and proposes a modularity function for finding
common communities in a dynamic heterogeneous network. Section 3 describes a
fast community detection method called DHNet that can efficiently estimate the
community label as well as the number of communities. Section 4 shows the
consistency property of DHNet under a dynamic heterogeneous stochastic block
model. Section 5 demonstrates the efficacy of DHNet through simulation studies
and Section 6 applies the proposed method to review data from Yelp. The paper
is concluded with a short discussion section.

##  2 Community Detection with Modularity

###  2.1 Notation

We write $[m]=\\{1,\ldots,m\\}$ for an integer $m>0$. To ease notation, we
start the introduction with a single heterogeneous networks with $L$ types of
nodes. Let $V^{[l]}=(v^{[l]}_{1},\ldots,v^{[l]}_{n_{l}})$ be the set
containing the $l$-th type of nodes for $l\in[L]$, where $n_{l}$ is the number
of $l$-th type nodes. Denote the heterogeneous network as
$\mathcal{G}=(\cup_{l=1}^{L}V^{[l]},\mathcal{E}\cup\mathcal{E}^{+})$, where
set $\mathcal{E}$ contains edges between nodes of the same type and set
$\mathcal{E}^{+}$ contains edges between nodes of different types. When
$\mathcal{E}=\emptyset$, $\mathcal{G}$ forms a multi-partite network, i.e.,
edges are only established between different types of nodes. Let $G^{[l]}$
denote the homogeneous network formed within node set $V^{[l]}$ with an
$n_{l}\times n_{l}$ adjacency matrix $A^{[l]}$, and
$G^{[l_{1}l_{2}]}=(V^{[l_{1}]}\cup V^{[l_{2}]},E^{[l_{1}l_{2}]})$ denote the
bi-partite network formed between node sets $V^{[l_{1}]}$ and $V^{[l_{2}]}$
with an $n_{l_{1}}\times n_{l_{2}}$ bi-adjacency matrix $A^{[l_{1}l_{2}]}$,
$l_{1},l_{2}\in[L]$. See Figure 2 for an example of a heterogeneous network
with $L=2$.

Figure 2: An illustrative example of a heterogeneous network with two types of
nodes.

Consider a dynamic heterogeneous network
$\\{\mathcal{G}(t),t\in\mathcal{T}\\}$ with $L$ types of nodes, where
$\mathcal{G}(t)=(\bigcup_{l=1}^{L}V^{[l]},\mathcal{E}(t)\cup\mathcal{E}^{+}(t))$
is a heterogeneous network at time $t$ defined as above. The network
$\mathcal{G}(t)$ at time $t$ can be uniquely represented by its adjacency
matrix $\mathcal{A}(t)$ defined as

$\mathcal{A}(t)=\left(\begin{array}[]{ccc}A^{[11]}(t)&\ldots&A^{[1L]}(t)\\\
\vdots&\ddots&\vdots\\\ A^{[L1]}(t)&\ldots&A^{[LL]}(t)\end{array}\right),$

where $A^{[l_{1}l_{2}]}(t)\in\mathbb{R}^{n_{l_{1}}\times n_{l_{2}}}$ is
defined as in Figure 2. Define
$\mathbf{d}^{[l]}(t)=(d^{[l]}_{1}(t),\ldots,d^{[l]}_{n_{1}}(t))$, where
$d^{[l]}_{i}(t)$ is the number of links incident to $v^{[l]}_{i}$ from
$V^{[l]}$ at time $t$, and
$\mathbf{d}^{[l_{1}l_{2}]}(t)=(d^{[l_{1}l_{2}]}_{1}(t),\ldots,d^{[l_{1}l_{2}]}_{n_{1}}(t))$,
where $d^{[l_{1}l_{2}]}_{i}$ is the number of links incident to
$v^{[l_{1}]}_{i}$ from $V^{[l_{2}]}$ at time $t$. Write the number of edges in
$A^{[l_{1}l_{2}]}(t)$ as
$m^{[l_{1}l_{2}]}(t)=\sum_{i,j}A_{ij}^{l_{1}l_{2}}(t)$ for $l_{1},l_{2}\in[L]$
.

###  2.2 Modularity function

The modularity function measures the strength of division of a network into
communities, and the maximum modularity function value is a metric frequently
used for quantifying the strength of community structure within a network
(Fortunato, 2010). The function was first defined in Newman and Girvan (2004)
for a simple network $G(V,E)$ with $n$ nodes, $m$ edges, adjacency matrix
$A_{n\times n}$ and a community assignment $\bm{e}=(e_{1},\ldots,e_{n})$,
where $e_{i}\in[K]$, as

$Q(\bm{e},G)=\frac{1}{2m}\sum_{1\leq i<j\leq
n}\left[A_{ij}-\mathbb{E}(A_{ij})\right]1(e_{i}=e_{j}),$ (1)

where $1(\cdot)$ is the indicator function. In (1), the expectation
$\mathbb{E}(A_{ij})$ is calculated under a null model for random networks with
no community structure. The most common choice for the null model is the
Chung-Lu model (Newman and Girvan, 2004; Newman, 2006). In a Chung-Lu model
(Chung et al., 2006), given the expected degrees for the nodes, the
probability of having an edge between nodes $i$ and $j$ depends only on the
their expected degrees. As noted by Newman (2006), the Chung–Lu model is the
only random graph model where the probability of having an edge between nodes
$i$ and $j$ is the product of separate functions of the expected degrees of
nodes $i$ and $j$, written as $f(d_{i})f(d_{j})$, where the functions must be
the same since the adjacency matrix is symmetric. Under the Chung–Lu model, it
has been shown that every network in the null space occurs with the same
probability and there is no preference for any particular graph configuration
(Zhang and Chen, 2017), which makes the model a desirable choice as the null.
It is seen that the modularity function in (1) measures the difference between
the observed number of intra-community edges and the expected number of intra-
community edges under the null with no community structure. Correspondingly,
the community label of a network is identified by maximizing the modularity
function with respect to $\bm{e}$.

To define the modularity function in a dynamic heterogeneous network, we first
describe the corresponding null model that characterizes a dynamic
heterogeneous network with no community structure. Consider the heterogeneous
network at time $t$,
$\mathcal{G}(t)=\left(\cup_{l=1}^{L}V^{[l]},\mathcal{E}(t)\cup\mathcal{E}^{+}(t)\right)$
with degree sequence
$\bm{D}(t)=\\{\mathbf{d}^{[l_{1}l_{2}]}(t),l_{1},l_{2}\in[L]\\}$. We define a
heterogeneous Chung-Lu model as the null. Specifically, under the null, we
assume that a heterogeneous network at time $t$ is generated with

$A_{ij}^{[l_{1}l_{2}]}(t)\sim\text{Bernoulli}\left(\frac{d_{i}^{\left[l_{1}l_{2}\right]}(t)d_{j}^{\left[l_{2}l_{1}\right]}(t)}{m^{\left[l_{1}l_{2}\right]}(t)}\right),\quad
l_{1},l_{2}\in[L],$ (2)

where all edges in $\mathcal{G}(t)$ are independent. Under (2), it is easy to
show that the expected degree sequence under the null is the same as the
observed degree sequence $\bm{D}(t)$. Following the same argument as in Zhang
and Chen (2017), it can be shown that under (2), every heterogeneous network
in the null space occurs with the same probability.

Next, we move to define the modularity matrix. At time $t\in\mathcal{T}$ and
given $\mathcal{A}(t)$, we write the
$(n_{1}+\dots+n_{L})\times(n_{1}+\dots+n_{L})$ modularity matrix
$\mathcal{M}(t)$ as

$\mathcal{M}(t)=\left(\begin{array}[]{ccc}\frac{M^{[11]}(t)}{m^{[11]}(t)}&\ldots&\frac{M^{[1L]}(t)}{m^{[1L]}(t)}\\\
\vdots&\ddots&\vdots\\\
\frac{M^{[L1]}(t)}{m^{[L1]}(t)}&\ldots&\frac{M^{[LL]}(t)}{m^{[LL]}(t)}\end{array}\right),$

where
$M^{[l_{1}l_{2}]}(t)=A^{[l_{1}l_{2}]}(t)-\mathbb{E}\left(A^{[l_{1}l_{2}]}(t)\right)$.
The modularity matrix $\mathcal{M}(t)$ measures the distance between the
observed network and the expected network under the null model at time $t$.
Given the dynamic heterogeneous networks
$\\{\mathcal{A}(t),\,t\in\mathcal{T}\\}$, the integrated modularity matrix
$\mathcal{M}$ is defined as

$\mathcal{M}=\left(\begin{array}[]{ccc}\mathcal{M}^{[11]}&\ldots&\mathcal{M}^{[1L]}\\\
\vdots&\ddots&\vdots\\\
\mathcal{M}^{[L1]}&\ldots&\mathcal{M}^{[LL]}\end{array}\right),\,\,\text{where}\,\,\mathcal{M}^{[l_{1}l_{2}]}=\frac{\int_{t\in\mathcal{T}}M^{[l_{1}l_{2}]}(t)}{\bar{m}^{[l_{1}l_{2}]}},$

$\bar{m}^{[l_{1}l_{2}]}=\int_{t\in\mathcal{T}}m^{[l_{1}l_{2}]}(t)$. The
integrated modularity matrix $\mathcal{M}$ measures the distance between the
observed network and the expected network under the null model over all
$t\in\mathcal{T}$.

We are now ready to define the modularity function. Write the community
assignment label as
$\bm{e}=\left(\mathbf{e}^{[1]},\dots,\mathbf{e}^{[L]}\right)$ with
$\mathbf{e}^{[l]}=\left(e_{1}^{[l]},\ldots,e_{n_{l}}^{[l]}\right),l\in[L]$,
the modularity function of the dynamic heterogeneous network is defined as

$\displaystyle Q(\bm{e},\\{\mathcal{G}(t)\\}_{t\in\mathcal{T}})$
$\displaystyle=\frac{1}{L^{2}}\sum_{1\leq l_{1},l_{2}\leq
L}\sum_{i,j}\mathcal{M}^{[l_{1}l_{2}]}_{ij}1(e^{[l_{1}]}_{i}=e^{[l_{2}]}_{j}).$
(3)

From the above definitions, it can be shown that
$Q(\bm{e},\\{\mathcal{G}(t)\\}_{t\in\mathcal{T}})\in[-1,1]$. This modularity
function measures the overall difference between the observed number of intra-
community edges and the expected number of intra-community edges under the
null model. When $Q(\bm{e},\\{\mathcal{G}(t)\\}_{t\in\mathcal{T}})$ approaches
1, the observed number of intra-community edges is greater than the expected
values, which indicates a strong community structure. In contrast, when
$Q(\bm{e},\\{\mathcal{G}(t)\\}_{t\in\mathcal{T}})$ approaches 0, the observed
number of intra-community edges is close to the expected values under the
null, which indicates no or weak community structure.

In practice, the networks are often only observed on a number of time points
$\mathcal{T}=\left\\{t_{1},t_{2},\ldots,t_{S}\right\\}$, where $S$ is the
total number of observations or snapshots. In this case, we can define

$\mathcal{M}=\left(\begin{array}[]{ccc}\sum_{s=1}^{S}M^{[11]}(t_{s})/\bar{m}^{[1]}&\ldots&\sum_{s=1}^{S}M^{[1L]}(t_{s})/\bar{m}^{[1L]}\\\
\vdots&\ddots&\vdots\\\
\sum_{s=1}^{S}M^{[L1]}(t_{s})/\bar{m}^{[L1]}&\ldots&\sum_{s=1}^{S}M^{[LL]}(t_{s})/\bar{m}^{[L]}\end{array}\right),$
(4)

where $\bar{m}^{[l_{1}l_{2}]}=\sum_{s=1}^{S}m^{[l_{1}l_{2}]}(t_{s})$,
$l_{1},l_{2}\in[L]$, and write the modularity function as

$Q(\bm{e},\\{\mathcal{G}(t_{s})\\}_{s\in[S]})=\frac{\sum_{s=1}^{S}m^{[l_{1}l_{2}]}\left(t_{s}\right)Q^{[l_{1}l_{2}]}\left(\bm{e},\mathcal{G}\left(t_{s}\right)\right)}{\sum_{s=1}^{S}m^{[l_{1}l_{2}]}\left(t_{s}\right)},$

where
$Q^{[l_{1}l_{2}]}\left(\bm{e},\mathcal{G}\left(t_{s}\right)\right)=\frac{1}{m^{[l_{1}l_{2}]}(t_{s})L^{2}}\sum_{l_{1},l_{2}=1}^{L}\sum_{i,j}M^{[l_{1}l_{2}]}_{ij}(t_{s})1\left(e^{[l_{1}]}_{i}=e^{[l_{2}]}_{j}\right)$.
The above modularity function can be considered as an averaged version of the
modularity in each graph $\mathcal{G}(t_{s})$, $s\in[S]$.

##  3 Modularity maximization

We aim to find the community assignment that maximizes the modularity function
(3), that is,

$\hat{\bm{c}}=\arg\max\limits_{\begin{subarray}{c}\bm{e}=(\mathbf{e}^{[1]},\ldots,\mathbf{e}^{[L]}),\\\
e_{i}^{[l]}\in\\{1,\ldots,K\\}\end{subarray}}{Q}(\bm{e},\\{\mathcal{G}(t)\\}_{t\in\mathcal{T}}).$
(5)

Finding the exact maximizer of (3) is challenging due to the combinatorial
nature of the problem and the fact that the number of communities $K$ is
generally unknown. Brandes et al. (2008) showed that finding the partition
that maximizes the modularity function for a simple graph is NP-hard. There
are a number of existing heuristic algorithmic solutions to maximizing the
modularity function, some of which are fast and hence feasible for very large
networks (Clauset et al., 2004; Wakita and Tsurumi, 2007; Blondel et al.,
2008), while some others could be more precise though restricted to graphs of
moderate sizes (Guimera et al., 2004; Massen and Doye, 2005).

In our approach, we adopt a fast Louvain-type maximization method. The Louvain
method was first proposed by Blondel et al. (2008) for modularity maximization
in simple graphs. In the Louvain method, small communities are first
identified by optimizing the modularity function locally on all nodes. Then
each small community is grouped into one “meta” node and the first step is
repeated. The Louvain method is fast to compute and enjoys a good empirical
performance. It has been successful applied to network analyses from various
scientific fields, permitting up to 100 million nodes and billions of edges.
Notably, the modularity maximum found by the Louvain method often compares
favorably with those found by alternative methods such as Clauset et al.
(2004) and Wakita and Tsurumi (2007); see Fortunato (2010).

Motivated by the Louvian algorithm, we propose a dynamic heterogeneous network
modularity maximization algorithm, referred to DHNet. To do so, we first
define a unit, which is set of nodes with at most one from each node type. For
example, a unit may contain one node of any type or $L$ nodes of different
types. A unit serves as the building block of a community in a heterogeneous
network. Next, given a heterogeneous $n\times n$ modularity matrix
$\mathcal{M}$ as in (4), we define a modularity network, which is a network of
$n$ nodes and the edge between nodes $(i,j)$ is $\mathcal{M}_{i,j}$. From (3)
and (5), it is easy to see that our optimization task is to find a partition
of the modularity network such that the within-community sum of edges from
$\mathcal{M}$ is maximized.

The algorithm DHNet starts with assigning each node to its own unit and then
each unit to its own community, leading to $n$ communities at the start of the
algorithm with each community containing only one node (or unit). The
optimization procedure is then carried out in two phases that are repeated
iteratively. In the first phase and for each unit $i$, DHNet removes this unit
from its current community and assigns it to its neighboring community
(communities to which unit $i$ is linked to), such that it leads to the
largest increase of the modularity in (3). If no move increases the
modularity, then unit $i$ remains in its current community. In the second
phase, the algorithm merges nodes of the same type in each community, such
that each community contains at most one node from each node type, and builds
a new modularity network. In the new modularity network, the units are
communities from the first phase and the edge between two nodes are given by
summing the edge weights connecting two corresponding sets of nodes from the
first phase. These steps are repeated iteratively until the modularity value
no longer increases. The algorithm can be summarized as Algorithm 1.

Algorithm 1 Dynamic Heterogeneous Network Modularity Maximization (DHNet)

Input: Dynamic heterogeneous networks $A_{1}$,…, $A_{S}$.

Step 1: Calculate the modularity matrix using (4).

Step 2: Assign each node to its own unit and assign each unit to its own
community.

Step 3: Repeat Steps 3.1-3.4 until the modularity value no longer increase.

Step 3.1: For each unit, place it into the neighboring community that leads to
the

largest modularity increase in (3). If no such move is possible, then this
unit stays

in its present community.

Step 3.2: Repeated apply Step 3.1 to all units until none can be moved.

Step 3.3: If the modularity is higher than that from the previous iteration,
merge

nodes of the same type in each community such that each community is regarded

as a unit and go to Step 3.4. If not, exit with the assignment from the
previous

iteration.

Step 3.4: Calculate the modularity matrix of the merged network.

Output: Community assignment and the corresponding modularity value.

Figure 3: A simple illustration of DHNet. Nodes of the same type are marked
using the same shape and nodes of the same color are in the same community.

Figure 3 shows an example of applying DHNet to a dynamic heterogeneous network
with two types of nodes. First, Step 1 calculates the modularity matrix
$\mathcal{M}$ from $S$ heterogeneous networks $A_{1}$,…, $A_{S}$ as in (4) and
Step 2 assigns each node to its own unit and assign each unit to its own
community. Then Step 3 is implemented once and the algorithm reaches
convergence. Specifically, in Figure 3 and after Steps 3.1-3.2, nodes
$\\{1,2,7\\}$, $\\{3,4,8\\}$ and $\\{5,6,9\\}$ are placed into three
communities, colored black, red and blue, respectively. Step 3.3 then produces
a new modularity network with 6 nodes and 3 units. The first unit has nodes
$\\{1,2\\}$ and 7, where nodes $1$ and $2$ are merged as they are of the same
type. Similarly, the second unit has $\\{3,4\\}$ and 8 and the third unit has
$\\{5,6\\}$ and 9. After this step, merging any of the three units cannot
further increase the modularity and thus DHNet returns three communities with
nodes $\\{1,2,7\\}$, $\\{3,4,8\\}$ and $\\{5,6,9\\}$, respectively.

Remark 1 (initialization). In Step 3.1, if there are multiple communities that
lead to the same maximum modularity increase, DHNet randomly selects a
community to assign the unit to. Hence, the result of DHNet may differ each
time the algorithm is implemented. Moreover, the result of the algorithm may
differ depending the node ordering in Step 2. That is, a node ordering of
$\\{1,2,3\\}$ or $\\{3,1,2\\}$ may give different results. We recommend
applying the Louvain method $\kappa$ times with random node orderings and
using the assignment with the largest modularity function value as the final
output. In our simulation studies and real data analysis, we set $\kappa=100$
and notice that the output from DHNet is not sensitive to node orderings.
Generally, it is recommended that $\kappa$ should increase with the size of
the network.

Figure 4: The computation time of DHNet with a varying network size $n$ and
number of communities $K$. In the left panel, we set the number of communities
$K=10$, and in the right panel, we set the network size $n=10,000$.

Remark 2 (time complexity). In DHNet, computing whether and where to move each
unit based on modularity changes is of time complexity $O(1)$. Hence, the
computation time at each iteration is roughly linear in the number of units,
which is less than or equal to the total number of nodes. Figure 4 provides
the computation time of DHNet with varying network size $n$ and number of
communities $K$. We set $S=20$ and generate these networks from DHSBM where
the inter- and intra-community connecting probabilities are 0.1 and 0.15,
respectively, in the homogeneous networks, and the inter- and intra-community
connecting probabilities are, respectively, 0.05 and 0.1 in the multi-partite
networks. All experiments are ran on an Intel(R) Xeon(R) with
$3.10\mathrm{GHz}$ and $192\mathrm{~{}GB}$ memory processor.

##  4 Consistency

In this section, we investigate the theoretical properties of DHNet for
finding common communities in a dynamic heterogeneous network. To do so, we
first propose a discrete-time heterogeneous stochastic block model with a
temporal correlation structure.

Dynamic Heterogeneous Stochastic Block Model (DHSBM)

  1. 1.

Dynamic heterogeneous network $\\{\mathcal{G}(t_{s}),\,s\in[S]\\}$ with $L$
node types has a latent community label $\bm{c}=$
$\left(\bm{c}^{[1]},\ldots,\bm{c}^{[L]}\right)$, where
$\bm{c}^{[l]}=(c_{1}^{[l]},\ldots,c_{n_{l}}^{[l]})$ and
$c_{i}^{[l]}\in\\{1,\dots,K\\}$ denotes the community that node $i$ of
type-$[l]$ belongs to, $l\in[L]$.

  2. 2.

The label $\bm{c}^{[l]}$ follows a multinomial distribution with $n_{l}$
trials and probability
$\bm{\pi}^{[l]}=\left(\pi_{1}^{[l]},\ldots,\pi_{K}^{[l]}\right)$, $l\in[L]$.

  3. 3.

Define the time-varying probability matrix $\bm{\Theta}(t_{s})$ such that

$\bm{\Theta}(t_{s})=\left(\begin{array}[]{ccc}\Theta^{[11]}(t_{s})&\ldots&\Theta^{[1L]}(t_{s})\\\
\vdots&\ddots&\vdots\\\
\Theta^{[L1]}(t_{s})&\ldots&\Theta^{[LL]}(t_{s})\end{array}\right),\,\,\text{where}\,\,\Theta^{[l_{1}l_{2}]}(t_{s})=\left(\begin{array}[]{ccc}\theta_{11}^{[l_{1}l_{2}]}(t_{s})&\cdots&\theta_{1K}^{[l_{1}l_{2}]}(t_{s})\\\
\vdots&\ddots&\vdots\\\
\theta_{K1}^{[l_{1}l_{2}]}(t_{s})&\cdots&\theta_{KK}^{[l_{1}l_{2}]}(t_{s})\end{array}\right)$

and $\theta_{k_{1}k_{2}}^{[l_{1}l_{2}]}(t_{s})$ is a function of $t_{s}$,
$l_{1},l_{2}\in[L]$, and $k_{1},k_{2}\in[K]$.

  4. 4.

Given $\bm{c}$, we treat $A_{ij}^{[l_{1}l_{2}]}\left(t_{s}\right)$’s as
independent Bernoulli random variables satisfying

$A_{ij}^{[l_{1}l_{2}]}\left(t_{s}\right)=uA^{[l_{1}l_{2}]}_{ij}\left(t_{s-1}\right)+(1-u)v^{[l_{1}l_{2}]},$

where $u\overset{iid}{\sim}\operatorname{Bernoulli}(\alpha)$, and given
$c^{[l_{1}]}_{i}=k_{1}$ and $c^{[l_{2}]}_{j}=k_{2}$,

$v^{[l_{1}l_{2}]}\overset{iid}{\sim}\text{ Bernoulli
}\left(\frac{\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)-\alpha\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)}{1-\alpha}\right),\quad
l_{1},l_{2}\in[L].$

In Assumption 4, it is possible to let
$u\overset{iid}{\sim}\operatorname{Bernoulli}(\alpha^{[l_{1}l_{2}]})$, though
we assume $\alpha^{[l_{1}l_{2}]}=\alpha$ to simplify notation. In our
theoretical analysis, we allow $\alpha^{[l_{1}l_{2}]}$ to vary with
$l_{1},l_{2}\in[L]$. Next, we require that $0\leq\alpha<1$,
$\alpha\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)\leq\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)$,
and
$\alpha\left(1-\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)\right)\leq
1-\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)$, so that the above
Bernoulli distribution is valid with the probability parameter in $[0,1]$.
Based on Assumption 4, some algebra shows that

$\mathbb{P}\left(A^{[l_{1}l_{2}]}_{ij}(t_{s})=1\right)=\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}(t_{s}),$

which shows that the marginal distribution of $A_{ij}^{[l_{1}l_{2}]}(t_{s})$
is
$\operatorname{Bernoulli}\left(\theta_{k_{1}k_{2}}^{[l_{1}l_{2}]}(t_{s})\right)$.
Hence, for a fixed $t_{s}$, $A^{[l_{1}l_{2}]}(t_{s})$ follows a stochastic
block model with a probability matrix $\Theta^{[l_{1}l_{2}]}(t_{s})$.
Additionally, under our DHSBM model, we have

$\operatorname{corr}\left(A^{[l_{1}l_{2}]}_{ij}\left(t_{s}\right),A^{[l_{1}l_{2}]}_{ij}\left(t_{s-1}\right)\right)=\alpha^{[l_{1}l_{2}]}\sqrt{\frac{\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)\left(1-\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)\right)}{\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)\left(1-\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)\right)}}.$

And for the special case $\alpha^{[l_{1}l_{2}]}=0$, $A_{ij}(t_{s}),s=1,\dots
S$, are independent. If $\Theta(t_{s})$ is constant over time, then
$\operatorname{corr}\left(A^{[l_{1}l_{2}]}_{ij}\left(t_{s}\right),A^{[l_{1}l_{2}]}_{ij}\left(t_{s-k}\right)\right)=(\alpha^{[l_{1}l_{2}]})^{k}$
for $k=1,2,\ldots$.

Next, we show the consistency property of the estimated assignment vector
$\bm{\hat{c}}$ under the DHSBM model when the network size $n$ and the number
of time points increases in that $nS\rightarrow\infty$. This regime is more
general and includes the results from Zhang and Cao (2017) and Zhang and Chen
(2020) as special cases. We say a label
$\bm{e}=\left(\bm{e}^{[1]},\ldots,\bm{e}^{[L]}\right)$ is consistent if it
satisfies

$\forall\epsilon>0,\quad
P\left[\frac{1}{n}\sum_{l=1}^{L}\sum_{i=1}^{n_{l}}I\left(e_{i}^{[l]}\neq
c_{i}^{[l]}\right)<\epsilon\right]\rightarrow 1\text{ as
}nS\rightarrow\infty,$

which stipulates that the misclassification ratio tends to zero. Here
$\hat{c}_{i}^{[l]}=c_{i}^{[l]}$ means that they belong to the same equivalent
class of label permutations. To allow sparsity, we reparameterize
$\bm{\Theta}(t_{s})$ as
$\tilde{\bm{\Theta}}(t_{s})=\rho_{n,S}\bm{\Theta}(t_{s})$, where
$\bm{\Theta}(t_{s})$ is fixed as $nS\rightarrow\infty$. This
reparameterization allows us to separate $\rho_{n,S}$, the sparsity parameter,
from the structure of the network.

###### Theorem 1

Consider a dynamic heterogeneous network
$\mathcal{G}\left(\bigcup_{i=1}^{L}V^{[i]},\mathcal{E}(t_{s})\cup\mathcal{E}^{+}(t_{s})\right)$
from the DHSBM with $\bm{c}$, $\bm{\pi}^{[l]}$’s, $\alpha$ and
$\bm{\Theta}(t_{s})$’s, and further assume that the community sizes are
balanced, i.e., $\min_{l}n_{l}/n$ is bounded away from zero. Define a $K\times
K$ matrix

$T_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})=\frac{\pi_{a}^{\left[l_{1}\right]}\pi_{b}^{\left[l_{2}\right]}\theta_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})}{\sum_{ab}\pi_{a}^{\left[l_{1}\right]}\pi_{b}^{\left[l_{2}\right]}\theta_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})}.$

Let
$W_{[ab]}^{\left[l_{1}l_{2}\right]}(t_{s})=T_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})-T_{a.}^{\left[l_{1}l_{2}\right]}(t_{s})T_{b.}^{\left[l_{1}l_{2}\right]}(t_{s})$
with
$T_{a.}^{\left[l_{1}l_{2}\right]}(t_{s})=\sum_{q=1}^{K}T_{aq}^{\left[l_{1}l_{2}\right]}(t_{s})$.
If the following assumptions hold

$\sum_{s=1}^{S}\sum_{l_{1},l_{2}}^{L}W_{aa}^{\left[l_{1}l_{2}\right]}(t_{s})>0\text{
and
}\sum_{s=1}^{S}\sum_{l_{1},l_{2}}^{L}W_{ab}^{\left[l_{l}l_{2}\right]}(t_{s})<0\quad\text{
for all }\quad a\neq b\in[K]$ (6)

and $nS\rho_{n,S}\rightarrow\infty$, then we have

$\forall\epsilon>0,\quad
P\left[\frac{1}{n}\sum_{l=1}^{L}\sum_{i=1}^{n_{l}}I\left(\hat{c}_{i}^{[l]}\neq
c_{i}^{[l]}\right)<\epsilon\right]\rightarrow 1\text{ as
}nS\rightarrow\infty,$

where $\hat{\bm{c}}$ is the maximizer of (3).

It is seen that the network is allowed to be highly sparse at each time point
$s$, e.g., the probability of forming an edge can be
$O\left(\frac{\log(nS)}{nS}\right)$. Denoting the average degree as
$\lambda=n\rho_{n,S}$, it is seen that consistency is achievable when $\lambda
S\rightarrow\infty$, while the single network case requires
$\lambda\rightarrow\infty$ to achieve community detection consistency (Zhang
and Chen, 2020). When $L=1$, the above result reduces to that in Zhang and Cao
(2017) and when $S=1$, the above result reduces to that in Zhang and Chen
(2020). We note that Zhang and Cao (2017) only considered the case where the
network size $n$ is fixed and their results require $\rho_{n,S}=O(1)$. In
comparison, our result in Theorem 1 allows $n$ and/or $S$ to diverge and only
requires $nS\rho_{n,S}\rightarrow\infty$ as $nS\rightarrow\infty$.

The condition in (6) requires that edges are on average more likely to be
established within communities than they are between communities, though
communities may not exist for all types of edges or at all time points. For
example, in the simulation setting in Section 5.2, the edges within type-[1]
nodes and/or type-[2] nodes have no community structure, while the edges
linking type-[1] and type-[2] nodes do at some time points. This type of
assortative condition, requiring more edges within communities than between
communities, is often required for algorithm-based community detection such as
modularity maximization. For the special case of $L=1$, $K=2$, and
$\bm{\Theta}(t_{s})$ is time homogeneous, the condition (6) can be simplified
as

$\theta_{11}^{[11]}\theta_{22}^{[11]}>\left(\theta_{12}^{[11]}\right)^{2}.$

When $L=2$, $K=2$ and $\bm{\Theta}(t_{s})$ is time-varying, the condition (6)
is satisfied if

$\sum_{s=1}^{S}\left(\theta_{11}^{[11]}(t_{s})+\theta_{11}^{[22]}(t_{s})+\theta_{11}^{[12]}(t_{s})+\theta_{11}^{[21]}(t_{s})\right)>\sum_{s=1}^{S}\left(\theta_{12}^{[11]}(t_{s})+\theta_{12}^{[22]}(t_{s})+\theta_{12}^{[12]}(t_{s})+\theta_{12}^{[21]}(t_{s})\right),$
$\sum_{s=1}^{S}\left(\theta_{22}^{[11]}(t_{s})+\theta_{22}^{[22]}(t_{s})+\theta_{22}^{[12]}(t_{s})+\theta_{22}^{[21]}(t_{s})\right)>\sum_{s=1}^{S}\left((\theta_{12}^{[11]}(t_{s})+\theta_{12}^{[22]}(t_{s})+\theta_{12}^{[12]}(t_{s})+\theta_{12}^{[21]}(t_{s})\right),$

which indicate that edges are more likely to form within communities than
between communities.

##  5 Simulation

In this section, we evaluate the clustering accuracy of DHNet and compare it
with several alternative solutions including:

Method 1: treat the dynamic heterogeneous network as a dynamic homogeneous
network without distinguishing the different node and edge types and apply a
dynamic network community detection method (Zhang and Cao, 2017).

Method 2: apply a heterogeneous community detection method (Zhang and Chen,
2020) to an aggregated matrix
$\mathcal{\bar{A}}=\left(\begin{array}[]{ccc}\bar{A}^{[11]}&\ldots&\bar{A}^{[1L]}\\\
\vdots&\ddots&\vdots\\\
\bar{A}^{[L1]}&\ldots&\bar{A}^{[LL]}\end{array}\right),$ where
${\bar{A}_{ij}^{[l_{1}l_{2}]}}=\max_{t}A_{ij}^{[l_{1}l_{2}]}(t)$, that is,
detect community based on a static summary heterogeneous graph.

Method 3: infer the community label from $\mathcal{G}\left(t_{s}\right)$ for a
randomly selected time point $t_{s}$ in $\\{t_{1},\ldots,t_{S}\\}$. That is,
community detection based on a single snapshot of the dynamic heterogeneous
network, which is the same as Zhang and Chen (2020).

Method 4: decompose the dynamic heterogeneous network with $L$ different types
of nodes into $L$ dynamic homogeneous networks and apply a dynamic network
community detection method (Zhang and Cao, 2017) to each separately, i.e.,
discard information from the edges linking different types of nodes.

We generate networks from the DHSBM proposed in Section 4 with $L$ types of
nodes, $K$ communities and $S$ equal-spaced observations within the time
interval $[0,1]$. We consider three different settings in our experiments
including a time-homogeneous DHSBM with independently sampled networks in
Section 5.1, a DHSBM with independently sampled networks in Section 5.2 and a
DHSBM with temporally correlated networks in Section 5.3. In each setting, we
consider dense and sparse networks. We set $L=2$, $K=3$, $n_{1}=300$,
$n_{2}=150$ and $\bm{\pi}^{[1]}=\bm{\pi}^{[2]}=(1/3,1/3,1/3)$. To evaluate the
clustering accuracy, we adopt the normalized mutual information (NMI) (Danon
et al., 2005), a commonly used metric in community detection experiments to
quantifies the difference between two clustering labels.

###  5.1 Simulation setting 1

We consider networks independently sampled from a DHSBM with a time-
homogeneous probability matrix defined as

$\bm{\Theta}(t)=\left(\begin{array}[]{ccc|ccc}\theta_{1}+r_{1}&\theta_{1}&\theta_{1}&\theta_{3}+r_{3}&\theta_{3}&\theta_{3}\\\
\theta_{1}&\theta_{1}+r_{1}&\theta_{1}&\theta_{3}&\theta_{3}+r_{3}&\theta_{3}\\\
\theta_{1}&\theta_{1}&\theta_{1}+r_{1}&\theta_{3}&\theta_{3}&\theta_{3}+r_{3}\\\
\hline\cr\theta_{3}+r_{3}&\theta_{3}&\theta_{3}&\theta_{2}+r_{2}&\theta_{2}&\theta_{2}\\\
\theta_{3}&\theta_{3}+r_{3}&\theta_{3}&\theta_{2}&\theta_{2}+r_{2}&\theta_{2}\\\
\theta_{3}&\theta_{3}&\theta_{3}+r_{3}&\theta_{2}&\theta_{2}&\theta_{2}+r_{2}\\\
\end{array}\right).$

In the type-$[1]$ (type-$[2]$) homogeneous network $G^{[1]}$ ($G^{[2]}$), the
parameter $\theta_{1}\left(\theta_{2}\right)$ represents the inter-community
connecting probability and $\theta_{1}+r_{1}\left(\theta_{2}+r_{2}\right)$
represents the intra-community connecting probability. In the type-$[12]$ bi-
partite network, $\theta_{3}$ describes the inter-community connecting
probability and $\theta_{3}+r_{3}$ describes the intra-community connecting
probability. The strength of the community structure is regulated by
$r_{1},r_{2}$ and $r_{3}$. We consider both dense and sparse networks in this
setting with scenarios 1 and 2 on dense and sparse networks, respectively.
Specifically, we consider

Scenario 1: $\theta_{1}=0.5$, $\theta_{2}=0.6$, $\theta_{3}=0.3$, $r_{1}=0$,
$r_{2}=0$,

Scenario 2: $\theta_{1}=0.1$, $\theta_{2}=0.2$, $\theta_{3}=0.05$, $r_{1}=0$,
$r_{2}=0$.

In Scenarios 1 and 2, neither $G^{[1]}$ or $G^{[2]}$ has a community
structures. We have also considered the case where $G^{[1]}$ has a weak
community structure while $G^{[2]}$ has no community structure. The results
are similar to those from Scenarios 1 and 2 and delayed to the supplement. We
set $S=20$ and vary $r_{3}$, i.e., the strength of the community structure in
$G^{[12]}$, from $0.05$ to $0.15$. Figure S1 summarizes the community
detection results averaged over 100 data replicates for Scenarios 1-2,
respectively.

(a) Type 1 nodes, Scenario 1 (b) Type 1 nodes, Scenario 2 (c) Type 2 nodes,
Scenario 1 (d) Type 2 nodes, Scenario 2

Figure 5: Average NMIs against the value of $r_{3}$ for different methods in
Setting 1 under Scenarios 1-2.

For dense networks in Scenario 1, it is seen from the left panel in Figure S1
that DHNet outperforms the other methods on all values of $r_{3}$. The NMIs
from Methods 1-3 are below 0.25 for both types of nodes. For Method 1, the
clustering output places nodes of the same type in the same community, leading
to an NMI close to zero. For Method 2, the aggregated network becomes very
dense and the number of inter-community edges are very similar to that of the
intra-community edges for each edge type, leading an NMI close to zero. Method
3 detects community based on a random snapshot of network, which contains
relatively weak structural information, leading to a lower NMI. In addition,
Method 4 ignores the edges linking different types of nodes and hence perform
well only when a strong community structure exists among the investigated type
of nodes.

For sparse networks in Scenario 2, it is seen from the right panel in Figure
S1 that the performance of Methods 1 and 3 increases notably with $r_{3}$, as
the community structure strength (i.e., the difference between the inter- and
intra- community connecting probability) is high in this scenario. Due to this
reason, Method 2 also performs better in the sparse case as the community
structure signal is strong in the aggregated network, with many more inter-
community edges than intra-community edges. Our method still outperforms most
of the other methods when the signal is weak, e.g., $r_{3}\leq 0.1$.

###  5.2 Simulation setting 2

We consider networks independently sampled from a DHSBM with a time-varying
probability matrix defined as

$\bm{\Theta}(t)=\left(\begin{array}[]{ccc|ccc}\theta_{1}+r_{1}&\theta_{1}&\theta_{1}&\theta_{3}+r_{31}(t)&\theta_{3}&\theta_{3}\\\
\theta_{1}&\theta_{1}+r_{1}&\theta_{1}&\theta_{3}&\theta_{3}+r_{32}(t)&\theta_{3}\\\
\theta_{1}&\theta_{1}&\theta_{1}+r_{1}&\theta_{3}&\theta_{3}&\theta_{3}+r_{33}(t)\\\
\hline\cr\theta_{3}+r_{31}(t)&\theta_{3}&\theta_{3}&\theta_{2}+r_{2}&\theta_{2}&\theta_{2}\\\
\theta_{3}&\theta_{3}+r_{32}(t)&\theta_{3}&\theta_{2}&\theta_{2}+r_{2}&\theta_{2}\\\
\theta_{3}&\theta_{3}&\theta_{3}+r_{33}(t)&\theta_{2}&\theta_{2}&\theta_{2}+r_{2}\\\
\end{array}\right).$

We set $\theta_{1}$, $\theta_{2}$, $\theta_{3}$, $r_{1}$ and $r_{2}$ the same
as those in the two scenarios in Simulation 1 and $r_{31}(t)$, $r_{32}(t)$ and
$r_{33}(t)$ as plotted in Figure 6. In this setting, at time $t=0$, community
1 in $G^{[12]}$ is active while communities 2-3 are inactive; at time $t=0.5$,
community 1 in $G^{[12]}$ becomes inactive while communities 2-3 are active;
at time $t=1$, community 2 in $G^{[12]}$ becomes inactive while communities 1
and 3 are both active. We consider $S$ in $[20,100]$, and Figure S2 summarizes
the community detection results averaged over 100 data replicates for
Scenarios 1-2, respectively.

(a) $r_{31}(t)$ (b) $r_{32}(t)$ (c) $r_{33}(t)$

Figure 6:  The three time-varying functions $r_{31}(t)$, $r_{32}(t)$ and
$r_{33}(t)$.

(a) Type 1 nodes, Scenario 1 (b) Type 1 nodes, Scenario 2 (c) Type 2 nodes,
Scenario 1 (d) Type 2 nodes, Scenario 2

Figure 7: Average NMIs against the value of $S$ for different methods in
Setting 2 under Scenarios 1-2.

It is seen that DHNet performs better than Methods 1-4 for all values of $S$,
regardless the sparsity of the networks. Interestingly, the performance of
Method 2 in Figure S2 is much worse than that in Figure S1 from Simulation 1
when $S$ is large. This is because the connecting probability is time varying
in Simulation 2, and the signal from communities that are active at different
time points may get ablated in an aggregated picture when $S$ is large.

###  5.3 Simulation setting 3

We consider temporally correlated network samples from a DHSBM with a time-
varying connecting probability matrix. Specifically, we adopt the time-varying
connecting probability $\bm{\Theta}(t)$ from Simulation setting 2. At time
$t_{s}$, the edge $A^{[l_{1}l_{2}]}_{ij}\left(t_{s}\right)$ is a Bernoulli
random variable with

$A_{ij}^{[l_{1}l_{2}]}\left(t_{s}\right)=uA^{[l_{1}l_{2}]}_{ij}\left(t_{s-1}\right)+(1-u)v^{[l_{1}l_{2}]},$

where $u\sim$ iid $\operatorname{Bernoulli}(\alpha)$ and

$v^{[l_{1}l_{2}]}\overset{iid}{\sim}\text{ Bernoulli
}\left(\frac{\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)-\alpha\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)}{1-\alpha}\right),\quad
l_{1},l_{2}\in[L].$

Given $\bm{\Theta}(t)$, a larger $\alpha$ leads to a higher correlation
between the networks at two adjacent time points. We set $S=100$,
$\alpha\in[0,0.4]$ for Scenario 1 and $\alpha\in[0,0.7]$ for Scenario 2, as to
keep the probability parameter
$\frac{\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s}\right)-\alpha\theta^{[l_{1}l_{2}]}_{k_{1}k_{2}}\left(t_{s-1}\right)}{1-\alpha}>0$.
Figure S3 summarizes the community detection results averaged over 100 data
replicates for Scenarios 1-2, respectively.

(a) Type 1 nodes, Scenario 1 (b) Type 1 nodes, Scenario 2 (c) Type 2 nodes,
Scenario 1 (d) Type 2 nodes, Scenario 2

Figure 8: Average NMIs against the value of $\alpha$ for different methods in
Setting 3 under Scenarios 1-2.

Similar conclusions as before can be drawn for all methods shown in Figure S3.
DHNet has the best performance out of the five methods. When $\alpha=0$,
$A_{ij}(t)$ is uncorrelated with the past observations, and the model is
equivalent to the model used in Simulation setting 2. In fact, as $\alpha$
increases, the effective sample size decreases, leading to a deteriorated
performance of DHNet. Method 3 relies only on a random snapshot of network and
as such, it is insensitive to changes in $\alpha$.

##  6 Yelp review network

Yelp is a well-known review website, founded in 2004 in the United States. It
collects reviews on a wide range of businesses such as restaurants, bars and
shops from many countries. On the Yelp platform, users can rate businesses,
submit reviews, and share experiences. We analyze the review data from the
Yelp Challenge (https://www.kaggle.com/yelp-dataset/yelp-dataset) during the
period from January 1st, 2006 to December 31, 2017. This dataset contains a
set of businesses and the category labels of each business (a business usually
has several labels), a set of users and the friendship information among these
users, and the reviews of these businesses by these users. The businesses and
users are anonymized and labeled with numerical identifiers. In this dataset,
the business-category (business is labeled with category) and user-user (user
is friend with user) information are not labeled by time (i.e., not time-
varying) while the user-business (business reviewed by user) interactions are
labeled by time, and a user may review a business several times.

Our analysis focuses on finding heterogeneous communities in the Yelp review
network and predicting interests for new users. Our results show improvements
both in terms of accuracy and interpretability over existing solutions, and
demonstrate the need to consider network heterogeneity and dynamics in
community detection.

###  6.1 Finding heterogeneous communities

To get a comprehensive view of the Yelp review network, we consider a
heterogeneous network with three types of nodes including business, user and
category, connected via three types of edges, user-user (user is friend with
user), user-business (business is reviewed by user), and business-category
(business is labeled with category); see Figure 1 for a simple illustration.
As discussed earlier, the user-user and business-category edges are not time-
varying but the user-business edges are. We focus on businesses that operated
continuously in the study period, business categories that had at least 10
occurrences and users that reviewed at least 20 times in the study period.
This gives a total of 3,566 businesses, 207 categories and 5,116 users, with
141,744 user-user and 17,280 business-category and 194,712 user-business
edges. Due to the high sparsity of user-business edges, we use year as the
time unit when constructing the dynamic network, that is, the network at time
$t$ summarizes the review activity between users and businesses in the $t$-th
year of the study period, and correspondingly $S=12$.

We applied DHNet to the constructed dynamic heterogeneous network with
$\kappa=200$ and identified 11 communities with a maximized modularity value
of 0.237. Table 1 shows the representative categories, number of users and
number of businesses in each identified community, along with a summarizing
theme. The complete list of categories in each community can be found in the
supplement. We found that each community identified by DHNet contains a
distinctive type of businesses. For example, Communities 3-4 are on dining and
Community 7 is mostly on Beauty & Medical. Users in Community 1 prefer
activities related to pets, users in Community 6 prefer bars and entertainment
and users in Community 9 show interests in traveling and sports. Community 11
is mostly on Auto and we did not identify users whose main review activity and
interests are in this type of businesses. These insights can help us
understand the life styles and interests of users in each community.

Table 1: Summary of the 11 communities identified by DHNet. | categories |  | # of
---
users
| # of
---
businesses
theme
1 |  Animal Shelters, Pet Groomers | 719 | 331 | Pets
Pet Services, Veterinarians
2 |  Tex-Mex, Southern | 939 | 408 | Tex-Mex
3 |  Tea, Fast Food | 2629 | 1036 | Casual Dining
Diners, Pizza, Restaurants
4 |  French, Pasta Shops, Steak House | 679 | 337 | Fine Dining
Professional Services, Seafood
5 |  Candy Stores, Farmers Market | 13 | 239 | Stores
Chocolatiers & Shops, Grocery | & Markets
6 |  Adult Entertainment, Bars | 46 | 225 | Bars,
Dance Clubs, Beer Bar | Entertainment
7 |  Beauty & Spas, Doctors | 79 | 441 | Beauty
Hair Salons, Health & Medical | & Medical
8 |  Home Services, Laundry Services | 5 | 289 | Shopping
Music & Video, Shopping Centers | & Life
9 |  Hotels & Travel, Venues & Event Spaces | 4 | 123 | Leisure
Landmarks & Historical Buildings, Tours | & Travel
10 | Buffets, Indian, Pakistani | 3 | 103 | Asian Fusion
11 |  Auto Parts & Supplies, Auto Repair | 0 | 34 | Auto
Automotive, Gas Stations, Tires

We had also applied Methods 1-3 from Section 5, though we did not implement
Method 4, which considers each homogeneous networks separately and discards
information from the edges linking different types of nodes, as there are no
business-business or category-category edges and the user-user edges are not
time-varying. The results from Method 1, which does not distinguish the
different node and edge types, are very difficult to interpret. For example,
one community contains only businesses and one community contains only users.
Method 2, which considers an aggregated heterogeneous network over time, also
identified 11 communities (see details of the communities in the supplement).
The community detection results from Method 2 are less interpretable compared
to DHNet and several communities contain mixed businesses themes. For example,
Hobby Shops is placed into Community 4 that is on fine dining, and Colleges &
Universities and Education are placed into Community 9 that is on leisure and
travel. Method 3, which considers a snapshot of the dynamic work, does not
perform well, as the network at each time point is highly sparse with a large
number of isolated nodes.

###  6.2 Prediction interests for new users

In this section, we aim to predict the interests, in terms of business
categories, for a new Yelp user based on Yelp activities of his/her friends, a
practically useful task in making recommendations and placing advertisements.
We focus on predicting interests in business categories as opposed to
individual businesses, as the number of businesses is large and user-business
interactions are highly sparse. For a new Yelp user, the platform can often
collect his/her friendship information with other existing Yelp users, by
accessing phone contacts, email contacts and Facebook friendship. In terms of
make recommendations, the Yelp activities of friends of a new user can help to
ease the “cold start” problem, the issue where personalized recommendations
cannot be made before a user interacts with the system (e.g., reviewing
businesses).

Consider training and testing datasets taken from two different time periods
(e.g., data from years 2006-2015 as training and years 2015-2017 as testing).
We are interested in making predictions for the new users in the testing set,
which are user accounts that did not exist in the training data. Specifically,
for a new user in the testing set, based on Yelp activities of his/her friends
in the training data, we predict the his/her interests over the business
categories and compare the prediction with the “true” measure calculated from
the testing set.

We compare two different prediction strategies. The first strategy utilizes
community detection results from DHNet in making the prediction and the second
strategy directly averages interests from the new user’s friends without using
any community information, referred to as the naive strategy. The naive
strategy is a commonly adopted practice in recommender systems (Tey et al.,
2021). Specifically, let $g_{i}$ denote the interest measure of the $i$-th new
user in the testing set, which is a probability distribution over all
categories; it is calculated using the appearance frequency of each category
in the businesses reviewed by this user. In the first strategy, we apply DHNet
to the training network data and find the category distribution of each
community, denoted as $f_{j}$ for the $j$-th community, calculated using the
appearance frequency of each category from the businesses in this community.
We then make prediction $g_{i}^{\texttt{DHNet}}$ using the weighted average of
$f_{j}$’s as below

$g_{i}^{\texttt{DHNet}}=\sum_{j}\frac{n_{ij}}{\sum_{j}n_{ij}}f_{j},$

where $n_{ij}$ denotes the number of friends that the $i$-th user has in the
$j$-th community. In this strategy, the prediction is a weighted average of
measures from all communities where the weight reflect the number of
connections the new user has to each community. In the second strategy, we
directly calculate the category distribution $g_{i}^{Naive}$ based on the
businesses that $i$-th user’s friends visited during the training period. This
strategy only focuses on the ego-centric network of the new user and does not
taken into the rich information in the network communities. To assess the
prediction accuracy, we use the Jensen–Shannon divergence (JSD) to compare the
estimated and observed category distributions, that is,

$\operatorname{JSD}(\hat{g_{i}}\|g_{i})=\frac{1}{2}D(\hat{g_{i}}\|m)+\frac{1}{2}D(g_{i}\|m),$

where $\hat{g}_{i}$ refers to the estimated category distribution,
$m=\frac{1}{2}(\hat{g}_{i}+g_{i}$) and $D(\hat{g}_{i}\|m)$ is the
Kullback–Leibler divergence between distributions $\hat{g}_{i}$ and $m$.

Table 2: Jensen–Shannon divergence of the DHNet and naive methods across the 6 moving windows. training years | 2006-2010 | 2007-2011 | 2008-2012 | 2009-2013 | 2010-2014 | 2011-2015
---|---|---|---|---|---|---
testing years | 2011-2012 | 2012-2013 | 2013-2014 | 2014-2015 | 2015-2016 | 2016-2017
$\operatorname{JSD}^{\texttt{DHNet}}$ | 0.122 | 0.133 | 0.132 | 0.136 | 0.135 | 0.138
$\operatorname{JSD}^{\texttt{Naive}}$ | 0.252 | 0.196 | 0.199 | 0.196 | 0.182 | 0.196

Table 2 compares the performance of the two strategies in 6 different sets of
training and testing periods, where

$\operatorname{JSD}^{\texttt{DHNet}}=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\operatorname{JSD}({g}_{i}^{\texttt{DHNet}}\|g_{i}),\quad\operatorname{JSD}^{\texttt{Naive}}=\frac{1}{n_{0}}\sum_{i=1}^{n_{0}}\operatorname{JSD}({g}_{i}^{\texttt{Naive}}\|g_{i})$

and $n_{0}$ denotes the number of new users in the testing period. It is seen
that DHNet outperforms the naive strategy in terms of predicting accuracy in
all training and testing datasets, demonstrating the advantage of utilizing
community structures when predicting user interests.

##  7 Discussion

Maximizing the modularity function as in (5) is not limited to the Louvain-
type method considered in DHNet. Other modularity maximization techniques
developed for a homogeneous network may be applied to (5) with some
modifications, such as the spectral method based on the eigen decomposition of
the modularity matrix or the stochastic optimization method in Massen and Doye
(2005). As noted in modularity maximization for other types of networks
(Fortunato, 2010; Zhang and Chen, 2020), we find that the Louvain-type method
is computationally much more efficient and yields a good performance in our
setting.

While the modularity function value increases at each step of DHNet and the
algorithm is guaranteed to converge, there is no guarantee that it will
converge to the global optimum. Since the modularity maximization problem is
NP-hard, most existing methods are heuristic methods that may only find local
optima and are not guaranteed to find the global optimum. A thorough
theoretical investigation of the local convergence of DHNet can be helpful and
we leave it as future work. Finally, our proposed method can be extended to
weighted and/or directed networks. To incorporate weighted and/or directed
edges into our framework, we need to define a null model for a weighted and/or
directed heterogeneous dynamic network, followed by calculating the
expectations under the null model. This is an interesting topic to investigate
next.

## References

  * Abbe (2017) Abbe, E. (2017), “Community detection and stochastic block models: recent developments,” The Journal of Machine Learning Research, 18, 6446–6531.
  * Blondel et al. (2008) Blondel, V. D., Guillaume, J.-L., Lambiotte, R., and Lefebvre, E. (2008), “Fast unfolding of communities in large networks,” Journal of Statistical Mechanics: Theory and Experiment, 2008, P10008.
  * Brandes et al. (2008) Brandes, U., Delling, D., Gaertler, M., Gorke, R., Hoefer, M., Nikoloski, Z., and Wagner, D. (2008), “On modularity clustering,” IEEE Transactions on Knowledge and Data Engineering, 20, 172–188.
  * Chung et al. (2006) Chung, F., Fan, R., Chung, F. R., Graham, F. C., Lu, L., Chung, K. F., et al. (2006), Complex Graphs and Networks, no. 107, American Mathematical Soc.
  * Clauset et al. (2004) Clauset, A., Newman, M. E., and Moore, C. (2004), “Finding community structure in very large networks,” Physical Review E, 70, 066111.
  * Danon et al. (2005) Danon, L., Diaz-Guilera, A., Duch, J., and Arenas, A. (2005), “Comparing community structure identification,” Journal of Statistical Mechanics: Theory and Experiment, 2005, P09008.
  * Fortunato (2010) Fortunato, S. (2010), “Community detection in graphs,” Physics Reports, 486, 75–174.
  * Guimera et al. (2004) Guimera, R., Sales-Pardo, M., and Amaral, L. A. N. (2004), “Modularity from fluctuations in random graphs and complex networks,” Physical Review E, 70, 025101.
  * Jiang et al. (2021) Jiang, S., Koch, B., and Sun, Y. (2021), “HINTS: Citation time series prediction for new publications via dynamic heterogeneous information network embedding,” in Proceedings of the Web Conference 2021, pp. 3158–3167.
  * Linden et al. (2003) Linden, G., Smith, B., and York, J. (2003), “Amazon. com recommendations: Item-to-item collaborative filtering,” IEEE Internet Computing, 7, 76–80.
  * Massen and Doye (2005) Massen, C. P. and Doye, J. P. (2005), “Identifying communities within energy landscapes,” Physical Review E, 71, 046101.
  * Moody and White (2003) Moody, J. and White, D. R. (2003), “Structural cohesion and embeddedness: A hierarchical concept of social groups,” American Sociological Review, 68, 103–127.
  * Newman (2006) Newman, M. E. (2006), “Finding community structure in networks using the eigenvectors of matrices,” Physical Review E, 74, 036104.
  * Newman and Girvan (2004) Newman, M. E. and Girvan, M. (2004), “Finding and evaluating community structure in networks,” Physical Review E, 69, 026113.
  * Sengupta and Chen (2015) Sengupta, S. and Chen, Y. (2015), “Spectral clustering in heterogeneous networks,” Statistica Sinica, 25, 1081–1106.
  * Sørlie et al. (2001) Sørlie, T., Perou, C. M., Tibshirani, R., Aas, T., Geisler, S., Johnsen, H., Hastie, T., Eisen, M. B., Van De Rijn, M., Jeffrey, S. S., et al. (2001), “Gene expression patterns of breast carcinomas distinguish tumor subclasses with clinical implications,” Proceedings of the National Academy of Sciences, 98, 10869–10874.
  * Sun et al. (2010) Sun, Y., Tang, J., Han, J., Gupta, M., and Zhao, B. (2010), “Community evolution detection in dynamic heterogeneous information networks,” in Proceedings of the Eighth Workshop on Mining and Learning with Graphs, pp. 137–146.
  * Tey et al. (2021) Tey, F. J., Wu, T.-Y., Lin, C.-L., and Chen, J.-L. (2021), “Accuracy improvements for cold-start recommendation problem using indirect relations in social networks,” Journal of Big Data, 8, 1–18.
  * Wakita and Tsurumi (2007) Wakita, K. and Tsurumi, T. (2007), “Finding community structure in mega-scale social networks,” in Proceedings of the 16th International Conference on World Wide Web, pp. 1275–1276.
  * Wang et al. (2022) Wang, X., Lu, Y., Shi, C., Wang, R., Cui, P., and Mou, S. (2022), “Dynamic heterogeneous information network embedding with meta-path based proximity,” IEEE Transactions on Knowledge and Data Engineering, 34, 1117 – 1132.
  * Xue et al. (2020) Xue, H., Yang, L., Jiang, W., Wei, Y., Hu, Y., and Lin, Y. (2020), “Modeling dynamic heterogeneous network for link prediction using hierarchical attention with temporal RNN,” arXiv preprint arXiv:2004.01024.
  * Yin et al. (2019) Yin, Y., Ji, L.-X., Zhang, J.-P., and Pei, Y.-L. (2019), “DHNE: Network representation learning method for dynamic heterogeneous networks,” IEEE Access, 7, 134782–134792.
  * Zhang and Cao (2017) Zhang, J. and Cao, J. (2017), “Finding common modules in a time-varying network with application to the Drosophila Melanogaster gene regulation network,” Journal of the American Statistical Association, 112, 994–1008.
  * Zhang and Chen (2017) Zhang, J. and Chen, Y. (2017), “A hypothesis testing framework for modularity based network community detection,” Statistica Sinica, 27, 437–456.
  * Zhang and Chen (2020) Zhang, J. and Chen, Y. (2020), “Modularity based community detection in heterogeneous networks,” Statistica Sinica, 30, 601–629.
  * Zhang et al. (2020) Zhang, J., Sun, W. W., and Li, L. (2020), “Mixed-effect time-varying network model and application in brain connectivity analysis,” Journal of the American Statistical Association, 115, 2022–2036.
  * Zhang et al. (2022) Zhang, Z., Huang, J., and Tan, Q. (2022), “Multi-view Dynamic heterogeneous information network embedding,” The Computer Journal, 65, 2016–2033.

Supplementary Materials for “Fast Community Detection in Dynamic and
Heterogeneous Networks"

Maoyu Zhang, Jingfei Zhang and Wenlin Dai

This supplementary material gives the proof of Theorem 1 in Section S1,
additional simulation results in Section S2 and additional real data analysis
results in Section S3.

##  S1 Proof of Theorem 1

First, we formalize the notations that will be used in the proof. Consider a
dynamic heterogeneous network
$\mathcal{G}\left(\bigcup_{i=1}^{L}V^{[i]},\mathcal{E}(t)\cup\mathcal{E}^{+}(t)\right)$,
let $G^{[l]}(t)$ denote the homogeneous network formed within node set
$V^{[l]}$ with an $n_{l}\times n_{l}$ adjacency matrix $A^{[l]}(t)$ and
$G^{[l_{1}l_{2}]}(t)=(V^{[l_{1}]}\cup V^{[l_{2}]},E^{[l_{1}l_{2}]}(t))$ denote
the bi-partite network formed between node sets $V^{[l_{1}]}$ and
$V^{[l_{2}]}$ with an $n_{l_{1}}\times n_{l_{2}}$ bi-adjacency matrix
$A^{[l_{1}l_{2}]}(t)$ at time $t$, $l_{1},l_{2}\in[L]$. Write the number of
edges in $A^{[l]}(t)$ and $A^{[l_{1}l_{2}]}(t)$ as
$m^{[l]}(t)=\sum_{i,j}A_{ij}^{l}(t)/2$ and
$m^{[l_{1}l_{2}]}(t)=\sum_{i,j}A_{ij}^{l_{1}l_{2}}(t)$, respectively. For a
dynamic heterogeneous network $\\{\mathcal{G}(t_{s}),s\in\mathcal{S}\\}$ from
the DHSBM model, each $A_{ij}^{[l]}(t_{s})$ and
$A_{ij}^{\left[l_{1}l_{2}\right]}(t_{s})$ are independent Bernoulli random
variables with

$E\left(A_{ij}^{[l]}(t_{s})\mid
c_{i}^{[l]}=k,c_{j}^{[l]}=h\right)=\theta_{kh}^{[l]}(t_{s}),\text{ and
}E\left(A_{ij}^{\left[l_{1}l_{2}\right]}(t_{s})\mid
c_{i}^{\left[l_{1}\right]}=k,c_{j}^{\left[l_{2}\right]}=h\right)=\theta_{kh}^{\left[l_{1}l_{2}\right]}(t_{s}).$

For a community assignment label
$\bm{e}=\left(\mathbf{e}^{[1]},\ldots,\mathbf{e}^{[L]}\right)$ with
$\mathbf{e}^{[l]}=\left(e_{1}^{[l]},\ldots,e_{n_{l}}^{[l]}\right),l\in[L]$,
define $K\times K$ matrices $O^{[l],s},l\in[L]$, and
$O^{\left[l_{1}l_{2}\right],s}$ $1\leq l_{1}\neq l_{2}\leq L$, such that

$\begin{gathered}O_{kh}^{[l],s}(\bm{e})=\sum_{ij}A_{ij}^{[l]}(t_{s})I\left(e_{i}^{[l]}=k,e_{j}^{[l]}=h\right),\\\
O_{kh}^{\left[l_{1}l_{2}\right],s}(\bm{e})=\sum_{ij}A_{ij}^{\left[l_{1}l_{2}\right]}(t_{s})I\left(e_{i}^{\left[l_{1}\right]}=k,e_{j}^{\left[l_{2}\right]}=h\right).\end{gathered}$

Define $O_{k}^{[l],s}=\sum_{h}O_{kh}^{[l],s}$ and
$O_{k}^{\left[l_{1}l_{2}\right],s}=\sum_{h}O_{kh}^{\left[l_{l}l_{2}\right],s},l\in[L],1\leq
l_{1}\neq l_{2}\leq L.$ Define $K\times K$ matrices
$R^{[l]}(\bm{e}),V^{[l]}(\bm{e}),l\in[L]$, such that

$R_{ab}^{[l]}(\bm{e})=\frac{1}{n}\sum_{i=1}^{n_{l}}I\left(e_{i}^{[l]}=a,c_{i}^{[l]}=b\right),\text{
and
}V_{ab}^{[l]}(\bm{e})=\frac{\sum_{i=1}^{n_{l}}I\left(e_{i}^{[l]}=a,c_{i}^{[l]}=b\right)}{\sum_{i=1}^{n_{l}}I\left(c_{i}^{[l]}=b\right)}.$

Write $\mathbb{O}(\bm{e})=\left(\mathcal{O}^{s}(\bm{e}),s=1,\ldots,S\right)$,
where
$\mathcal{O}^{s}(\bm{e})=\left\\{O^{{[l]},s},O^{{\left[l_{1}l_{2}\right]},s},l\in[L],1\leq
l_{1}\neq l_{2}\leq L\right\\}$, and
$\mathcal{R}=\left\\{R^{[1]},\ldots,R^{[L]}\right\\}$.

The modularity function $Q(\bm{e},\\{\mathcal{G}(t_{s})\\}_{s\in[S]})$ can be
expressed as

$\frac{1}{2\bar{m}^{[l]}L^{2}}\sum_{s=1}^{S}\sum_{l=1}^{L}\sum_{k=1}^{K}\left(O_{kk}^{[l],s}-\frac{(O_{k}^{[l],s})^{2}}{\sum_{kh}O_{kh}^{[l],s}}\right)+\frac{1}{\bar{m}^{[l_{1}l_{2}]}L^{2}}\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}\sum_{k=1}^{K}\left(O_{kk}^{{\left[l_{1}l_{2}\right]},s}-\frac{O_{k}^{{\left[l_{1}l_{2}\right]},s}O_{k}^{{\left[l_{2}l_{1}\right]},s}}{\sum_{kh}O_{kh}^{{\left[l_{1}l_{2}\right]},s}}\right),$

where $\bar{m}^{[l]}=\sum_{s=1}^{S}m^{[l]}(t_{s})$ and
$\bar{m}^{[l_{1}l_{2}]}=\sum_{s=1}^{S}m^{[l_{1}l_{2}]}(t_{s})$.

Define $\mu_{n,S}=n^{2}S\rho_{n,S}$, we have:

$\displaystyle\frac{1}{\mu_{n,S}}E\left(O_{kh}^{{\left[l_{1}l_{2}\right]},s}(\bm{e})\mid\bm{c}\right)$
$\displaystyle=\frac{1}{\mu_{n,S}}E\left(\sum_{ij}A_{ij}^{\left[l_{1}l_{2}\right]}(t_{s})I\left(e_{i}^{\left[l_{1}\right]}=k,e_{j}^{\left[l_{2}\right]}=h\right)\mid\bm{c}\right)$
$\displaystyle=\frac{1}{n^{2}S}\sum_{ij}\sum_{ab}\theta_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})I\left(e_{i}^{\left[l_{1}\right]}=k,c_{i}^{\left[l_{1}\right]}=a\right)I\left(e_{j}^{\left[l_{2}\right]}=h,e_{j}^{\left[l_{2}\right]}=b\right).$

Define
$H^{{\left[l_{1}l_{2}\right]},s}(\mathcal{R}(\bm{e}))=\frac{1}{\mu_{n,S}}E\left(O^{{\left[l_{1}l_{2}\right]},s}(\bm{e})\mid\bm{c}\right)$,
we have

$H^{{\left[l_{1}l_{2}\right]},s}(\mathcal{R}(\bm{e}))=\frac{1}{S}R^{\left[l_{1}\right]}(\bm{e})\theta^{\left[l_{1}l_{2}\right]}(t_{s})R^{\left[l_{2}\right]}(\bm{e})^{\prime},\quad
1\leq l_{1}\neq l_{2}\leq L.$

Similarly, we can define
$H^{{[l]},s}(\mathcal{R}(\bm{e}))=\frac{1}{\mu_{n,S}}E\left(O^{[l],s}(\bm{e})\mid\bm{c}\right)$
and write

$H^{{[l]},s}(\mathcal{R}(\bm{e}))=\frac{1}{S}R^{[l]}(\bm{e})\theta^{[l]}(t_{s})R^{[l]}(\bm{e})^{\prime},\quad
l\in[L],s\in[S].$

Write
$\mathcal{H}=\left\\{H^{{[l]},1},H^{{\left[l_{1}l_{2}\right]},S},\dots,H^{{[l]},1},H^{{\left[l_{1}l_{2}\right]},S}l\in[L],1\leq
l_{1}\neq l_{2}\leq L\right\\}$.

Consider a community label
$\bm{e}=\left(\mathbf{e}^{[1]},\ldots,\mathbf{e}^{[L]}\right)$ with
$\mathbf{e}^{[l]}=\left(e_{1}^{[l]},\ldots,e_{n_{l}}^{[l]}\right),l\in[L]$.
Further, define

$\mathcal{J}(\mathbb{O}(\bm{e}))=\frac{1}{S}\sum_{s=1}^{S}J\left(\mathcal{O}^{s}(\bm{e})\right),$

where

$J\left(\mathcal{O}^{s}(\bm{e})\right)=\frac{1}{L^{2}}\left[\sum_{l=1}^{L}\sum_{k=1}^{K}\left(O_{kk}^{[l],s}-\frac{(O_{k}^{[l],s})^{2}}{\sum_{kh}O_{kh}^{[l],s}}\right)+\sum_{l_{1}\neq
l_{2}}^{L}\sum_{k=1}^{K}\left(O_{kk}^{{\left[l_{1}l_{2}\right]},s}-\frac{O_{k}^{{\left[l_{1}l_{2}\right]},s}O_{k}^{{\left[l_{2}l_{1}\right]},s}}{\sum_{kh}O_{kh}^{{\left[l_{1}l_{2}\right]},s}}\right)\right].$

Here we suppress the argument $\bm{e}$ for brevity. Then for convenient, we
write

$J(\mathcal{O}^{s})=\sum_{l=1}^{L}J_{1}\left(O^{{[l]},s}\right)+\sum_{l_{1}\neq
l_{2}}^{L}J_{2}\left(O^{{\left[l_{1}l_{2}\right]},s},O^{{\left[l_{2}l_{1}\right]},s}\right),$

where

$J_{1}\left(O^{{[l]},s}\right)=\sum_{k=1}^{K}\left(O_{kk}^{{[l]},s}-\frac{(O_{k}^{{[l]},s})^{2}}{\sum_{kh}O_{kh}^{[l],s}}\right),$

and

$J_{2}\left(O^{\left[l_{1}l_{2}\right],s},O^{\left[l_{2}l_{1}\right],s}\right)=\sum_{k=1}^{K}\left(O_{kk}^{\left[l_{1}l_{2}\right],s}-\frac{O_{k}^{\left[l_{1}l_{2}\right],s}O_{k}^{\left[l_{2}l_{1}\right],s}}{\sum_{kh}O_{kh}^{\left[l_{1}l_{2}\right],s}}\right).$

Showing the $\hat{\bm{e}}$ that maximizes the
$Q(\bm{e},\\{\mathcal{G}(t_{s})\\}_{s\in[S]})$ is consistent is equivalent to
showing the $\hat{\bm{e}}$ that maximizes the
$\mathcal{J}(\mathbb{O}(\bm{e}))$ is consistent. We show consistency by
showing that there exists $\delta_{n,S}\rightarrow 0$, such that

$P\left(\max_{\bm{e}:\eta(\bm{e},\bm{c})\geq\delta_{n,S}}\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)<\mathcal{J}\left(\frac{\mathbb{O}(\bm{c})}{\mu_{n,S}}\right)\right)\rightarrow
1\text{ as }nS\rightarrow\infty,$

where
$\eta(\bm{e},\bm{c})=\sum_{l=1}^{L}\sum_{ab}\left|V_{ab}^{[l]}(\bm{e})-V_{ab}^{[l]}(\bm{c})\right|$.

Since $\mathcal{J}(.)$ is Lipschitz in all its arguments, we have

$\displaystyle\left|\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)-\mathcal{J}(\mathcal{H}(\mathcal{R}))\right|$
$\displaystyle\leq
M_{1}\max_{l}\left\|\sum_{s=1}^{S}\frac{O^{[l],s}(\mathcal{E})}{\mu_{n,S}}-\sum_{s=1}^{S}H^{[l],s}(\mathcal{R})\right\|_{\infty}$
$\displaystyle+M_{1}\max_{l_{1}\neq
l_{2}}\left\|\sum_{s=1}^{S}\frac{O^{\left[l_{1}l_{2}\right],s}(\mathcal{E})}{\mu_{n,S}}-\sum_{s=1}^{S}H^{\left[l_{1}l_{2}\right],s}(\mathcal{R})\right\|_{\infty}.$

Here $\|X\|_{\infty}=\max_{kh}\left|X_{kh}\right|$. To continue with the
proof, we need to use the Bernstein’s inequality, Lemma A.1 of
zhao2012consistency.
Bernstein’s inequality Let $X_{1},\dots,X_{n}$ be independent variables.
Suppose that $|X_{i}|\leq M$ for all $i$. Then, for all positive $t$

$P\left(\left|\sum_{i=1}^{n}X_{i}-\sum_{i=1}^{n}E\left(X_{i}\right)\right|>t\right)\leq
2\exp\left(-\frac{t^{2}/2}{\sum\operatorname{var}\left(X_{i}\right)+Mt/3}\right).$

Define $\tau=\max_{ij,s}\operatorname{var}\left(A_{ij}^{[l]}(t_{s})\right)$.
For any $\epsilon<3\tau$, if we write $\omega=\epsilon n^{2}S\rho_{n,S}$, we
have

$\displaystyle
P\left(\left|\sum_{s=1}^{S}\frac{O_{kh}^{[l],s}(\mathcal{E})}{\mu_{n,S}}-\sum_{s=1}^{S}H_{kh}^{[l],s}(\mathcal{R})\right|>\epsilon\right)$
$\displaystyle\leq
2\exp\left(-\frac{\omega^{2}/2}{\operatorname{var}\left(\sum_{s=1}^{S}O_{kh}^{[l],s}(\mathcal{E})\right)+2\omega/3}\right)$
$\displaystyle\leq
2\exp\left(-\frac{\epsilon^{2}n^{4}S^{2}\rho_{n}^{2}}{8n^{2}S\rho_{n}\tau}\right)$
$\displaystyle=2\exp\left(-\frac{\epsilon^{2}\mu_{n,S}}{8\tau}\right)$

Notice that
$\operatorname{var}\left(\sum_{s=1}^{S}O_{kh}^{[l],s}(\mathcal{E})\right)\leq
2n^{2}S\max_{ij}\operatorname{var}\left(A_{ij}^{[l]}(t_{s})\right)$.

The left hand side of the inequality converges to $0$ in probability uniformly
over $\bm{e}$ as $nS\rho_{n,S}\rightarrow\infty$. Following similar arguments,
we can show that

$P\left(\left|\sum_{s=1}^{S}\frac{O_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{E})}{\mu_{n,S}}-\sum_{s=1}^{S}H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{R})\right|>\epsilon\right)\rightarrow
0\text{ as }nS\rightarrow\infty.$

Therefore $\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)$ is
uniformly close to $\mathcal{J}(\mathcal{H}(\mathcal{R}(\bm{e})))$, i.e.,
there exists $\epsilon_{n,S}\rightarrow 0$ such that

$P\left(\max_{\bm{e}}|\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)-\mathcal{J}(\mathcal{H}(\mathcal{R}(\bm{e})))|<\epsilon_{n,S}\right)\rightarrow
1\text{ as }nS\rightarrow\infty.$ (S1)

To show that there exists $\delta_{n,S}\rightarrow 0$, such that

$P\left(\max_{\bm{e}:\eta(\bm{e},\bm{c})\geq\delta_{n,S}}\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)<\mathcal{J}\left(\frac{\mathbb{O}(\bm{c})}{\mu_{n,S}}\right)\right)\rightarrow
1\text{ as }nS\rightarrow\infty.$

Next we show that $\mathcal{J}(\mathcal{H}(\mathcal{R}(\bm{e})))$ is uniquely
maximized over $\left\\{\mathcal{R}:R^{[l]}\geq
0,R^{[l]^{\prime}}\mathbf{1}=\pi^{[l]},l=\right.$ $1,\ldots,L\\}$ at
$\mathcal{S}=\mathcal{R}(\bm{c})$. Since
$\mathcal{J}(\mathcal{H}(\mathcal{R}))$ is the population version of
$\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)$, if
$\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)$ is maximized by
the true community label $\bm{c},J(\mathcal{H}(\mathcal{R}))$ should also be
maximized by the true assignment $\mathcal{S}$. Define

$\triangle_{kh}=\left\\{\begin{array}[]{rll}1&\text{ for }&k=h\\\ -1&\text{
for }&k\neq h\end{array}\right.$

Using the equalities

$\sum_{k}\left(H_{kk}^{[l],s}-\frac{(H_{k}^{[l],s})^{2}}{\sum_{kh}H_{kh}^{[l],s}}\right)+\sum_{k\neq
h}\left(H_{kh}^{[l],s}-\frac{H_{k}^{[l],s}H_{h}^{[l],s}}{\sum_{kh}H_{kh}^{[l],s}}\right)=0,\quad
l\in[L],s=1,\dots,S,$

and

$\sum_{k}\left(H_{kk}^{\left[l_{1}l_{2}\right],s}-\frac{H_{k}^{\left[l_{1}l_{2}\right],s}H_{k}^{\left[l_{2}l_{1}\right],s}}{\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}}\right)+\sum_{k\neq
h}\left(H_{kh}^{\left[l_{1}l_{2}\right],s}-\frac{H_{k}^{\left[l_{1}l_{2}\right],s}H_{h}^{\left[l_{2}l_{1}\right],s}}{\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}}\right)=0,\quad
1\leq l_{1}\neq l_{2}\leq L,s=1,\dots,S.$

We have

$\displaystyle\mathcal{J}(\mathcal{H}(\mathcal{R}))=$
$\displaystyle\sum_{s=1}^{S}\sum_{l=1}^{L}J_{1}\left(H^{[l],s}(\mathcal{R})\right)+\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}J_{2}\left(H^{\left[l_{1}l_{2}\right],s}(\mathcal{R}),H^{\left[l_{2}l_{1}\right],s}(\mathcal{R})\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{s=1}^{S}\sum_{l=1}^{L}\sum_{kh}\triangle_{kh}\left(H_{kh}^{[l],s}(\mathcal{R})-\frac{H_{k}^{[l],s}(\mathcal{R})H_{h}^{[l],s}(\mathcal{R})}{\sum_{kh}H_{kh}^{[l],s}(\mathcal{R})}\right)+$
$\displaystyle\frac{1}{2}\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}\sum_{kh}\triangle_{kh}\left(H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{R})-\frac{H_{k}^{\left[l_{1}l_{2}\right],s}(\mathcal{R})H_{h}^{\left[l_{2}l_{1}\right],s}(\mathcal{R})}{\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{R})}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2S}\sum_{s=1}^{S}\sum_{l=1}^{L}\sum_{kh}\triangle_{kh}\left(\sum_{ab}\theta_{ab}^{[l]}(t_{s})R_{ka}^{[l]}(\bm{e})R_{hb}^{[l]}(\bm{e})-\frac{\left(\sum_{aq}\theta_{aq}^{[l]}(t_{s})R_{ka}^{[l]}(\bm{e})\pi_{q}^{[l]}\right)\left(\sum_{br}\theta_{br}^{[l]}(t_{s})R_{hb}^{[l]}(\bm{e})\pi_{r}^{[l]}\right)}{\sum_{kh}H_{kh}^{[l],s}(\mathcal{R})}\right)$
$\displaystyle+$ $\displaystyle\frac{1}{2S}\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}\sum_{kh}\triangle_{kh}\left(\sum_{ab}\theta_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})R_{ka}^{\left[l_{l}\right]}(\bm{e})R_{hb}^{\left[l_{2}\right]}(\bm{e})\right.$
$\displaystyle\hskip
180.67499pt-\left.\frac{\left(\sum_{aq}\theta_{aq}^{\left[l_{1}l_{2}\right]}(t_{s})R_{ka}^{\left[l_{1}\right]}(\bm{e})\pi_{q}^{\left[l_{2}\right]}\right)\left(\sum_{br}\theta_{br}^{\left[l_{2}l_{1}\right]}(t_{s})R_{hb}^{\left[l_{2}\right]}(\bm{e})\pi_{r}^{\left[l_{1}\right]}\right)}{\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{R})}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{2S}\sum_{s=1}^{S}\sum_{l=1}^{L}\sum_{kh}\sum_{ab}\triangle_{kh}R_{ka}^{[l]}(\bm{e})R_{hb}^{[l]}(\bm{e})\left(\theta_{ab}^{[l]}(t_{s})-\frac{\left(\sum_{q}\theta_{aq}^{[l]}(t_{s})\pi_{q}^{[l]}\right)\left(\sum_{r}\theta_{br}^{[l]}(t_{s})\pi_{r}^{[l]}\right)}{\sum_{kh}H_{kh}^{[l]}(\mathcal{R})}\right)$
$\displaystyle+$ $\displaystyle\frac{1}{2S}\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}\sum_{kh}\sum_{ab}\triangle_{kh}R_{ka}^{\left[l_{1}\right]}(\bm{e})R_{hb}^{\left[l_{2}\right]}(\bm{e})\left(\theta_{ab}^{\left[l_{l}l_{2}\right]}(t_{s})-\frac{\left(\sum_{q}\theta_{aq}^{\left[l_{1}l_{2}\right]}(t_{s})\pi_{q}^{\left[l_{2}\right]}\right)\left(\sum_{r}\theta_{br}^{\left[l_{2}l_{1}\right]}(t_{s})\pi_{r}^{\left[l_{1}\right]}\right)}{\sum_{kh}H_{kh}^{\left[l_{l}l_{2}\right],s}(\mathcal{R})}\right)$
$\begin{aligned}
\leq&\frac{1}{2S}\sum_{s=1}^{S}\sum_{l=1}^{L}\sum_{kh}\sum_{ab}\triangle_{ab}R_{ka}^{[l]}(\bm{e})R_{hb}^{[l]}(\bm{e})\left(\theta_{ab}^{[l]}(t_{s})-\frac{\left(\sum_{q}\theta_{aq}^{[l]}(t_{s})\pi_{q}^{[l]}\right)\left(\sum_{r}\theta_{br}^{[l]}(t_{s})\pi_{r}^{[l]}\right)}{\sum_{kh}H_{kh}^{[l]}(\mathcal{R})}\right)\\\
+&\frac{1}{2S}\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}\sum_{kh}\sum_{ab}\triangle_{ab}R_{ka}^{\left[l_{1}\right]}(\bm{e})R_{hb}^{\left[l_{2}\right]}(\bm{e})\left(\theta_{ab}^{\left[l_{l}l_{2}\right]}(t_{s})-\frac{\left(\sum_{q}\theta_{aq}^{\left[l_{1}l_{2}\right]}(t_{s})\pi_{q}^{\left[l_{2}\right]}\right)\left(\sum_{r}\theta_{br}^{\left[l_{2}l_{1}\right]}(t_{s})\pi_{r}^{\left[l_{1}\right]}\right)}{\sum_{kh}H_{kh}^{\left[l_{l}l_{2}\right],s}(\mathcal{R})}\right)\\\
=&\frac{1}{2S}\sum_{s=1}^{S}\sum_{l=1}^{L}\sum_{ab}\triangle_{ab}\pi_{a}^{[l]}\pi_{b}^{[l]}\left(\theta_{ab}^{[l]}(t_{s})-\frac{\left(\sum_{q}\theta_{aq}^{[l]}(t_{s})\pi_{q}^{[l]}\right)\left(\sum_{r}\theta_{br}^{[l]}(t_{s})\pi_{r}^{[l]}\right)}{\sum_{kh}H_{kh}^{[l],s}(\mathcal{S})}\right)\\\
+&\frac{1}{2S}\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}^{L}\sum_{ab}\triangle_{ab}\pi_{a}^{\left[l_{1}\right]}\pi_{b}^{\left[l_{2}\right]}\left(\theta_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})-\frac{\left(\sum_{q}\theta_{aq}^{\left[l_{2}l_{1}\right]}(t_{s})\pi_{q}^{\left[l_{2}\right]}\right)\left(\sum_{r}\theta_{br}^{\left[l_{1}l_{2}\right]}(t_{s})\pi_{r}^{\left[l_{1}\right]}\right)}{\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{S})}\right)\\\
=&\sum_{s=1}^{S}\sum_{l=1}^{L}J_{1}\left(H^{[l],s}(\mathcal{S})\right)+\sum_{s=1}^{S}\sum_{l_{1}\neq
l_{2}}J_{2}\left(H^{\left[l_{1}l_{2}\right],s}(\mathcal{S}),H^{\left[l_{2}l_{1}\right],s}(\mathcal{S})\right)=\mathcal{J}(\mathcal{H}(\mathcal{S}))\end{aligned}.$

Here we used the conditions in Theorem 1 for the inequality, and the
relationship that

$\sum_{kh}H_{kh}^{[l],s}(\mathcal{R})=\frac{1}{S}\sum_{kh}\sum_{ab}\theta_{ab}^{[l]}(t_{s})R_{ka}^{[l]}(\bm{e})R_{hb}^{[l]}(\bm{e})=\frac{1}{S}\sum_{ab}\theta_{ab}^{[l]}(t_{s})\pi_{a}^{[l]}\pi_{b}^{[l]}=\sum_{kh}H_{kh}^{[l],s}(\mathcal{S}),$

and

$\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{R})=\frac{1}{S}\sum_{kh}\sum_{ab}\theta_{ab}^{\left[l_{1}l_{2}\right]}(t_{s})R_{ka}^{\left[l_{1}\right]}(\bm{e})R_{hb}^{\left[l_{2}\right]}(\bm{e})=\frac{1}{S}\sum_{ab}\theta_{ab}^{\left[l_{l}l_{2}\right]}(t_{s})\pi_{a}^{\left[l_{1}\right]}\pi_{b}^{\left[l_{2}\right]}=\sum_{kh}H_{kh}^{\left[l_{1}l_{2}\right],s}(\mathcal{S}).$

We have shown that $\mathcal{S}$ is a maximizer of
$\mathcal{J}(\mathcal{H}(\mathcal{R}))$.

Next we need to show that $\mathcal{S}$ is the unique maximizer of
$\mathcal{J}(\mathcal{H}(\mathcal{R}))$. This can be shown using Lemma $3.2$
in bickel2009nonparametric. Since the inequality
$\mathcal{J}(\mathcal{H}(\mathcal{R}))\leq\mathcal{J}(\mathcal{H}(\mathcal{S}))$
holds only if $\triangle_{kh}=\triangle_{ab}$ whenever
$R_{ka}^{[l]}(\bm{e})R_{hb}^{[l]}(\bm{e})>0,l\in[L]$, and $\triangle$ does not
have two identical columns, using the results in Lemma 3.2, we have
$\mathcal{S}$ uniquely maximizes $J(\mathcal{H}(\mathcal{R}))$. Now that we
have shown that $\mathcal{J}(\mathcal{H}(\mathcal{R}))$ is uniquely maximized
by $\mathcal{S}$. By the continuity of $\mathcal{J}(.)$ in the neighborhood of
$\mathcal{S}$, there exists $\delta_{n,S}\rightarrow\infty$, such that

$J(\mathcal{H}(\mathcal{R}))-J(\mathcal{H}(\mathcal{S}))\geq
2\epsilon_{n,S}\quad\text{ for }\quad\eta(\bm{e},\bm{c})\geq\delta_{n,S}.$

Here we used the fact that

$\displaystyle\eta(\mathcal{R}(\bm{e}),\mathcal{S})$
$\displaystyle=\sum_{l=1}^{L}\sum_{ab}\left|\pi_{b}^{[l]}V_{ab}^{[l]}(\bm{e})-\pi_{b}^{[l]}V_{ab}^{[l]}(\bm{c})\right|$
$\displaystyle\geq\left(\min_{l,b}\pi_{b}^{[l]}\right)\times\sum_{l=1}^{L}\sum_{ab}\left|V_{ab}^{[l]}(\bm{e})-V_{ab}^{[l]}(\bm{c})\right|=\left(\min_{l,b}\pi_{b}^{[l]}\right)\times\eta(\bm{e},\bm{c}).$

Thus, with (S1), we have that

$\displaystyle
P\left(\max_{\bm{e}:\eta(\bm{e},\bm{c})\geq\delta_{n,s}}\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)<\mathcal{J}\left(\frac{\mathbb{O}(\bm{c})}{\mu_{n,S}}\right)\right)$
$\displaystyle\geq
P\left(\left|\max_{\bm{e}:\eta(\bm{e},\bm{c})\geq\delta_{n,S}}\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right)-\max_{\bm{e}:\eta(\bm{e},\bm{c})\geq\delta_{n,S}}\mathcal{J}(\mathcal{H}(\mathcal{R}))\right|<\epsilon_{n,S},\left|\mathcal{J}\left(\frac{\mathbb{O}(\bm{c})}{\mu_{n,S}}\right)-\mathcal{J}(\mathcal{H}(\mathcal{S}))\right|\leq\epsilon_{n,S}\right)\rightarrow
1,$

and this implies that

$P\left(\eta(\hat{\bm{c}},\bm{c})\leq\delta_{n,S}\right)\rightarrow 1,$

where

$\hat{\bm{c}}=\arg\max_{\bm{e}}\mathcal{J}\left(\frac{\mathbb{O}(\bm{e})}{\mu_{n,S}}\right),$

since

$\displaystyle\frac{1}{n}\sum_{l=1}^{L}\sum_{i=1}^{n_{l}}I\left(\hat{c}_{i}^{[l]}\neq
c_{i}^{[l]}\right)=\sum_{l=1}^{L}\sum_{k}\pi_{k}^{[l]}\left(1-V_{kk}^{[l]}(\hat{\bm{c}})\right)$
$\displaystyle\leq\sum_{l}^{L}\sum_{k}\left(1-V_{kk}^{[l]}(\hat{\bm{c}})\right)$
$\displaystyle=\frac{1}{2}\sum_{l=1}^{L}\left(\sum_{k}\left(1-V_{kk}^{[l]}(\hat{\bm{c}})\right)+\sum_{k\neq
h}V_{kh}^{[l]}(\hat{\bm{c}})\right)$
$\displaystyle=\eta(\hat{\bm{c}},\bm{c})/2.$

We have thus established the consistency property of $\hat{\bm{c}}$.

##  S2 Additional simulation results

In this section, we provide some additional simulation results for simulation
settings 1-3, where the network generation is the same as in the simulation
section of the text, except that we have also considered the case where
$G^{[1]}$ has a weak community structure while $G^{[2]}$ has no community
structure, leading to $r_{1}=0.05$ in Scenarios S1 and S2.
Scenario S1: $\theta_{1}=0.5$, $\theta_{2}=0.6$, $\theta_{3}=0.3$,
$r_{1}=0.05$, $r_{2}=0$,
Scenario S2: $\theta_{1}=0.1$, $\theta_{2}=0.2$, $\theta_{3}=0.05$,
$r_{1}=0.05$, $r_{2}=0$.

(a) Type 1 nodes, Scenario S1 (b) Type 1 nodes, Scenario S2

(c) Type 2 nodes, Scenario S1 (d) Type 2 nodes, Scenario S2

Figure S1: Average NMIs against the value of $r_{3}$ for different methods in
Setting 1 under Scenarios S1-S2.

(a) Type 1 nodes, Scenario S1 (b) Type 1 nodes, Scenario S2

(c) Type 2 nodes, Scenario S1 (d) Type 2 nodes, Scenario S2

Figure S2: Average NMIs against the value of $S$ for different methods in
Setting 2 under Scenarios S1-S2.

(a) Type 1 nodes, Scenario 1 (b) Type 1 nodes, Scenario 2

(c) Type 2 nodes, Scenario 1 (d) Type 2 nodes, Scenario 2

Figure S3: Average NMIs against the value of $\alpha$ for different methods in
Setting 3 under Scenarios 1-2.

##  S3 Additional real data results

Table S1: The detected communities by the proposed method. community |  category |  number
---|---|---
1 |  "Animal Shelters","Community Service/Non-Profit", "Greek", "Mediterranean", "Middle Eastern", "Pet Groomers", "Pet Services", "Pet Stores", "Pets", "Turkish", "Veterinarians" |  11
2 |  "Barbeque", "Cajun/Creole" "Southern", "Tex-Mex" |  4
3 |  "Bagels" "Bakeries" "Breakfast & Brunch" "Bubble Tea" "Burgers" "Cafes" "Cantonese" "Caribbean" "Caterers" "Chicken Wings" "Chinese" "Coffee & Tea" "Delis", "Desserts", "Dim Sum" "Diners" "Donuts" "Fast Food" "Food Delivery Services" "Gluten-Free" "Hot Dogs" "Ice Cream & Frozen Yogurt" "Italian" "Juice Bars & Smoothies" "Kosher" "Latin American" "Mexican" "Noodles" "Pizza" "Restaurants" "Salad" "Sandwiches" "Soup" "Tea Rooms" "Thai" "Vegan" "Vegetarian" "Vietnamese" |  38
4 |  "Custom Cakes" "French" "Pasta Shops" "Professional Services" "Seafood" "Steakhouses" "Taiwanese" |  7
5 |  "Beer" "Candy Stores" "Cheese Shops" "Chocolatiers & Shops" "Convenience Stores" "Drugstores" "Ethnic Food" "Farmers Market" "Florists" "Food" "Fruits & Veggies" "Grocery" "Health Markets" "Meat Shops" "Modern European" "Organic Stores" "Seafood Markets" "Specialty Food" "Wine & Spirits" |  19
6 |  "Adult Entertainment" "Bars" "Beer Bar" "Breweries" "Brewpubs" "British" "Cocktail Bars" "Comfort Food" "Dance Clubs" "Dive Bars" "Gay Bars" "German" "Irish" "Irish Pub" "Jazz & Blues" "Karaoke" "Lounges" "Music Venues" "Nightlife" "Pool Halls" "Pubs" "Spanish" "Sports Bars" "Tapas Bars" "Tapas/Small Plates" "Wine Bars" |  26
7 |  "Barbers" "Beauty & Spas" "Cosmetics & Beauty Supply" "Day Spas" "Doctors" "Fitness & Instruction", "Gyms", "Hair Removal" "Hair Salons" "Health & Medical" "Makeup Artists" "Massage" "Nail Salons", "Optometrists", "Skin Care" "Trainers", "Waxing", "Yoga" |  18
8 |  "Accessories" "Antiques" "Appliances" "Art Galleries" "Art Supplies" "Arts & Crafts" "Bike Rentals" "Bike Repair/Maintenance" "Bikes" "Books" "Bookstores" "Cards & Stationery" "Children’s Clothing" "Colleges & Universities" "Computers" "Department Stores" "Dry Cleaning & Laundry" "Education" "Electronics" "Fashion" "Flowers & Gifts", "Furniture Stores" "Gift Shops" "Hardware Stores" "Hobby Shops" "Home & Garden" "Home Decor" "Home Services" "Jewelry" "Kitchen & Bath" "Laundry Services" "Local Services" "Mags" "Men’s Clothing" "Music & DVDs" "Music & Video" "Nurseries & Gardening" "Outdoor Gear" "Real Estate" "Shoe Stores" "Shopping" "Shopping Centers" "Specialty Schools" "Sporting Goods" "Sports Wear" "Thrift Stores" "Toy Stores" "Used" "Vintage & Consignment" "Vinyl Records" "Women’s Clothing" |  51
9 |  "Active Life" "Amusement Parks" "Arcades" "Arts Entertainment" "Botanical Gardens" "Bowling" "Cinema" "Event Planning & Services" "Hotels" "Hotels & Travel" "Landmarks & Historical Buildings" "Local Flavor" "Museums" "Parks" "Party & Event Planning" "Public Services & Government" "Public Transportation" "Tours" "Transportation" "Venues & Event Spaces" |  21
10 |  "Asian Fusion", "Buffets", "Indian", "Japanese", "Korean", "Pakistani", "Sushi Bars" |  7
11 |  "Auto Parts & Supplies" "Auto Repair" "Automotive" "Gas Stations" "Tires" |  5
Table S2: The detected communities by Method 2. community |  category
---|---
1 |  "Animal Shelters" "Community Service/Non-Profit", "Greek" "Mediterranean", "Middle Eastern", "Pet Groomers", "Pet Services", "Pet Stores" "Pets" "Turkish", "Vegan" "Vegetarian" "Veterinarians"
2 |  "Barbeque", "Dim Sum", "Southern", "Tapas/Small Plates", "Tex-Mex"
3 |  "Bagels" "Bakeries" "Breakfast & Brunch" "Bubble Tea" "Burgers" "Cafes" "Cajun/Creole" "Cantonese" "Caribbean" "Caterers" "Chicken Wings" "Chinese" "Coffee & Tea" "Delis" "Desserts" "Diners" "Donuts" "Fast Food" "Food Delivery Services" "German" "Hot Dogs" "Ice Cream & Frozen Yogurt" "Italian" "Kosher" "Latin American" "Mexican" "Noodles""Pizza" "Restaurants" "Salad" "Sandwiches" "Seafood" "Soup" "Spanish" "Steakhouses" "Tea Rooms" "Thai"
4 |  "Asian Fusion" "Custom Cakes" "French" "Hobby Shops" "Japanese" "Korean" "Sushi Bars" "Taiwanese" "Vietnamese"
5 |  ""Beer" "Candy Stores" "Cheese Shops" "Chocolatiers & Shops" "Convenience Stores" "Drugstores" "Ethnic Food" "Farmers Market" "Florists" "Flowers & Gifts" "Food" Fruits & Veggies" "Gluten-Free" "Grocery" "Health Markets" "Juice Bars & Smoothies" "Meat Shops" "Organic Stores" "Pasta Shops" "Seafood Markets" "Specialty Food" "Wine & Spirits"
6 |  "Adult Entertainment" "Arts & Entertainment" "Bars" "Beer Bar" "Breweries" "Brewpubs" "British" "Cinema" "Cocktail Bars" "Comfort Food" "Dance Clubs" "Dive Bars" "Gay Bars" "Irish" "Irish Pub"n"Jazz & Blues" "Karaoke" "Lounges" "Modern European" "Music Venues" "Nightlife" "Performing Arts" "Pool Halls" "Pubs" "Sports Bars" "Tapas Bars" "Wine Bars"
7 |  "Barbers" "Beauty & Spas" "Cosmetics & Beauty Supply" "Day Spas" "Doctors" "Hair Removal" "Hair Salons" "Health & Medical" "Makeup Artists" "Massage" "Nail Salons" "Optometrists" "Skin Care" "Waxing"
8 |  "Accessories" Antiques" "Appliances" "Art Galleries" "Art Supplies" "Arts & Crafts" "Bike Rentals" "Bike Repair/Maintenance" "Bikes" "Books" "Bookstores" "Cards & Stationery" "Children’s Clothing" "Computers" "Department Stores" "Dry Cleaning & Laundry" "Electronics" "Fashion" "Furniture Stores" "Gift Shops" "Hardware Stores" "Home & Garden" "Home Decor" "Home Services" "Jewelry" "Kitchen & Bath" "Laundry Services" "Local Services" "Mags" "Men’s Clothing" "Music & DVDs" "Music & Video" "Nurseries & Gardening" "Outdoor Gear" "Professional Services" "Real Estate" "Shoe Stores" "Shopping" "Shopping Centers" "Specialty Schools" "Sporting Goods" "Sports Wear" "Thrift Stores" "Toy Stores" "Used" "Vintage & Consignment" "Vinyl Records" "Women’s Clothing"
9 |  "Active Life" "Amusement Parks" "Arcades" "Botanical Gardens" "Bowling" "Colleges & Universities" "Education" "Event Planning & Services" "Fitness & Instruction" "Gyms" "Hotels" "Hotels & Travel" "Landmarks & Historical Buildings" "Local Flavor" "Museums" "Parks" "Party & Event Planning" "Public Services & Government" "Public Transportation" "Tours" "Trainers" "Transportation" "Venues & Event Spaces" "Yoga"
10 |  "Buffets" "Indian" "Pakistani"
11 |  "Auto Parts & Supplies" "Auto Repair" "Automotive" "Gas Stations" "Tires"