Estimating mixed memberships in multi-layer networks

Multi-layer networks have emerged as a powerful tool for describing complex real-world systems. Unlike single-layer networks, these structures consist of multiple layers, where nodes represent entities and edges represent their interactions [30, 24, 5]. In this paper, we focus on multi-layer networks with the same nodes set in all layers and nodes solely interacting within their respective layers, excluding any cross-layer connections. Such networks are ubiquitous, spanning from social networks to transportation systems, biological networks, and trade networks. For instance, in social networks, individuals tend to form connections within different platforms (e.g., Facebook, Twitter, WeChat, LinkedIn), while cross-platform connections are typically not allowed. In transportation networks, layers might represent different modes of transportation (roads, railways, airways), with nodes corresponding to locations and edges indicating the availability of a specific mode between two locations [5]. The absence of cross-layer connections indicates the inability to directly switch modes without intermediate stops. Biological networks also exhibit rich multi-layer structures. In gene regulatory networks, layers could represent different types of gene interactions, with nodes representing genes and edges depicting their specific interactions [31, 4, 46]. The absence of cross-layer connections underscores the specialized nature of interactions within each layer and their vital role in governing cellular functions. The FAO-trade multi-layer network gathers trade relationships between countries across various products sourced from the Food and Agriculture Organization of the United Nations [10]. Figure 1 illustrates the networks of the first three products within this dataset.

Community detection in multi-layer networks is a crucial analytical tool, revealing latent structures within networks [23, 18]. A community (also known as a group, cluster, or block) typically comprises nodes more densely interconnected than those outside [35, 36, 11, 12, 19]. In social networks, individuals often form distinct communities based on interests, occupations, or locations. For instance, on Facebook, users might belong to hobby-based communities like photography or hiking, identifiable by dense interaction patterns among members. In transportation networks, communities emerge due to geographical proximity or functional similarity. For instance, cities might form communities linked by tightly integrated railways. In practical applications, a node can simultaneously belong to multiple communities. For example, in social networks, an individual may be a member of several different social groups. Similarly, in transportation networks, a specific location can function as a hub for multiple modes of transportation, bridging diverse communities. Likewise, in biological networks, a gene can belong to several overlapping communities, participating in a wide range of processes.

In the last few years, community detection in multi-layer networks, where each node belongs exclusively to a single community, has attracted considerable attention. For instance, several studies have been proposed under the multi-layer stochastic block model (MLSBM), a model that assumes the network for each layer can be modeled by the well-known stochastic block model (SBM) [17]. MLSBM also assumes that the community information is common to all layers. Under the framework within MLSBM, [16] studied the consistent community detection for maximum likelihood estimate and a spectral method designed by running the K-means algorithm on a few eigenvectors of the sum of adjacency matrices when only the number of layers grows. [37] studied the consistency of co-regularized spectral clustering, orthogonal linked matrix factorization, and the spectral method in [16] as both the number of nodes and the number of layers grow within the MLSBM context. Moreover, [26] established consistency results for a least squares estimation of community memberships within the MLSBM framework. Recently, [27] studied the consistency of a bias-adjusted (i.e., debiased) spectral clustering method using a novel debiased sum of squared adjacency matrices under MLSBM. Their numerical experiments demonstrated that this approach significantly outperforms the spectral method proposed by [16]. Also see [39, 21, 9, 43] for other recent works exploring variants of MLSBM.

However, the MLSBM model has one significant limitation, being tailored exclusively for multi-layer networks with non-overlapping communities. The mixed membership stochastic block (MMSB) model [1] is a popular statistical model for capturing mixed community memberships in single-layer networks. Under MMSB, some methods have been developed to estimate mixed memberships for single-layer networks generated from MMSB, including variational expectation-maximization [1, 15], nonnegative matrix factorization [40, 45], tensor decomposition [2], and spectral clustering [29, 20]. In this paper, we consider the problem of estimating nodes’ common community memberships in multi-layer networks generated from the multi-layer mixed membership stochastic block (MLMMSB) model, a multi-layer version of MMSB. The main contributions of this paper are summarized as follows:

• 1.

  We introduce three spectral methods for estimating mixed memberships in multi-layer networks generated from MLMMSB. These methods employ vertex-hunting algorithms on a few selected eigenvectors of three aggregate matrices: the sum of adjacency matrices, the debiased sum of squared adjacency matrices, and the sum of squared adjacency matrices.

• 2.

  We establish per-node error rates for these three methods under mild conditions on network sparsity, demonstrating their consistent mixed membership estimation as the number of nodes and/or layers increases within the MLMMSB framework. Our theoretical analysis reveals that the method utilizing the debiased sum of squared adjacency matrices consistently outperforms the method using the sum of squared adjacency matrices in terms of error rate. Additionally, both methods generally exhibit lower error rates than the method based on the sum of adjacency matrices. This underscores the advantage of debiased spectral clustering in mixed membership community detection for multi-layer networks. To the best of our knowledge, this is the first work to estimate mixed memberships using the above three aggregate matrices and the first work to establish consistency results within the MLMMSB framework.

• 3.

  To assess the quality of mixed membership community detection in multi-layer networks, we introduce two novel modularity metrics: fuzzy sum modularity and fuzzy mean modularity. The first is derived from computing the fuzzy modularity of the sum of adjacency matrices, while the second is obtained by averaging the fuzzy modularity of adjacency matrices across individual layers. To the best of our knowledge, our two modularity metrics are the first to measure the quality of mixed membership community detection in multi-layer networks.

