Limit Results for Estimation of Connectivity Matrix in Multi-layer Stochastic Block Models

The multi-layer network has gained widespread attention for its ability to represent multiple relationships among the same entities of interest, thereby enhancing the understanding of complex network data (Mucha et al., 2010; Holme and Saramäki, 2012; Kivelä et al., 2014; Han et al., 2015). For example, the Worldwide Food and Agricultural Trade (WFAT) multi-layer trade network (De Domenico et al., 2015), encompassing trade relationships among the same countries (nodes) across different commodities (layers), offers a richer insight than a single-layer trade network.

The multi-layer stochastic block model (multi-layer SBM) (Han et al., 2015; Paul and Chen, 2016) has been popularly used for the community detection task in the multi-layer network, where each layer corresponds to a stochastic block model (SBM) (Holland et al., 1983). In particular, given common unobserved communities (blocks), the edges in each layer are generated according to a layer-wise block probability matrix or the so-called connectivity matrix. Various methods have been proposed to detect the communities of multi-layer networks under the multi-layer SBM, see, e.g., Han et al. (2015); Bhattacharyya and Chatterjee (2018); Paul and Chen (2020); Jing et al. (2021); Arroyo et al. (2021); Lei and Lin (2023).

Most of the aforementioned literature on community detection of multi-layer SBMs is devoted to the statistical consistency. However, the asymptotic distributional properties are also vital and particularly useful in subsequent inference tasks. In the context of single-layer SBMs or other related models, various versions of asymptotically Gaussian behavior emerged for different purposes (Agterberg and Cape, 2023). Specifically, Tang and Priebe (2018) provided a central limit theorem for the eigenvectors of the normalized Laplacian matrix under the random dot product graph (RDPG) model (Young and Scheinerman, 2007). Rubin-Delanchy et al. (2022) extended these asymptotic results to the generalized RDPG and applied the theory to SBMs. Bickel et al. (2013) and Tang et al. (2022) established the asymptotic normality results for the estimation of the connectivity matrix using the maximum likelihood method and spectral embedding, respectively. Fan et al. (2022b) tested whether two nodes from a mixed-membership model (Airoldi et al., 2008) have the same membership parameters by establishing the asymptotic distributional properties. Jin et al. (2023) proposed a goodness-of-fit approach to select the number of communities and the asymptotic null distribution of the test statistic is normal. Also see Agterberg and Cape (2023) for a comprehensive overview of asymptotic normality in SBMs.

In the multi-layer SBM, the asymptotic distributional properties are crucial for subsequent inferences especially for those tailored for the multi-layer networks. However, as far as we are aware, the asymptotic properties under multi-layer SBMs are largely unexplored. An important inference task in multi-layer SBMs is to test whether networks across different layers are generated from the same SBM, essentially the same connectivity matrix given the consensus community memberships. For instance, in analyzing the WFAT multi-layer network, the goal extends beyond the community detection to ascertaining whether trading patterns across different commodities are similar. Asymptotic properties play a fundamental role in such inference. In the literature, Arroyo et al. (2021) and Zheng and Tang (2022) showed that the connectivity matrices estimated using the spectral embedding are asymptotically normal under the common subspace independent edge graph (COSIE) model. These studies necessitate that the connectivity matrices for all network layers are full rank. However, in practical multi-layer networks like WFAT, not all layers carry comprehensive cluster information, potentially leading to missing clusters, which results in the rank-deficiency of individual connectivity matrices (Su et al., 2024). Moreover, the methods in Arroyo et al. (2021) and Zheng and Tang (2022) treat multi-layer SBMs as special case of their COSIE model, rather than being specifically tailored for them.

Motivated by the above problems, we study the asymptotic properties in the estimation of connectivity matrices under multi-layer SBMs, without imposing the full rank assumption on their populations. We develop a simple and efficient spectral clustering-based method for the scaled connectivity matrix estimation, where the scaled connectivity matrix is also called the score matrix (Arroyo et al., 2021). Under proper conditions, we establish the asymptotic normality of the estimated scaled connectivity matrices under multi-layer SBMs and its variants, the degree-corrected multi-layer SBM. We emphasize that the asymptotic properties can be used to perform various statistical inferences, such as interval estimation and hypothesis testing.

The main contributions of this paper are threefold. First, we systematically study and specifically tailor the asymptotic normality of the scaled connectivity matrix estimation under the multi-layer SBM and its variant, the multi-layer degree-corrected SBM. To the best of our knowledge, the asymptotic normality under the degree-corrected SBM is less explored, let alone the multi-layer degree-corrected SBM. Second, the conditions under which our results hold are relatively mild. Compared to previous work which necessitated full rank in each layer (Arroyo et al., 2021; Zheng and Tang, 2022), we allow individual connectivity matrices to be rank-deficient and only require their squared summation to be full rank. This adaptability makes our approach particularly suitable for analyzing multi-layer networks where each layer only captures partial underlying communities. Third, we apply the established asymptotic normality to two statistical inference tasks, one for interval estimation of the scaled connectivity matrices, and the other for testing whether the population matrices corresponding to different layers are the same. The numerical results show superior performance of the proposed method over compared methods in these tasks. The analysis of a real-world trade network dataset provides interpretable results.

