Spectral Estimation of Large Stochastic Blockmodels with Discrete Nodal Covariates

Abstract

In many applications of network analysis, it is important to distinguish between observed and unobserved factors affecting network structure. We show that a network model with discrete unobserved link heterogeneity and binary (or discrete) covariates corresponds to a stochastic blockmodel (SBM). We develop a spectral estimator for the effect of covariates on link probabilities, exploiting the correspondence of SBMs and generalized random dot product graphs (GRDPG). We show that computing our estimator is much faster than standard variational expectation–maximization algorithms and scales well for large networks. Monte Carlo experiments suggest that the estimator performs well under different data generating processes. Our application to Facebook data shows evidence of homophily in gender, role and campus-residence, while allowing us to discover unobserved communities. Finally, we establish asymptotic normality of our estimators.

KEYWORDS:

• Asymptotic normality
• Conditionally independent links
• Generalized random dot product graphs
• Large networks
• Spectral estimation
• Stochastic blockmodel

Acknowledgments

We are grateful to Cong Mu and Jipeng Zhang for excellent research assistance. We thank the editor, associate editor and two referees, Avanti Athreya, Eric Auerbach, Federico Bandi, Stephane Bonhomme, Youngser Park and Eleonora Patacchini for comments and suggestions.

Funding

Screeplots (upper left and center left), Estimated latent positions $\widehat{\mathit{Y}}$ (right, only 2 dimensions out of 4 per plot) and estimated latent positions $\widehat{\mathit{X}}$, that is $\widehat{p}$ and $\widehat{q}$ (bottom left, up to orthogonal transformation) for n = 2000.

Notes

1 Recent advances use further approximations and parallelization to improve computational efficiency (Roy, Atchade, and Michailidis 2019; Vu, Hunter, and Schweinberger 2013). We do not pursue such extensions in this article.

2 It must be noted that the support ${\mathit{X}}_d$ of F, is a subset of ${\mathbb{R}}^d$ such that ${\mathit{x}}^{\top }{\mathit{I}}_{d_1,d_2}\mathit{y}\in [0,1]$ for all $\mathit{x},\mathit{y}\in {\mathit{X}}_d$.

3 Alternatively, we can think of a network where the vectors ${\mathit{X}}_i$ are drawn from a discrete mixture with mass centered at $\mathit{\nu }$, that is,

${\mathit{X}}_i\sim {\pi }_1{\delta }_{{\mathit{\nu }}_1}+{\pi }_2{\delta }_{{\mathit{\nu }}_2}+\dots +{\pi }_K{\delta }_{{\mathit{\nu }}_K}.$ (4)

4 The extensions to multiple binary or discrete covariates is shown later in the article.

5 In the Bernoulli case, ${\mathit{E}}_{ij}$ is a shifted Bernoulli variable, with values ${\mathit{E}}_{ij}=1-{\mathit{P}}_{ij}$ with probability ${\mathit{P}}_{ij}$ and ${\mathit{E}}_{ij}={\mathit{P}}_{ij}$ with probability $1-{\mathit{P}}_{ij}$.

6 A possible alternative is to exploit the overidentification of β in the matrix ${\widehat{\mathit{\theta }}}_Z$ to develop a minimum distance estimator. We do not pursue this direction here, as we focus on a fully spectral estimator.

7 The notation $n{\rho }_n=\omega \left(\sqrt{n}\right)$ means that for any real constant a > 0 there exists an $n_0\ge 1$ such that ${\rho }_n>a/\sqrt{n}\ge 0$ for every integer $n\ge n_0$.

8 We do not compare the spectral estimator to the variational EM estimator, because the latter is too slow for a Monte Carlo with 1000 repetitions, even after parallelizing the execution.

9 The time of estimation reported in the table includes the following steps: (a) compute the ASE from the adjacency matrix; (b) compute the matrix of latent positions; (c) cluster latent positions to recover blocks; (d) compute matrix ${\widehat{\mathit{B}}}_Z$; (e) cluster diagonal entries of matrix ${\widehat{\mathit{B}}}_Z$ to recover the unobservable block structure; (f) estimate $\mathit{\beta }$ using the information on the block structure and the entries of matrix ${\widehat{\mathit{B}}}_Z$; (g) compute simple mean and weighted mean of the estimated $\widehat{\mathit{\beta }}$ according to (41) and (43). The simulation takes a little longer because we need to generate the data and the adjacency matrices for the Monte Carlo. Code is available in GitHub.

10 The entire dataset is available at https://archive.org/details/oxford-2005-facebook-matrix.

11 Roles include students, faculty, staff, alumni. We focus on students because they are the ones mostly using the platform in 2005.

12 This is a standard procedure in the literature on SBMs (Athreya et al. 2018; Abbe 2018; Roy, Atchade, and Michailidis 2019).

13 Before we proceed to estimation, we regularize the adjacency matrix using the standard method proposed in (Le, Levina, and Vershynin 2017). This regularization step avoids numerical issues with the spectral decomposition arising from significant node degree heterogeneity.

14 Multiple methods exist for selecting the embedding dimension in practice, and this remains a topic of current research. In the context of networks, choosing a dimension smaller than the true d will introduce bias in the estimated latent positions; on the other hand, using an embedding dimension larger than the true d will increase the variance of the estimated latent positions. In this tradeoff we prefer to err on the side of overestimating d. Specifically, we choose the value one plus the location of the first elbow in the screeplot.

Additional information

Funding

Funding from the Institute of Data Intensive Engineering and Science (IDIES) at Johns Hopkins University and NSF grant SES-1951005 is gratefully acknowledged. Joshua Cape also gratefully acknowledges support from NSF grant DMS-1902755.