• 4.

  We conduct extensive simulations to validate our theoretical findings and demonstrate the practical effectiveness of our methods and metrics through real-world multi-layer network applications.

The remainder of this paper is structured as follows: Section 2 introduces the model, followed by Section 3 outlining the methods developed. Section 4 presents the theoretical results. Subsequently, Section 5 includes experimental results on computer-generated multi-layer networks, while Section 6 focuses on real-world multi-layer networks. Lastly, Section 7 concludes the paper and technical proofs are provided in A.

Notation. We employ the notation $[m]$ to denote the set $\{1,2,\ldots,m\}$. Furthermore, $I_{m\times m}$ stands for the $m$-by-$m$ identity matrix. For a vector $x$, $\|x\|_{q}$ represents its $l_{q}$ norm. When considering a matrix $M$ and an index set $s$ that is a subset of $[m]$, $M(s,:)$ refers to the sub-matrix of $M$ comprising the rows indexed by $s$. Additionally, $M^{\prime}$ denotes the transpose of $M$, $\|M\|_{F}$ is its Frobenius norm, $\|M\|_{\infty}$ represents the maximum absolute row sum, $\|M\|_{2\rightarrow\infty}$ signifies the maximum row-wise $l_{2}$ norm, $\mathrm{rank}(M)$ gives its rank, and $\lambda_{k}(M)$ stands for the $k$-th largest eigenvalue of $M$ in magnitude. The notation $\mathbb{E}[\cdot]$ is used to denote expectation, while $\mathbb{P}(\cdot)$ represents probability. Finally, $e_{i}$ is defined as the indicator vector with a 1 in the $i$-th position and zeros elsewhere.

Throughout, we consider undirected and unweighted multi-layer networks with $n$ common nodes and $L$ layers. For the $l$-th layer, let the $n\times n$ symmetric matrix $A_{l}$ be its adjacency matrix such that $A_{l}(i,j)=A_{l}(j,i)=1$ if there is an edge connecting node $i$ and node $j$ and $A_{l}(i,j)=A_{l}(j,i)=0$ otherwise, for all $i\in[n],j\in[n],l\in[L]$, i.e., $A_{l}=A^{\prime}_{l}\in\{0,1\}^{n\times n}$. Additionally, we allow for the possibility of self-edges (loops) in this paper. We assume that the multi-layer network consists of $K$ common communities

Throughout, we assume that the number of communities $K$ is known in this paper. The challenging task of theoretically estimating $K$ in multi-layer networks exceeds the scope of this paper. Let $\Pi\in[0,1]^{n\times K}$ be the membership matrix such that $\Pi(i,k)$ represents the “weight” that node $i$ belongs to the $k$-th community $\mathcal{C}_{k}$ for $i\in[n],k\in[K]$. Suppose that $\Pi$ satisfies the following condition

We call node $i$ a “pure” node if one entry of $\Pi(i,:)$ is 1 while the other $(K-1)$ elements are zeros and a “mixed” node otherwise. Assume that

Based on Equation (3), let $\mathcal{I}$ be the index set of pure nodes such that $\mathcal{I}=\{p_{1},p_{2},\ldots,p_{K}\}$ with $p_{k}\in\{1,2,\ldots,n\}$ being an arbitrary pure node in the $k$-th community $\mathcal{C}_{k}$ for $k\in[K]$. Similar to [29], without loss of generality, we reorder the nodes such that $\Pi(\mathcal{I},:)=I_{K\times K}$.

Suppose that all the layers share a common mixed membership matrix $\Pi$ but with possibly different edge probabilities. In particular, we work with a multi-layer version of the popular mixed membership stochastic block model (MMSB) [1]. To be precise, we define the multi-layer MMSB below:

Since $B_{l}$ can vary for different $l$, $A_{l}$ generated by Equation (4) may have different expectation $\Omega_{l}\equiv\mathbb{E}[A_{l}]$ for $l\in[L]$. Notably, when $L=1$, the MLMMSB model degenerates to the popular MMSB model. When all nodes are pure, MLMMSB reduces to the MLSBM studied in previous works [16, 37, 26, 27]. Define $\mathcal{B}=\{B_{1},B_{2},\ldots,B_{L}\}$. Equation (4) says that MLMMSB is parameterized by $\Pi,\rho$, and $\mathcal{B}$. For brevity, we denote MLMMSB defined by Equation (4) as $\mathrm{MLMMSB}(\Pi,\rho,\mathcal{B})$. Since $\mathbb{P}(A_{l}(i,j)=1)=\Omega_{l}(i,j)=\rho\Pi(i,:)B_{l}\Pi^{\prime}(j,:)\leq\rho$, we see that decreasing the value of $\rho$ results in a sparser multi-layer network, i.e., $\rho$ controls the overall sparsity of the multi-layer network. For this reason, we call $\rho$ the sparsity parameter in this paper. We allow the sparsity parameter $\rho$ to go to zero by either increasing the number of nodes $n$, the number of layers $L$, or both simultaneously. We will study the impact of $\rho$ on the performance of the proposed methods by considering $\rho$ in our theoretical analysis. The primal goal of community detection within the MLMMSB framework is to accurately estimate the common mixed membership matrix $\Pi$ from the observed $L$ adjacency matrices $\{A_{l}\}^{L}_{l=1}$. This estimation task is crucial for understanding the underlying community structure in multi-layer networks.