The rest of the paper is organized as follows. Section 2 introduces the estimation method of individual connectivity matrices in the multi-layer SBM. Section 3 provides the asymptotic normality of the estimates. Section 4 extends the method to multi-layer degree-corrected SBMs and presents the corresponding asymptotic normality results. Sections 5 and 6 include the simulation and real data experiment. Section 7 concludes the paper and provides possible extensions. Technical proofs are included in the Appendix.

In this section, we first present and reparameterize the multi-layer stochastic block model. Then we propose the method for estimating the individual connectivity matrices. The following notation is used throughout the paper.

In this section, we present the asymptotic properties of the individual estimates of the scaled connectivity matrices $M_{l}(l\in[L])$ in the multi-layer SBM.

To establish the asymptotic normality of $M_{l}$, we need the following assumptions.

Assumption 1 also indicates that there exist positive constants $c_{2},c_{3}$, and an orthogonal matrix $Z\in\mathbb{O}(K)$, such that $c_{2}/\sqrt{n}\leq|(UZ)_{is}|\leq c_{3}/\sqrt{n}$ for all $i\in[n]$ and $s\in[K]$. This is referred to as eigenvector delocalization (Rudelson and Vershynin, 2015; He et al., 2019). Indeed, by recalling the relationship $U=\Theta\Delta^{-1/2}$ and choosing the orthogonal matrix $Z$ with $c_{2}/\sqrt{K}\leq Z_{st}\leq c_{3}/\sqrt{K}$ for each $K$ and each pair $1\leq s,t\leq K$, the delocalization of $U$ follows by Assumption 1.

Assumption 2 specifies the growth rate of the minimum eigenvalue of $\sum_{l=1}^{L}B_{l}^{2}$ in order to make the theoretical bound of the estimated common eigenspace $\hat{U}$ around $U$ concrete. The linear growth rate is reasonable because when the connectivity matrices are common, Assumption 2 holds naturally.

Assumption 3 is critical to ensure the Lindeberg conditions for the Central Limit Theorem are met. For balanced community sizes, Assumption 3 simplifies to $\sum_{s=1}^{K}\sum_{t=1}^{K}\rho B_{l,st}(1-\rho B_{l,st})=\omega(1/n^{2})$, which suffices if $\rho=\omega(1/n^{2})$.

Before stating the next assumption, we first define the following $K(K+1)/2\times K(K+1)/2$ matrix $\Sigma^{l}$,

$$\Sigma_{\frac{2s+t(t-1)}{2},\frac{2s^{\prime}+t^{\prime}(t^{\prime}-1)}{2}}^{l}:=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}(U_{is}U_{jt}+U_{js}U_{it})\;(U_{is^{\prime}}U_{jt^{\prime}}+U_{js^{\prime}}U_{it^{\prime}})\;Q_{l,ij}(1-Q_{l,ij})$$

(3)

for each $1\leq s\leq t\leq K$ and $1\leq s^{\prime}\leq t^{\prime}\leq K$. Also, define $G_{l}=A_{l}-P_{l}$ with $P_{l}=\mathbb{E}(A_{l})=Q_{l}-\text{diag}(Q_{l})$ to be the noise matrix. As we will see, $\Sigma^{l}$ actually serves as the covariance matrix of the vectorized $U^{T}G_{l}U$, denoted by ${\rm vec}(U^{T}G_{l}U)\in\mathbb{R}^{K(K+1)/2}$, where we focus on the upper triangular part of the symmetric $U^{T}G_{l}U$. Specifically, for any $1\leq s\leq t\leq K$, the $\frac{2s+t(t-1)}{2}$th entry in ${\rm vec}(U^{T}G_{l}U)$ corresponds to $(U^{T}G_{l}U)_{s,t}$.

Similar to Assumption 3, Assumption 4 contributes to verifying the Lindeberg conditions for the Central Limit Theorem. By Assumption 3 and the discussions after Assumption 1, we can conclude that each entry of $\Sigma_{l}$ is $\omega(1/n^{2})$. Hence, Assumption 4 is stronger than Assumption 3.

With these assumptions, we obtain the following entry-wise and vector-wise asymptotic properties of the estimated scaled connectivity matrices under the multi-layer SBM.

In this section, we extend the proposed method and the corresponding asymptotic normality results to multi-layer degree-corrected stochastic block models, which is a counterpart of the multi-layer SBM but in each layer the network is assumed to be generated from a degree-corrected SBM (DCSBM) (Karrer and Newman, 2011).