In this section, we propose three spectral methods designed to estimate the mixed membership matrix $\Pi$ within the MLMMSB model for multi-layer networks in which nodes may belong to multiple communities with different weights. We recall that $\Omega_{l}$ is the expectation of the $l$-th observed adjacency matrix $A_{l}$ for $l\in[L]$ and $\Pi$ contains the common community memberships for all nodes. To provide intuitions of the designs of our methods, we consider the oracle case where $\{\Omega_{l}\}^{L}_{l=1}$ are directly observed. First, we define two distinct aggregate matrices formed by $\{\Omega_{l}\}^{L}_{l=1}$: $\Omega_{\mathrm{sum}}\equiv\sum_{l\in[L]}\Omega_{l}$ and $\tilde{S}_{\mathrm{sum}}\equiv\sum_{l\in[L]}\Omega^{2}_{l}$, where the former represents the sum of all expectation adjacency matrices, the latter is the sum of all squared expectation adjacency matrices, and both matrices contain the information about the mixed memberships of nodes in the multi-layer network since $\Omega_{\mathrm{sum}}=\rho\Pi(\sum_{l\in[L]}B_{l})\Pi^{\prime}$ and $\tilde{S}_{\mathrm{sum}}=\rho^{2}\Pi(\sum_{l\in[L]}B_{l}\Pi^{\prime}\Pi B_{l})\Pi^{\prime}$. Since $\mathrm{rank}(\Pi)=K$, it is easy to see that $\mathrm{rank}(\Omega_{\mathrm{sum}})=K$ if $\mathrm{rank}(\sum_{l\in[L]}B_{l})=K$ and $\mathrm{rank}(\tilde{S}_{\mathrm{sum}})=K$ if $\mathrm{rank}(\sum_{l\in[L]}B^{2}_{l})=K$. Assuming that the number of communities $K$ is much smaller than the number of nodes $n$, we observe that $\Omega_{\mathrm{sum}}$ and $\tilde{S}_{\mathrm{sum}}$ possess low-dimensional structure. This low-rank property is crucial for the development of our spectral methods, as it allows us to efficiently extract meaningful information from the high-dimensional data. Building on these insights, we present the following lemma that characterizes the geometries of the two aggregate matrices $\Omega_{\mathrm{sum}}$ and $\tilde{S}_{\mathrm{sum}}$. This lemma forms the theoretical foundation for our methods, enabling us to develop accurate and computationally efficient algorithms for estimating the mixed membership matrix $\Pi$ within the MLMMSB model.

Recall that $\Pi$ satisfies Equations (2) and (3), by $U=\Pi U(\mathcal{I},:)$ according to Lemma 1, the rows of $U$ forms a $K$-simplex in $\mathbb{R}^{K}$, referred to as the Ideal Simplex of $U$ ($IS_{U}$ for brevity). Notably, the $K$ rows of $U(\mathcal{I},:)$ serve as its vertices. Given that $\Pi$ satisfies Equations (2) and (3), and $U=\Pi U(\mathcal{I},:)$, it follows that $U(i,:)=\Pi(i,:)U(\mathcal{I},:)$ for $i\in[n]$. This implies that $U(i,:)$ is a convex linear combination of the $K$ vertices of $IS_{U}$ with weights determined by $\Pi(i,:)$. Consequently, a pure row of $U$ (call $U(i,:)$ pure row if node $i$ is a pure node and mixed row otherwise) lies on one of the $K$ vertices of the simplex, while a mixed row occupies an interior position within the simplex. Such simplex structure is also observed in mixed membership estimation for community detection in single-layer networks [29, 20, 41] and in topic matrix estimation for topic modeling [22]. Notably, the mixed membership matrix $\Pi$ can be precisely recovered by setting $\Pi=UU^{-1}(\mathcal{I},:)$ provided that the corner matrix $U(\mathcal{I},:)$ is obtainable. Thanks to the simplex structure in $U=\Pi U(\mathcal{I},:)$, the vertices of the simplex can be found by a vertex-hunting technique such as the successive projection (SP) algorithm [3, 13, 14]. Applying SP to all rows of $U$ with $K$ clusters enables the exact recovery of the $K\times K$ corner matrix $U(\mathcal{I},:)$. Details of the SP algorithm are provided in Algorithm 1 [14]. SP can efficiently find the $K$ distinct pure rows of $U$. Similarly, $V=\Pi V(\mathcal{I},:)$ also exhibits a simplex structure, leading to the recovery of $\Pi$ via $\Pi=VV^{-1}(\mathcal{I},:)$ by Lemma 1. Consequently, applying the SP algorithm to $V$ with $K$ clusters also yields an exact reconstruction of the mixed membership matrix $\Pi$.

For real-world multi-layer networks, the $L$ adjacency matrices $A_{1},A_{2},\ldots,A_{L}$ are observed instead of their expectations. For our first method, set $A_{\mathrm{sum}}=\sum_{l\in[L]}A_{l}$. We have $\mathbb{E}[A_{\mathrm{sum}}]=\Omega_{\mathrm{sum}}$ under MLMMSB. Subsequently, let $\hat{U}\hat{\Sigma}\hat{U}^{\prime}$ be the top $K$ eigen-decomposition of $A_{\mathrm{sum}}$ where $\hat{U}$ is an orthogonal matrix with $\hat{U}^{\prime}\hat{U}=I_{K\times K}$ and $\hat{\Sigma}$ is a diagonal matrix with its $k$-th diagonal element representing the $k$-th largest eigenvalue of $A_{\mathrm{sum}}$ in magnitude for $k\in[K]$. Given that the expectation of $A_{\mathrm{sum}}$ is $\Omega_{\mathrm{sum}}$, it follows that $\hat{U}$ can be interpreted as a slightly perturbed version of $U$. To proceed, we apply the SP algorithm to all rows of $\hat{U}$ with $K$ clusters, resulting in an estimated index set of pure nodes, denoted as $\hat{\mathcal{I}}$. We infer that $\hat{U}(\hat{\mathcal{I}},:)$ should closely approximate the corner matrix $U(\mathcal{I},:)$. Finally, we estimate the mixed membership matrix $\Pi$ by computing $\hat{U}\hat{U}^{-1}(\hat{\mathcal{I}},:)$. Our first algorithm, called “Successive projection on the sum of adjacency matrices” (SPSum, Algorithm 1) summarises the above analysis. It efficiently estimates the mixed membership matrix $\Pi$ by executing the SP algorithm on the top $K$ eigenvectors of $A_{\mathrm{sum}}$.

In our second methodological approach, we define the aggregate matrix $S_{\mathrm{sum}}$ as the sum of the squared adjacency matrices, adjusted for bias, specifically as $S_{\mathrm{sum}}=\sum_{l\in[L]}(A_{l}^{2}-D_{l})$. Here, $D_{l}$ represents a diagonal matrix, with its $i$-th diagonal entry being the degree of node $i$ in layer $l$, i.e., $D_{l}(i,i)=\sum_{j\in[n]}A_{l}(i,j)$ for $i\in[n]$ and $l\in[L]$. The matrix $S_{\mathrm{sum}}$ was first introduced in [27] within the context of MLSBM, serving as a debiased estimate of the sum of squared adjacency matrices. This debiasing is necessary as $\sum_{l\in[L]}A_{l}^{2}$ alone provides a biased approximation of $\sum_{l\in[L]}\Omega_{l}^{2}$. By subtracting the diagonal matrix $D_{l}$ from each squared adjacency matrix, we can effectively remove this bias, rendering $S_{\mathrm{sum}}$ a reliable estimator of $\tilde{S}_{\mathrm{sum}}$ as demonstrated in [27]. Subsequently, let $\hat{V}\hat{\Lambda}\hat{V}^{\prime}$ be the top $K$ eigen-decomposition of $S_{\mathrm{sum}}$. Given that $\hat{V}$ provides a close approximation of $V$, applying the SP algorithm to all rows of $\hat{V}$ yields a reliable estimate of the corner matrix $V(\mathcal{I},:)$. In summary, our second spectral method, centered on $S_{\mathrm{sum}}$, is outlined in Algorithm 2. We refer to this method as “Successive projection on the debiased sum of squared adjacency matrices” (SPDSoS for brevity).

The third method utilizes the summation of squared adjacency matrices, expressed as $\sum_{l\in[L]}A^{2}_{l}$, as a substitute for $S_{\mathrm{sum}}$ in Algorithm 2. Notably, this approach excludes the crucial bias-removal step inherent in $S_{\mathrm{sum}}$. We call this method “Successive projection on the sum of squared adjacency matrices”, abbreviated as SPSoS. In the subsequent section, we will demonstrate that SPDSoS consistently outperforms SPSoS, and both methods are generally theoretically superior to SPSum.

In this section, we demonstrate the consistency of our methods by presenting their theoretical upper bounds for per-node error rates as the number of nodes $n$ and/or the number of layers $L$ increases within the context of the MLMMSB model. Assumption 1 provides a prerequisite lower bound for the sparsity parameter $\rho$ to ensure the theoretical validity of SPSum.

In Assumption 1, the choice of $\tau$ plays a crucial role in determining the level of sparsity regime. Since $\sum_{l\in[L]}A_{l}(i,j)\leq L$, $\tau$ has an upper bound $L$. However, to consider an even sparser scenario, we introduce a constant positive value $\beta$ that is strictly less than $L$ and remains fixed regardless of the growth of $n$ or $L$. This ensures that the aggregate number of edges connecting any two nodes $i$ and $j$ across all layers remains below $\beta$ even in the limit as $n$ or $L$ approaches infinity. Under these conditions, Assumption 1 is revised to reflect a more stringent requirement: $\rho\geq\frac{\beta^{2}\mathrm{log}(n+L)}{nL}$. This revised assumption characterizes the most challenging scenario for community detection. To establish the theoretical bounds for SPSum, we need the following requirement on $\mathcal{B}$.

Assumption 2 is mild as $|\lambda_{K}(\sum_{l\in[L]}B_{l})|$ represents the smallest singular value resulting from the summation of $L$ matrices, and it is reasonable to presume that $|\lambda_{K}(\sum_{l\in[L]}B_{l})|$ is on the order of $O(L)$. To simplify our theoretical analysis, we also introduce the following condition:

Condition 1 is mild as $K=O(1)$ implies that the number of communities remains constant, while $\lambda_{K}(\Pi^{\prime}\Pi)=O(\frac{n}{K})$ ensures the balance of community sizes, where the size of the $k$-th community $\mathcal{C}_{k}$ is defined as $\sum_{i\in[n]}\Pi(i,k)$ for $k\in[K]$. Our main result for SPSum offers an upper bound on its per-node error rate in terms of $\rho,n$, and $L$ under the MLMMSB framework.