The SBM can not capture the degree heterogeneity inherent in the networks. To address this, the DCSBM extends the standard SBM by incorporating node specific parameters, allowing degrees to vary within the same community.

We now introduce the multi-layer DCSBM, where the layer-wise networks are generated from the DCSBM. Without specification, the notes and notation are the same as those in Section 2. Define $\psi\in\mathbb{R}^{n}_{+}$ to be the degree heterogeneity parameter which measures the propensity of a node in forming edges with other nodes and is consensus among layers. For identifiability, we assume $\max_{i\in C_{k}}\psi_{i}=1$ for all $k\in[K]$ (Lei and Rinaldo, 2015; Zhang et al., 2022). Given $\psi$, the community assignments $g_{i}\in[K]$, the connectivity matrix $B_{l}(l\in[L])$, and the network sparsity parameter $\rho$, the element $A_{l,ij}(i<j)$ is generated independently as follows

$$A_{l,ij}\sim{\rm Bernoulli}(\rho\psi_{i}\psi_{j}B_{l,g_{i}g_{j}}),$$

and $A_{l,ij}=A_{l,ji}$ and $A_{l,ii}=0$. The population adjacency matrix of $A_{l}$ is then

$$Q_{l}=\rho{\rm diag}(\psi)\Theta B_{l}\Theta^{T}{\rm diag}(\psi).$$

To facilitate further analysis, we now give some additional notation. For each $k\in[K]$, define $G_{k}=\{1\leq i\leq n:g_{i}=k\}$ as the set of nodes whose community membership is $k$. Let $\phi_{k}$ be an $n\times 1$ vector that matches $\psi$ on the index set $G_{k}$ and is zero elsewhere. Define $\Psi={\rm diag}(\|\phi_{1}\|_{2}^{2},\ldots,\|\phi_{K}\|_{2}^{2})$, where each diagonal element $\|\phi_{k}\|_{2}^{2}$ represents the effective community size of community $k$. This allows us to reformat $Q_{l}$ as

$$Q_{l}=\rho{\rm diag}(\psi)\Theta\Psi^{-1/2}\Psi^{1/2}B_{l}\Theta^{T}\Psi^{1/2}\Psi^{-1/2}{\rm diag}(\psi):=UM_{l}U^{T},$$

where $U={\rm diag}(\psi)\Theta\Psi^{-1/2}$ and $M_{l}=\rho\Psi^{1/2}B_{l}\Psi^{1/2}$. It can be shown that $U^{T}U=I$.

As described in the case of multi-layer SBMs, we regard $U$ as the eigenspace of the following $\sum_{l=1}^{L}Q_{l}^{2}$,

$$\sum_{l=1}^{L}Q_{l}^{2}=\sum_{l=1}^{L}\rho^{2}{\rm diag}(\psi)\Theta B_{l}\Theta^{T}{\rm diag}^{2}(\psi)\Theta B_{l}\Theta^{T}{\rm diag}(\psi)=U\sum_{l=1}^{L}\rho^{2}\Psi^{1/2}B_{l}\Psi B_{l}\Psi^{1/2}U^{T}.$$

Suppose $\sum_{l=1}^{L}B_{l}^{2}$ is of full rank, then $\sum_{l=1}^{L}\rho^{2}\Psi^{1/2}B_{l}\Psi B_{l}\Psi^{1/2}$ is of full rank as well. Similar to the case of multi-layer SBMs, we perform the eigendecomposition on the bias-adjusted sum of squares of the adjacency matrices $\sum_{l=1}^{L}(A_{l}^{2}-D_{l})$, where recall $D_{l}$ is an $n\times n$ diagonal matrix with $D_{l,ii}=\sum_{j=1}^{n}A_{l,ij}$ and the resulting eigenvectors are denoted by $\hat{U}$. We further obtain $\hat{M_{l}}$ with the help of $\hat{U}$ through (2).

The following Assumptions E1, E3, and E4 are needed for establishing the asymptotic normality of $\hat{M}_{l}$, which are counterparts of the assumptions under the multi-layer SBM.

Assumption E1 requires that the node propensity parameters restricted to each community have the same order of Euclidean norm. This condition is frequently imposed in the analysis of the DCSBM; see Su et al. (2020); Jin et al. (2023) and Agterberg et al. (2022). Generally, Assumptions E3 and E4 are more stringent than Assumptions 3 and 4 because $\|\psi\|_{2}\lesssim\sqrt{n}$; while in the special case of multi-layer SBMs where $\psi=\mathbf{1}$ and $\|\psi\|_{2}\asymp\sqrt{n}$, Assumptions E3 and E4 reduce to Assumptions 3 and 4, respectively. With these assumptions, the entry-wise and vector-wise asymptotic normality hold for the estimated scaled connectivity matrix $\hat{M}_{l}$ under the multi-layer DCSBM.