According to Theorem 1, it becomes evident that as $n$ (and/or $L$) approaches infinity, SPSum’s error rate diminishes to zero, highlighting the advantage of employing multiple layers in community detection and SPSum’s consistent community detection. Additionally, Theorem 1 says that to attain a sufficiently low error rate for SPSum, the sparsity parameter $\rho$ must decrease at a rate slower than $\frac{\mathrm{log}(n+L)}{nL}$ (i.e., $\rho$ must be significantly greater than $\frac{\mathrm{log}(n+L)}{nL}$). This finding aligns with the sparsity requirement stated in Assumption 1 for SPSum. It is noteworthy that the theoretical upper bound for SPSum’s per-node error rate, given as $O(\sqrt{\frac{\mathrm{log}(n+L)}{\rho nL}})$ in Theorem 1, is independent of how we select the value of $\tau$ in Assumption 1. This underscores the generality of our theoretical findings.

The assumption 3 stated below establishes the lower bound of the sparsity parameter $\rho$ that ensures the consistency of both SPDSoS and SPSoS.

Analogous to Assumption 1, when considering an extremely sparse scenario where $\tilde{\tau}$ does not exceed a positive constant $\tilde{\beta}$, Assumption 3 dictates that $\rho$ must satisfy $\rho\geq\frac{1}{n}\sqrt{\frac{\tilde{\beta}^{2}\mathrm{log}(n+L)}{L}}$. This requirement aligns with the sparsity requirement in Theorem 1 of [27]. The subsequent assumption serves a similar purpose to Assumption 2.

Theorems 2 and 3 are the main results of SPDSoS and SPSoS, respectively.

According to Theorems 2 and 3, it is evident that both SPDSoS and SPSoS enjoy consistent mixed membership estimation. This consistency arises from the fact that their per-node error rates tend towards zero as the number of nodes $n$ (and/or the number of layers $L$) approaches infinity. Furthermore, similar to Theorem 1, the theoretical bounds outlined in Theorems 2 and 3 are independent of the selection of $\tilde{\tau}$ in Assumption 3.

By the proof of Theorem 3, we know that SPSoS’s theoretical upper bound of per-node error rate is the summation of SPDSoS’s theoretical upper bound of per-node error rate and $O(\frac{1}{\rho n})$. Therefore, SPDSoS’s error rate is always smaller than that of SPSoS and this indicates the benefit of the bias-removal step in $S_{\mathrm{sum}}$. Furthermore, to compare the theoretical performances between SPSum and SPSoS, we consider the following two simple cases:

• 1.

  When SPSoS’s error rate is at the order $\frac{1}{\rho n}$, if SPSum significantly outperforms SPSoS, we have $\sqrt{\frac{\mathrm{log}(n+L)}{\rho nL}}\ll\frac{1}{\rho n}\Leftrightarrow L\gg\rho n\mathrm{log}(n+L)$.

• 2.

  When SPSoS’s error rate is at the order $\sqrt{\frac{\mathrm{log}(n+L)}{\rho^{2}n^{2}L}}$, if SPSum significantly outperforms SPSoS, we have $\sqrt{\frac{\mathrm{log}(n+L)}{\rho nL}}\ll\sqrt{\frac{\mathrm{log}(n+L)}{\rho^{2}n^{2}L}}\Leftrightarrow\rho n\ll 1\Leftrightarrow\rho\ll\frac{1}{n}$. Since $\sqrt{\frac{\mathrm{log}(n+L)}{\rho nL}}\ll 1$, we have $\rho\gg\frac{\mathrm{log}(n+L)}{nL}$. Combine $\rho\gg\frac{\mathrm{log}(n+L)}{nL}$ with $\rho\ll\frac{1}{n}$, we have $\frac{\mathrm{log}(n+L)}{nL}\ll\frac{1}{n}\Leftrightarrow L\gg\mathrm{log}(n+L)$.

The preceding points indicate that SPSum’s superior performance over SPSoS is limited to scenarios where the number of layers $L$ is exceptionally large ($L\gg\rho n\mathrm{log}(n+L)$ or $L\gg\mathrm{log}(n+L)$). However, this requirement is unrealistic for the majority of real-world multi-layer networks. Given that SPDSoS consistently outperforms SPSoS, it follows that SPDSoS also prevails over SPSum in most scenarios. Additionally, even in the rare scenarios where SPSum significantly surpasses SPDSoS for a significantly large $L$, Theorem 2 assures that SPDSoS’s error rate is negligible. Therefore, among these three methods, we arrive at the following conclusions: (a) SPDSoS consistently outperforms SPSoS; (b) SPSum’s significant superiority over SPDSoS and SPSoS necessitates an unreasonably high number of layers, $L$, which is impractical for most real-world multi-layer networks. Consequently, we can confidently assert that SPDSoS and SPSoS virtually always surpass SPSum.

In this section, we evaluate the performance of our proposed methods through the utilization of computer-generated multi-layer networks. For these simulated networks, we possess knowledge of the ground-truth mixed membership matrix $\Pi$. To quantify the performance of each method, we employ two metrics: the Hamming error and the Relative error. The Hamming error is defined as $\mathrm{Hamming~{}error}=\mathrm{min}_{\mathcal{P}\in\mathcal{S}}\frac{\|\hat{\Pi}-\Pi\mathcal{P}\|_{1}}{n}$, while the relative error is given by $\mathrm{Relative~{}error}=\mathrm{min}_{\mathcal{P}\in\mathcal{S}}\frac{\|\hat{\Pi}-\Pi\mathcal{P}\|_{F}}{\|\Pi\|_{F}}$. Here, $\mathcal{S}$ represents the set of all $K$-by-$K$ permutation matrices, accounting for potential label permutations. For each parameter setting in our simulation examples, we report the average of each metric for each proposed approach across 100 independent repetitions.

For all simulations conducted below, we set $K=3$ and let $n_{0}$ be the number of pure nodes within each community. For mixed nodes, say node $i$, we formulate its mixed membership vector $\Pi(i,:)$ in the following manner. Initially, we generate two random values $r^{(1)}_{i}=\mathrm{rand}(1)$ and $r^{(2)}_{i}=\mathrm{rand}(1)$, where $\mathrm{rand}(1)$ represents a random value drawn from a Uniform distribution over the interval $[0,1]$. Subsequently, we define the mixed membership vector as $\Pi(i,:)=(\frac{r^{(1)}_{i}}{2},\frac{r^{(2)}_{i}}{2},1-\frac{r^{(1)}_{i}}{2}-\frac{r^{(2)}_{i}}{2})$. Regarding the connectivity matrices, for the $l$-th matrix $B_{l}$, we assign $B_{l}(k,\tilde{k})=B_{l}(\tilde{k},k)=\mathrm{rand}(1)$ for all $1\leq k\leq\tilde{k}\leq K$ and $l\in[L]$. Finally, the number of nodes $n$, the sparsity parameter $\rho$, the number of layers $L$, and the number of pure nodes $n_{0}$ within each community are independently set for each simulation.

Experiment 1: Changing $\rho$. Fix $(n,L,n_{0})=(500,100,100)$ and let $\rho$ range in$\{0.02,0.04,\ldots,0.2\}$. The multi-layer network becomes denser as $\rho$ increases. The results shown in panel (a) and panel (b) of Figure 2 demonstrate that SPDSoS performs slightly better than SPSoS, and both methods significantly outperform SPSum. Meanwhile, the error rates of SPDSoS and SPSoS decrease rapidly as the sparse parameter $\rho$ increases while SPSum’s error rates decrease slowly.

Experiment 2: Changing $L$. Fix $(n,\rho,n_{0})=(500,0.1,100)$ and let $L$ range in$\{10,20,\ldots,100\}$. More layers are observed as $L$ increases. The results are presented in panel (c) and panel (d) of Figure 2. It is evident that SPDSoS and SPSoS exhibit comparable performance, both significantly surpassing SPSum. Furthermore, as $L$ increases, SPDSoS and SPSoS demonstrate improved performance, whereas SPSum’s performance remains relatively unchanged throughout this experiment.

Experiment 3: Changing $n$. Fix $(L,\rho)=(40,0.1)$, let $n$ range in $\{200,400,\ldots,2000\}$, and let $n_{0}=\frac{n}{4}$ for each choice of $n$. Panel (e) and panel (f) of Figure 2 display the results, which demonstrate that SPDSoS and SPSoS enjoy similar performances and that both methods perform better than SPSum. Meanwhile, we also observe that the error rates of all methods decrease as $n$ increases here.

Experiment 4: Changing $n_{0}$. Fix $(n,L,\rho)=(600,50,0.1)$ and let $n_{0}$ range in $\{20,40,\ldots,200\}$. The number of pure nodes increases as $n_{0}$ grows. The results are displayed in the final two panels of Figure 2. Our observations indicate that SPDSoS and SPSoS exhibit notably superior performance compared to SPSum. Furthermore, both SPDSoS and SPSoS exhibit improved performance in scenarios with a higher number of pure nodes, whereas SPSum demonstrates weaker performance in this experiment.

In this section, we demonstrate the application of our methods to multi-layer networks in the real world. The application of a mixed membership estimation algorithm to such networks consistently yields an estimated mixed membership matrix, denoted as $\hat{\Pi}$. Notably, $\hat{\Pi}$ may vary depending on the specific algorithm used. Consequently, accurately assessing the quality of the estimated mixed membership community partition becomes a crucial problem. To address this challenge, we introduce two modularity metrics in this paper, designed to quantitatively evaluate the quality of mixed membership community detection in real-world multi-layer networks.

Recall that $A_{\text{sum}}=\sum_{l\in[L]}A_{l}$ represents the summation of all adjacency matrices. This summation effectively quantifies the weight between nodes. Given that nodes sharing similar community memberships tend to exhibit stronger connections than those with different memberships, $A_{\text{sum}}$ can be interpreted as the weighted adjacency matrix of an assortative network [33, 34]. Our fuzzy sum modularity, $Q_{\text{fsum}}$, is defined as follows:

where fsum stands for “fuzzy sum”, $d_{\text{sum}}(i)=\sum_{j\in[n]}A_{\text{sum}}(i,j)$ for $i\in[n]$, and $m_{\text{sum}}=\sum_{i\in[n]}d_{\text{sum}}(i)$. Notably, when $L=1$, our fuzzy sum modularity $Q_{\text{fsum}}$ simplifies to the fuzzy modularity introduced in Equation (14) of [32]. Furthermore, when $L=1$ and all nodes are pure, our modularity metric reduces to the classical Newman-Girvan modularity [36].