In this section, we conduct simulations to investigate the finite sample performance of the proposed method. In Section 5.1, we verify the vanishing behavior of the bias term $E_{l}$ in the formulation of asymptotic normality. In Sections 5.2 and 5.3, we test the efficacy of the derived asymptotic normality in two downstream statistical inference tasks, one for the interval estimation of the entries of the scaled connectivity matrices and the other for hypothesis testing to infer whether the multi-layer network exhibits homogeneity.

The proposed method is denoted by SCCE, namely, the abbreviation for spectral clustering-based method for connectivity matrix estimation. Since studies on the asymptotic properties under multi-layer SBMs are relatively limited, we compare SCCE with the Multiple Adjacency Spectral Embedding (MASE) method studied in Arroyo et al. (2021) and Zheng and Tang (2022).

In this section, we use the proposed method SCCE to measure the homogeneity across all layers of a real multi-layer network, specifically the WFAT dataset (De Domenico et al., 2015). For this purpose, we simultaneously test the $\binom{L}{2}$ hypotheses

$$H_{kl,0}:M_{k}=M_{l},\quad\quad 1\leq k<l\leq L,$$

using a Holm type step-down procedure. Next, we provide the data description and the results for the real dataset.

The WFAT dataset, sourced from the Food and Agriculture Organization (FAO) of the United Nations, contains annual trade records for over 400 food and agricultural products across all countries globally. Our analysis focuses specifically on 2020 trade data for cereals, a category defined by the FAO and detailed in Table A3. We constructed a multi-layer network using this data. Each layer represents a different type of cereal, with nodes representing countries and edges within each layer representing the trade relationships with respect a specific cereal. In particular, a trade link between countries is established if the trade value for a cereal exceeds $\$2000$. To ensure network connectivity, we excluded countries with a total degree over 37 layers less than 23. This guarantees that each node is linked to at least one other node in at least half of the layers. The resulting network consists of 37 layers, each containing 114 common nodes, denoted by $A_{l},l\in\{1,\ldots,37\}$. The five continents, including America, Africa, Asia, Europe, and Oceania, suggest the choice $K=5$ (Noroozi and Pensky, 2024).

Figure 4 shows the results of the multiple comparisons conducted using a Holm type step-down procedure, where olive green indicates the acceptance and white denotes the rejection of the individual hypotheses. We control the Holm procedure to ensure that the family-wise error rate is no greater than $\alpha=0.05$. Additionally, the $p$-value matrix for individual hypothesis tests across all network layers is shown. In Figure 4(a), several distinct blocks are observed where all individual hypotheses within each block are simultaneously accepted, demonstrating the presence of homogeneity within the multi-layer network. For instance, layers 5 to 14, corresponding to ‘Germ of maize’, ‘Triticale’, ‘Cereals n.e.c.’, ‘Millet’, ‘Sorghum’, ‘Bran of maize’, ‘Bran of cereals n.e.c.’, ‘Flour of mixed grain’, ‘Flour of rice’, and ‘Rye’, exhibit no rejections of the individual hypotheses in any pairwise comparisons. This outcome suggests a significant similarity in trade patterns among these cereal products. Notably, most of these cereals are consumed primarily by Asia, which might lead to simultaneous acceptance of the individual hypotheses. This observation suggests the need for an integrated analysis of the trade patterns of these cereals.

In this paper, we focused on the asymptotic normality in the estimation of the scaled connectivity matrix under multi-layer SBMs and multi-layer DCSBMs. To this end, we first leveraged the leading eigenvectors of the bias-adjusted sum of squares of the adjacency matrices to estimate the common eigenspaces across layers (Lei and Lin, 2023). Building on this, we further derived efficient estimates for the layer-wise scaled connectivity matrices. These estimates are shown to satisfy asymptotic normality under mild additional assumptions. To corroborate the theoretical results, we conducted a series of simulations to apply the derived asymptotic normality to statistical inference tasks, including interval estimation and hypothesis testing. Finally, we applied the method to a real-world network dataset to infer the population homogeneity within the multi-layer network.

There are many ways to extend the content of this paper. First, we assume the connectivity matrices are different among layers. However, in practice, the network layers are often grouped such that connectivity matrices are the same within each group (Jing et al., 2021; Fan et al., 2022a). It is thus meaningful to make use of all the network layers within each group to estimate the connectivity matrices, which would lead to improved estimation accuracy and a weakened requirement for network sparsity. Second, we focused on the undirected multi-layer networks. It is of great importance to extend the method and theory to multi-layer directed networks (Zheng and Tang, 2022; Su et al., 2024). Third, we assumed the number of communities is known and fixed. In the future, it is important to study the estimation of the number of communities under multi-layer SBMs.