Our second modularity metric is the fuzzy mean modularity, defined as the average of fuzzy modularities across all layers. Here’s how it is formulated:

where fmean stands for “fuzzy mean” and $m_{l}=\sum_{i\in[n]}D_{l}(i,i)$ represents the sum of diagonal elements in the degree matrix $D_{l}$ for layer $l$ for $l\in[L]$. When all nodes are pure, our fuzzy mean modularity $Q_{\text{fmean}}$ simplifies to the multi-normalized average modularity introduced in Equation (3.1) of [38]. When there is only a single layer ($L=1$), $Q_{\text{fmean}}$ reduces to the fuzzy modularity described in [32]. Furthermore, if both conditions hold ($L=1$ and all nodes are pure), $Q_{\text{fmean}}$ degenerates to the classical Newman-Girvan modularity [36].

Analogous to the Newman-Girvan modularity [36], the fuzzy modularity [32], and the multi-normalized average modularity [38], a higher value of the fuzzy sum modularity $Q_{\text{fsum}}$ indicates a better community partition. Consequently, we consistently favor a larger value of $Q_{\text{fsum}}$. Similarly, a larger $Q_{\text{fmean}}$ also indicates a better community partition. To the best of our knowledge, our fuzzy sum modularity $Q_{\text{fsum}}$ and fuzzy mean modularity $Q_{\text{fmean}}$ are the first metrics designed to evaluate the quality of mixed membership community detection in multi-layer networks.

For real-world multi-layer networks, where the number of communities $K$ is usually unknown, we adopt the strategy introduced in [32] to estimate $K$. Specifically, we determine $K$ by selecting the one that maximizes the fuzzy sum modularity $Q_{\text{fsum}}$ (or the fuzzy mean modularity $Q_{\text{fmean}}$). This strategy ensures that we obtain an optimal community partition based on the chosen modularity metric.

In this paper, we consider the following real-world multi-layer networks, which can be accessed at https://manliodedomenico.com/data.php:

• 1.

  Lazega Law Firm: This data is a multi-layer social network with $n=71$ nodes and $L=3$ layers [25, 42]. For this data, nodes represent company partners and layers denote different cooperative relationships (Advice, Friendship, and Co-work) among partners.

• 2.

  C.Elegans: This data is from biology and it collects the connection of Caenorhabditis elegans [7]. It has $n=279$ nodes and $L=3$ layers with nodes denoting Caenorhabditis elegans and layers being different synaptic junctions (Electric, Chemical Monadic, and Chemical Polyadic).

• 3.

  CS-Aarhus: This data is a multi-layer social network with $n=61$ nodes and $L=5$ layers, where nodes represent employees of the Computer Science department at Aarhus and layers denote different relationships (Lunch, Facebook, Coauthor, Leisure, Work) [28].

• 4.

  FAO-trade: This data collects different types of trade relationships among countries from FAO (Food and Agriculture Organization of the United Nations) [10]. It has $n=214$ nodes and $L=364$ layers, where nodes represent countries, layers denote products, and edges denote trade relationships among countries.

The estimated number of communities and the corresponding modularity of each method for the real data analyzed in this paper are comprehensively presented in Table 1 and Table 2. After analyzing the results in these tables, we arrive at the following conclusions:

• 1.

  $Q_{\text{fsum}}$ and $Q_{\text{fmean}}$ consistently demonstrate similar values for each method, indicating a high degree of agreement. Specifically, SPSum achieves the highest modularity score using both metrics across multiple datasets, including Lazega Law Firm, CS-Aarhus, and FAO-trade. It is worth mentioning that this observation (SPSum surpassing the other two methods) contrasts with both theoretical and numerical discoveries. A plausible explanation for this discrepancy could be that $Q_{\text{fsum}}$ and $Q_{\text{fmean}}$ are computed directly from $\{A_{l}\}^{L}_{l=1}$ rather than $\{A^{2}_{l}\}^{L}_{l=1}$. Similarly, for SPSoS, its $Q_{\text{fsum}}$ and $Q_{\text{fmean}}$ scores for CS-Aarhus rank second among all methods, with scores that are closely aligned. This consistency across measures suggests that both $Q_{\text{fsum}}$ and $Q_{\text{fmean}}$ effectively capture similar aspects of community structure, providing reliable and comparable assessments of different methods.

• 2.

  While SPSum may not perform as well in numerical studies, it nevertheless exhibits superior performance in terms of modularity for real-world datasets compared to SPDSoS and SPSoS. The only exception is the C.Elegans network, where SPSum’s modularity score is slightly lower than SPDSoS’s.

• 3.

  The CS-Aarhus network exhibits a more distinct community structure compared to the other three real multi-layer networks, as evidenced by the higher modularity scores obtained by our methods for CS-Aarhus. Conversely, the FAO-trade network possesses the least discernible community structure among all real multi-layer networks, as the modularity scores achieved by all proposed methods for FAO-trade are lower than those of the other three networks.

• 4.

  The results presented in Table 1 and Table 2 indicate that the optimal number of communities for the four real-world multi-layer networks, namely Lazega Law Firm, CS-Aarhus, CS-Aarhus, and FAO-trade, are 3, 2, 5, and 2, respectively. This determination is based on selecting the value of $K$ that yields the highest modularity score across all methods for each dataset.

To simplify our analysis and consider the fact that SPSum generates modularity scores that are either higher or comparable to those of SPDSoS and SPSoS for the four real-world multi-layer networks under consideration, we focus our further analysis on the SPSum method. Let $\hat{\Pi}$ represent the estimated mixed membership matrix obtained by SPSum for a given real-world multi-layer network. For a node $i$ within the network, we define its estimated home base community as the one that corresponds to the maximum value in the $i$-th row of $\hat{\Pi}$, denoted as $\tilde{k}=\mathrm{argmax}_{k\in[K]}\hat{\Pi}(i,k)$. Furthermore, we categorize nodes based on their membership distribution: a node is considered highly mixed if $\mathrm{max}_{k\in[K]}\hat{\Pi}(i,k)\leq 0.6$, highly pure if $\mathrm{max}_{k\in[K]}\hat{\Pi}(i,k)\geq 0.9$, and neutral otherwise. We introduce two additional metrics, $\varsigma_{mixed}$ and $\varsigma_{pure}$, which represent the proportions of highly mixed and highly pure nodes, respectively. Additionally, we define the balanced parameter $\upsilon$ as the ratio of the minimum to the maximum $l_{1}$ norm of the columns of $\hat{\Pi}$, i.e., $\upsilon=\frac{\mathrm{min}_{k\in[K]}\|\hat{\Pi}(:,k)\|_{1}}{\mathrm{max}_{k\in[K]}\|\hat{\Pi}(:,k)\|_{1}}$. A higher value of $\upsilon$ indicates a more balanced multi-layer network. Table 3 presents the values of these three indices obtained by applying SPSum to the real-world multi-layer networks studied in this paper. Based on the results in Table 3, we draw the following conclusions:

• 1.

  Across all networks, the number of highly pure nodes significantly exceeds the number of highly mixed nodes, indicating that most nodes are strongly associated with a single community, while only a few exhibit mixed membership.

• 2.

  In the Lazega Law Firm network, approximately 10 nodes are highly mixed, 35 nodes are highly pure, and 26 are neutral. Meanwhile, this network exhibits the highest proportion of highly mixed nodes among the four real-world multi-layer networks analyzed.

• 3.

  The C.Elegans network has the lowest proportion of highly mixed nodes and the highest proportion of highly pure nodes. Its balanced parameter, $\upsilon$, is 0.9512, the highest among all networks, indicating that the sizes of the two estimated communities in C.Elegans are nearly identical.

• 4.

  The CS-Aarhus network exhibits indices that are comparable to those of the Lazega Law Firm network.

• 5.

  FAO-trade displays a proportion of highly mixed (and pure) nodes similar to C.Elegans. However, its balanced parameter is the lowest among all networks, indicating that FAO-trade exhibits the most unbalanced community structure among the four real-world datasets.

Figure 3 presents a ternary diagram that visualizes the estimated community membership matrix $\hat{\Pi}$ obtained from the SPSum algorithm for the Lazega Law Firm network when there are three communities. In this diagram, we observe that node $i$ is positioned closer to one of the triangle’s vertices compared to node $j$ if $\mathrm{max}_{k\in[3]}\hat{\Pi}(i,k)>\mathrm{max}_{k\in[3]}\hat{\Pi}(j,k)$ for $i\in[71],j\in[71]$. A pure node is located at one of the triangle’s vertices, and a neutral node is closer to a vertex than a highly mixed node. Therefore, Figure 3 indicates the purity of each node for the Lazega Law Firm network.

Figures 4-7 present the communities estimated by our SPSum method for the four real multi-layer networks. In these figures, nodes sharing the same color represent nodes belonging to the same home base community, while black squares denote highly mixed nodes. From these figures, we can clearly find the communities detected by SPSum in each layer for each multi-layer network.

This paper considers the problem of estimating community memberships of nodes in multi-layer networks under the multi-layer mixed membership stochastic block model, a model that permits nodes to belong to multiple communities simultaneously. We have developed spectral methods leveraging the eigen-decomposition of some aggregate matrices, and have provided theoretical guarantees for their consistency, demonstrating the convergence of per-node error rates as the number of nodes and/or layers increases under MLMMSB. To the best of our knowledge, this is the first work to estimate the mixed community memberships for multi-layer networks by using these aggregate matrices. Our theoretical analysis reveals that the algorithm designed based on the debiased sum of squared adjacency matrices always outperforms the algorithm using the sum of squared adjacency matrices, while they generally outperform the method using the sum of adjacency matrices in the task of estimating mixed memberships for multi-layer networks. Such a result is new for mixed membership estimation in multi-layer networks to the best of our knowledge. Extensive simulated studies support our theoretical findings, validating the efficiency of our method using the debiased sum of adjacency matrices. Additionally, the proposed fuzzy modularity measures offer a novel perspective for evaluating the quality of mixed membership community detection in multi-layer networks.

For future research, first, developing methods with theoretical guarantees for estimating the number of communities in MLMMSB remains a challenging and meaningful task. Second, accelerating our methods for detecting mixed memberships in large-scale multi-layer networks is crucial for practical applications. Third, exploring more efficient algorithms for estimating mixed memberships would further enrich our understanding of community structures in multi-layer networks. Finally, extending our framework to directed multi-layer networks would broaden the scope of our work and enable the analysis of even more complex systems.