Bayesian estimation of the latent dimension and communities in stochastic blockmodels

TL;DR

This paper introduces a Bayesian model that automatically determines the latent space dimension and number of communities in stochastic blockmodels, improving community detection in networks.

Contribution

A novel Bayesian approach for simultaneous, automatic estimation of latent dimension and community count in stochastic blockmodels, extending to directed and bipartite graphs.

Findings

Effective in recovering community structure from simulated data

Shows promising results on real-world network data

Outperforms traditional methods requiring pre-specified parameters

Abstract

Spectral embedding of adjacency or Laplacian matrices of undirected graphs is a common technique for representing a network in a lower dimensional latent space, with optimal theoretical guarantees. The embedding can be used to estimate the community structure of the network, with strong consistency results in the stochastic blockmodel framework. One of the main practical limitations of standard algorithms for community detection from spectral embeddings is that the number of communities and the latent dimension of the embedding must be specified in advance. In this article, a novel Bayesian model for simultaneous and automatic selection of the appropriate dimension of the latent space and the number of blocks is proposed. Extensions to directed and bipartite graphs are discussed. The model is tested on simulated and real world network data, showing promising performance for recovering latent community structure.

Equations

$$L(\mathbf{A}|\mathbf{z},\mathbf{B})=\prod _{1\le i<j\le n}{B}_{z_i{z}_j}^{A_{ij}}(1-B_{z_i{z}_j}{)}^{1-A_{ij}}.$$

$$\mathbf{A}=\hat{\mathbf{\Gamma }}\hat{\mathbf{\Lambda }}{\hat{\mathbf{\Gamma }}}^{\top }+{\hat{\mathbf{\Gamma }}}_{\perp }{\hat{\mathbf{\Lambda }}}_{\perp }{\hat{\mathbf{\Gamma }}}_{\perp }^{\top },$$

$$\hat{\mathbf{X}}=\hat{\mathbf{\Gamma }}|\hat{\mathbf{\Lambda }}{|}^{1/2}\in {\mathbb{R}}^{n\times d},$$

$$\mathbf{L}={\mathbf{D}}^{-1/2}\mathbf{A}{\mathbf{D}}^{-1/2},\mathbf{D}=\mathrm{diag}\left({\sum }_{j=1}^n{A}_{ij}\right),$$

$$\mathbf{L}=\tilde{\mathbf{\Gamma }}\tilde{\mathbf{\Lambda }}{\tilde{\mathbf{\Gamma }}}^{\top }+{\tilde{\mathbf{\Gamma }}}_{\perp }{\tilde{\mathbf{\Lambda }}}_{\perp }{\tilde{\mathbf{\Gamma }}}_{\perp }^{\top }.$$

$$\tilde{\mathbf{X}}=\tilde{\mathbf{\Gamma }}|\tilde{\mathbf{\Lambda }}{|}^{1/2}.$$

$$\mathbf{I}(p,q)=\mathrm{diag}(\underset{p}{\underbrace{1,\dots ,1}},\underset{q}{\underbrace{-1,\dots ,-1}}).$$

$$\mathbb{P}(A_{ij}=1)={\mathbf{x}}_i^{\top }\mathbf{I}(d_{+},d_{-}){\mathbf{x}}_j.$$

$${\mathbf{Q}}_m{\hat{\mathbf{x}}}_i\stackrel{d}{⟶}\mathbb{N}\{{\mathbf{\mu }}_{z_i},m^{-1/2}\mathbf{\Sigma }({\mathbf{\mu }}_{z_i})\},m\to \infty ,i=1,\dots ,n,$$

$${\mathbf{x}}_i|d,z_i,{\mathbf{\mu }}_{z_i},{\mathbf{\Sigma }}_{z_i},{\mathbf{\sigma }}_{z_i}^2$$

$$({\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)|d$$

$${\sigma }_{k,j}^2$$

$$d|\mathbf{z}$$

$$z_i|\mathbf{\theta }$$

$$\mathbf{\theta }|K$$

$$K$$

$$d\stackrel{d}{\sim }\mathrm{Geometric}(\delta ),$$

$$L(\mathbf{X})=\prod _{i=1}^n\left\{\varphi ({\mathbf{x}}_{i,:d};{\mathbf{\mu }}_{z_i},{\mathbf{\Sigma }}_{z_i})\prod _{j=d+1}^m\varphi (x_{i,j};0,{\sigma }_{z_i,j}^2)\right\},$$

$$p(\{\bm{\mu}_{k}\},\{\bm{\Sigma}_{k}\},\{\bm{\sigma}^{2}_{k}\},\bm{z},\bm{\theta},K,d|\mathbf{X})\propto L(\mathbf{X})\prod_{k=1}^{K}\bigg{\{}p(\bm{\mu}_{k},\bm{\Sigma}_{k}|d)\prod_{j=d+1}^{m}p(\sigma^{2}_{k,j}|d)\bigg{\}}\prod_{i=1}^{n}p(z_{i}|\bm{\theta})p(K)p(d).$$

$$p({\mathbf{\mu }}_k|\mathbf{X},\mathbf{z},d)=t_{{\nu }_{n_k}}\left({\mathbf{\mu }}_k|{\mathbf{m}}_{:d}^{(k)},{\mathbf{D}}_{:d}^{(k)}/({\kappa }_{n_k}{\nu }_{n_k})\right),$$

$${\nu }_{n_k}$$

$${\kappa }_{n_k}$$

$${\mathbf{m}}_{:d}^{(k)}$$

$${\mathbf{D}}_{:d}^{(k)}$$

$$p({\mathbf{X}}_{:d}^{(k)}|d,\mathbf{z})={\pi }^{-n_{k}d/2}\frac{{\kappa }_0^{d/2}|{\mathbf{\Delta }}_d{|}^{({\nu }_0+d-1)/2}}{{\kappa }_{n_k}^{d/2}|{\mathbf{D}}_{:d}^{(k)}{|}^{({\nu }_{n_k}+d-1)/2}}\prod _{i=1}^d\frac{\Gamma \{({\nu }_{n_k}+d-i)/2\}}{\Gamma \{({\nu }_0+d-i)/2\}},$$

$$\lambda_{n_{k}}=\lambda_{0}+n_{k},\quad s^{(k)}_{j}=\left\{\lambda_{0}\sigma^{2}_{0}+\sum\nolimits_{i:z_{i}=k}x_{ij}^{2}\right\}\Big{/}\lambda_{n_{k}}.$$

$$p(\mathbf{X}^{(k)}_{d:}|d,\bm{z})=\pi^{-n_{k}(m-d)/2}\bigg{\{}\frac{\Gamma(\lambda_{n_{k}}/2)}{\Gamma(\lambda_{0}/2)}\bigg{\}}^{m-d}\prod_{j=d+1}^{m}\frac{(\lambda_{0}\sigma^{2}_{0})^{\lambda_{0}/2}}{(\lambda_{n_{k}}s^{(k)}_{j})^{\lambda_{n_{k}}/2}}.$$

$$p(\mathbf{X}|\mathbf{z},K)=\sum _{d=1}^{K}p(d|K)\prod _{k=1}^{K}p({\mathbf{X}}_{:d}^{(k)}|d,\mathbf{z})p({\mathbf{X}}_{d:}^{(k)}|d,\mathbf{z}),$$

$$p(\mathbf{z}|K)=\frac{\Gamma (\alpha ){\prod }_{k=1}^K\Gamma (n_k+\alpha /K)}{\Gamma (\alpha /K{)}^K\Gamma (n+\alpha )}.$$

$$p(z_i=k|{\mathbf{z}}_{-i},\mathbf{X},d,K)\propto p(z_i=k|{\mathbf{z}}_{-i},d,K)p({\mathbf{x}}_i|\{{\mathbf{x}}_j{\}}_{j\ne i:z_j=k},d).$$

$$p(z_i=k|{\mathbf{z}}_{-i},d,K)\propto \frac{n_k^{-i}+\alpha /K}{n-1+\alpha }.$$

$$p({\mathbf{x}}_i|\{{\mathbf{x}}_j{\}}_{j\ne i:z_j=k},d)=\frac{p({\mathbf{x}}_i,\{{\mathbf{x}}_j{\}}_{j\ne i:z_j=k}|d)}{p(\{{\mathbf{x}}_j{\}}_{j\ne i:z_j=k}|d)}.$$

$$p({\mathbf{x}}_{i,:d}|\{{\mathbf{x}}_{j,:d}{\}}_{j\ne i:z_j=k},d)=t_{{\nu }_{n_k^{-i}}}({\mathbf{x}}_{i,:d}\left|{\mathbf{m}}_{:d}^{(k)},\frac{{\kappa }_{n_k^{-i}}+1}{{\kappa }_{n_k^{-i}}{\nu }_{n_k^{-i}}}{\mathbf{D}}_{:d}^{(k),-i}\right).$$

$$p({\mathbf{x}}_{i,d:}|\{{\mathbf{x}}_{j,d:}{\}}_{j\ne i:z_j=k},d)=\prod _{j=d+1}^m{t}_{{\lambda }_{n_k^{-i}}}(x_{ij}\left|0,s_j^{(k),-i}\right).$$

$$\alpha (K^{\star },{\mathbf{z}}^{\star }|K,\mathbf{z})=\min \left\{1,\frac{p(\mathbf{X}|d,K^{\star },{\mathbf{z}}^{\star })p(d|{\mathbf{z}}^{\star },K^{\star })p({\mathbf{z}}^{\star }|K^{\star })p(K^{\star })}{p(\mathbf{X}|d,K,\mathbf{z})p(d|\mathbf{z},K)p(\mathbf{z}|K)p(K)q(K^{\star },{\mathbf{z}}^{\star }|K,\mathbf{z})}\right\}.$$

Full text

1 Introduction

A network can be represented as a graph $\mathbb{G}=(V,E)$, where $V$ is a set of nodes and $E\subseteq V\times V$ is a set of edges indicating the pairs of nodes which have interacted. The graph can be characterised by the adjacency matrix $\mathbf{A}\in\{0,1\}^{n\times n}$, where $n=|{V}|$ and for $1\leq i,j\leq n$, $A_{ij}=\mathds{1}_{E}\{(i,j)\}$, such that $A_{ij}=1$ if a link between the nodes $i$ and $j$ exists, and $A_{ij}=0$ otherwise. The graph is said to be undirected if $(i,j)\in E\iff(j,i)\in E$ and $\mathbf{A}$ is constrained to be symmetric; otherwise, the graph is said to be directed. It will be assumed that a node cannot link to itself, implying $\mathbf{A}$ is a hollow matrix.

Latent space models (Hoff et al., 2002; Handcock et al., 2007; Krivitsky et al., 2009) represent a flexible approach to statistical analysis of networks: each node $i$ is assigned a latent position $\bm{x}_{i}$ in a $d$-dimensional latent space $\mathbb{X}$, and edges between pairs of nodes are typically generated independently (independent-edge graph), with the probability of observing a link between nodes $i$ and $j$ obtained through a kernel function $\psi:\mathbb{X}\times\mathbb{X}\to[0,1]$ of the respective latent positions: $\mathbb{P}(A_{ij}=1)=\psi(\bm{x}_{i},\bm{x}_{j})$. Different ideas and techniques for embedding observed graphs into low dimensional spaces are explored in the literature (for a survey, see, for example Cai et al., 2018). Random dot product graphs (RDPGs) (Nickel, 2006; Young and Scheinerman, 2007, 2008) are a popular class of latent position models, where $\mathbb{X}\subseteq\mathbb{R}^{d}$, and the function $\psi(\cdot)$ is an inner product $\langle\cdot,\cdot\rangle$ on $\mathbb{X}\times\mathbb{X}$. RDPGs are analytically tractable and have therefore been extensively studied; a survey of the existing statistical inference techniques is presented in Athreya et al. (2017).

The stochastic blockmodel (SBM) (Holland et al., 1983) is the classical statistical model for clustering graphs (Snijders and Nowicki, 1997; Nowicki and Snijders, 2001). Assuming $K$ communities, each node is assigned a community membership $z_{i}\in\{1,\dots,K\}$ with probabilities $\bm{\theta}$, from the $K-1$ probability simplex. The probability of a link only depends on the community allocations $z_{i}$ and $z_{j}$ of the two nodes. Given a symmetric matrix $\mathbf{B}\in[0,1]^{K\times K}$ of inter-community probabilities, then independently $\mathbb{P}(A_{ij}=1)=B_{z_{i}z_{j}}$. The likelihood for an observed symmetric adjacency matrix $\mathbf{A}$ is therefore

[TABLE]

Stochastic blockmodels have appealing statistical properties, and can well approximate any independent-edge network model if the number of communities is sufficiently large (Bickel and Chen, 2009; Wolfe and Olhede, 2013). Limits and bounds for community detection in stochastic blockmodels, from an information-theoretic approach, have been extensively studied and explored in the literature (Abbe et al., 2016; Abbe, 2018). SBMs have also been successfully extended in different directions: allowing for multiple overlapping communities (Airoldi et al., 2008), incorporating degree corrections (Karrer and Newman, 2011), capturing node popularity (Sengupta and Chen, 2018), Bayesian models (Peixoto, 2018; van der Pas and van der Vaart, 2018) and dynamically evolving models (Xu and Hero III, 2013; Matias and Miele, 2017; Pensky and Zhang, 2017; Ludkin et al., 2018). Stochastic blockmodels can also be easily represented as random dot product graphs: each community is assigned a latent position, which is common to all the nodes belonging to the cluster, and $\mathbf{B}$ is obtained from the inner products of those positions. Hence, in this framework, $d=\operatorname{rank}(\mathbf{B})\leq K$.

Spectral clustering (von Luxburg, 2007) provides a consistent statistical estimation procedure for the latent positions in SBMs (Rohe et al., 2011; Sussman et al., 2012; Fishkind et al., 2013; Lyzinski et al., 2014; Lei and Rinaldo, 2015; Lyzinski et al., 2017) and more generally in random dot product graphs (Tang et al., 2013; Sussman et al., 2014). Community memberships in the stochastic block model can be consistently estimated using $K$-means clustering of the rows (Rohe et al., 2011). Rubin-Delanchy et al. (2017) directly links spectral embedding to the generalised random dot product graph (GRDPG), an extension of the RDPG, and advocates for Gaussian mixture modelling (GMM) of the rows of the embedding, especially when the Laplacian matrix is used. Alternatives to spectral clustering include variational methods (Celisse et al., 2012) and pseudo-likelihood approaches (Amini et al., 2013). SBMs have been extended to the directed case (Wang and Wong, 1987; Rohe et al., 2016), and appropriate embeddings for co-clustering, in most cases based on the singular value decomposition (SVD), have been derived in the literature (Dhillon, 2001; Malliaros and Vazirgiannis, 2013; Zheng and Skillicorn, 2015).

One of the practical issues of spectral embedding, and in general all graph embedding algorithms (Cai et al., 2018), is that it requires a suitable prespecified latent dimensionality $d$ (usually $d\ll|{V}|$) as input, and, subsequently, a suitable number of clusters $K$, conditionally, crucially, on the previous choice of $d$. For a practical example of this procedure on a real world network, see Priebe et al. (2019). In general, in spectral clustering, similarly to what practitioners do in principal component analysis (PCA), the investigator examines the scree-plot of the eigenvalues and chooses the dimension based on the location of elbows in the plot (Jolliffe, 2002), or uses the eigengap heuristic (see, for example, von Luxburg, 2007). Automatic methods for thresholding have also been suggested (Zhu and Ghodsi, 2006; Chatterjee, 2015). A relevant body of literature is also devoted to methods for the selection of the number of communities in stochastic blockmodels (Zhao et al., 2011; Bickel and Sarkar, 2016; Lei, 2016; Newman and Reinert, 2016; Franco Saldaña et al., 2017; Riolo et al., 2017; Wang and Bickel, 2017; Chen and Lei, 2018; Rastelli et al., 2018). Often, practitioners simply set $d=K$, for some $d$, assuming that $\mathbf{B}$ has full rank in the stochastic blockmodel framework. Under the full rank assumption, one may estimate $d=K$ as the number of eigenvalues of $\mathbf{A}$ which are larger than $\sqrt{n}$ (Chatterjee, 2015; Lei, 2016). In this setting, hypothesis tests have also been proposed (Zhao et al., 2011; Bickel and Sarkar, 2016; Lei, 2016). In this work, the problem of selecting $d$ is approached from the perspective of variable selection in model based clustering (Fraley and Raftery, 2002; Lau and Green, 2007), which is widely studied in the literature (Fowlkes et al., 1988; Law et al., 2004; Raftery and Dean, 2006; Maugis et al., 2009; Andrews and McNicholas, 2014).

Similarly, the problem of correctly selecting the number of clusters is also common in $K$-means or GMMs, since it is usually required to specify a number of components in the mixture. Usually the parameter is chosen by minimising information criteria (for example, AIC or BIC). Numerous other techniques have been proposed for GMMs with unknown number of components (Mengersen and Robert, 1996; Richardson and Green, 1997; Stephens, 2000; Nobile, 2004; Zhang et al., 2004; Dellaportas and Papageorgiou, 2006; Nobile and Fearnside, 2007; Miller and Harrison, 2018).

Clearly, the sequential approach in estimating $d$ and $K$ is suboptimal, and it is desirable to jointly estimate the two parameters. This article addresses the problem in a Bayesian framework, proposing a novel methodology to automatically select $d$ and $K$, simultaneously. Techniques for selection of $K$ in GMMs will be incorporated within the spectral embedding framework, allowing for $K$ and $d$, the number of communities and latent dimension of the latent positions, to be random and learned from the data. A structured Bayesian model for simultaneously inferring the dimension of the latent space, the number of communities, and the community allocations is proposed. The model is based on asymptotic results (Tang and Priebe, 2018; Rubin-Delanchy et al., 2017) on the leading and informative components of spectral embeddings, obtained for $d$ fixed and known. The asymptotic theory is combined with realistic assumptions about the remaining components of the embedding, empirically tested and justified on simulated data. Furthermore, extensions to the directed and bipartite case will be discussed. The proposed model has multiple advantages: the latent dimension $d$ and number of communities $K$ are modelled separately, and the Bayesian framework allows for automatic selection of the two parameters. The model also allows estimation of $d$ even when $d<K$, and gives insights on the goodness-of-fit of the stochastic blockmodel on observed network data, based on the embedding. The method is tested on simulated data and applied to real world computer and transportation networks. It should be noted that Yang et al. (2019) have simultaneously and independently proposed a similar inferential procedure within a frequentist framework.

The article is organised as follows: Section 2 introduces adjacency and Laplacian spectral embeddings and the GRDPG. The novel Bayesian model for selection of the appropriate dimension of the latent space is discussed in Section 3, followed by careful illustration of posterior inference procedures in Section 4. Section 5 discusses the effects of the curse of dimensionality on the model, and suggests a remedy. Extensions of the model are presented in Section 6, and results and applications are finally discussed in Section 7.

2 GRDPG and spectral embeddings

The adjacency matrix of an undirected graph is usually embedded into a latent space of dimension $d$ using one of two different procedures: the adjacency spectral embedding or the Laplacian spectral embedding. Suppose $\mathbf{A}\in\{0,1\}^{n\times n}$ is a symmetric adjacency matrix of an undirected graph with $n$ nodes.

Definition 1** (Adjacency spectral embedding – ASE).**

For $d\in\{1,\ldots,n\}$, consider the spectral decomposition

[TABLE]

where $\hat{\bm{\Lambda}}$ is a $d\times d$ diagonal matrix containing the top $d$ eigenvalues in magnitude, in decreasing order, $\hat{\bm{\Gamma}}$ is a $n\times d$ matrix containing the corresponding orthonormal eigenvectors, and the matrices $\hat{\bm{\Lambda}}_{\perp}$ and $\hat{\bm{\Gamma}}_{\perp}$ contain the remaining $n-d$ eigenvalues and eigenvectors. The adjacency spectral embedding $\hat{\mathbf{X}}=[\hat{\bm{x}}_{1},\dots,\hat{\bm{x}}_{n}]^{\top}$ of $\mathbf{A}$ in $\mathbb{R}^{d}$ is

[TABLE]

where the operator $|\cdot|$ applied to a matrix returns the absolute value of its entries.

Definition 2** (Laplacian spectral embedding – LSE).**

For $d\in\{1,\ldots,n\}$, consider the (modified) normalised Laplacian matrix

[TABLE]

and its spectral decomposition

[TABLE]

The Laplacian spectral embedding $\tilde{\mathbf{X}}=[\tilde{\bm{x}}_{1},\dots,\tilde{\bm{x}}_{n}]^{\top}$ of $\mathbf{A}$ in $\mathbb{R}^{d}$ is

[TABLE]

The modified Laplacian $\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}$ (Rohe et al., 2011) is preferred to the more common $\mathbf{I}_{n}-\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}$ since the eigenvalues of the former lie in $(-1,1)$, providing a convenient interpretation for disassortative networks (Rubin-Delanchy et al., 2016). Also, using the Laplacian embedding is preferable in sparse regimes, whereas the adjacency embedding should be used for relatively dense graphs (Cape et al., 2018). A case-study with discussion on the differences between community structures obtained from ASE and LSE is presented in Priebe et al. (2019).

In this work, the stochastic block model will be interpreted as a specific case of a generalised random dot product graph (GRDPG) (Rubin-Delanchy et al., 2017). The GRDPG, a generic latent position model for graphs, is defined below.

Definition 3** (Generalised random dot product graph, GRDPG).**

Let $d_{+},d_{-}$ be non-negative integers such that $d=d_{+}+d_{-}$. Let $\mathbb{X}\subseteq\mathbb{R}^{d}$ such that $\forall\ \bm{x},\bm{x}^{\prime}\in\mathbb{X}$, $0\leq\bm{x}^{\top}{\mathbf{I}}(d_{+},d_{-})\bm{x}^{\prime}\leq 1$, where

[TABLE]

Let $\mathcal{F}$ be a probability measure on $\mathbb{X}$, $\mathbf{A}\in\{0,1\}^{n\times n}$ be a symmetric matrix and $\mathbf{X}=(\bm{x}_{1},\dots,\bm{x}_{n})^{\top}\in\mathbb{X}^{n}$. Then $(\mathbf{A},\mathbf{X})\sim\mathrm{GRDPG}_{d_{+},d_{-}}(\mathcal{F})$ if $\bm{x}_{1},\dots,\bm{x}_{n}\overset{iid}{\sim}\mathcal{F}$ and for $i<j$, independently

[TABLE]

To represent the $K$-community stochastic blockmodel as a GRDPG model, $\mathcal{F}$ can be chosen to have mass concentrated at $\bm{\mu}_{1},\dots,\bm{\mu}_{K}\in\mathbb{R}^{d}$ such that $\bm{\mu}_{i}^{\top}{\mathbf{I}}(d_{+},d_{-})\bm{\mu}_{j}=B_{ij}\ \forall\ i,j\in\{1,\dots,K\}$. For estimation of the latent positions in a SBM, interpreted as a GRDPG, Rubin-Delanchy et al. (2017) suggest the following algorithm, based on Gaussian mixture modelling (see, for example, Fraley and Raftery, 2002).

Intuitively, the algorithm approximately holds because, taking a graph with $m$ nodes, and restricting the attention to the first $n$ nodes, with $n<m$:

[TABLE]

where $\mathbf{Q}_{m}$ is a matrix from the indefinite orthogonal group $\mathbb{O}(d_{+},d_{-})$ and $\bm{\Sigma}(\bm{\mu}_{z_{i}})$ can be analytically computed (Rubin-Delanchy et al., 2017). The result holds for $d$ fixed and known, but in this work it is of interest to treat $d$ as a random, unknown parameter. If a $m$-dimensional embedding is considered, with $m>d$, then asymptotic theory implies an approximate normal distribution with non-zero means and an elliptic covariance within each cluster for the top-$d$ components of the embedding; but, to the best of our knowledge, no theoretical results have been obtained for the remaining $m-d$ columns. It is therefore necessary to propose a model for the remaining part of the embedding, which will be carefully described in Section 3, and empirically justified and assessed in Section 7.1.

It should be noted that it is only possible to estimate the vectors $\{\bm{\mu}_{j}\}$ up to an orthogonal transformation; specifically, for any matrix $\mathbf{Q}\in\mathbb{O}(d_{+},d_{-})$, $(\mathbf{Q}\bm{\mu}_{z_{i}})^{\top}{\mathbf{I}}(d_{+},d_{-})(\mathbf{Q}\bm{\mu}_{z_{j}})=\bm{\mu}_{z_{i}}^{\top}{\mathbf{I}}(d_{+},d_{-})\bm{\mu}_{z_{j}}$, which implies that the likelihood (1.1) is invariant to any such transformation. This is inconsequential for much the same reason, however, as knowledge of the true transformation would not change any inferences about the network.

A second source of non-identifiability in the GRDPG interpretation of the SBM is the "uniqueness up to artificial dimension blow-up" (Cape et al., 2018): for $(\mathbf{A},\mathbf{X})\sim\mathrm{GRDPG}_{d_{+},d_{-}}(\mathcal{F})$, there exists $\check{\mathcal{F}}$ on $\mathbb{R}^{\check{d}}$, with $\check{d}>d$, such that $(\mathbf{A},\mathbf{X})\overset{d}{=}(\check{\mathbf{A}},\check{\mathbf{X}})$ with $(\check{\mathbf{A}},\check{\mathbf{X}})\sim\mathrm{GRDPG}_{\check{d}_{+},\check{d}_{-}}(\check{\mathcal{F}})$. In the stochastic blockmodel setting, this essentially means that any matrix $\mathbf{B}\in[0,1]^{K\times K}$ with rank $d$ can be obtained as an inner product between latent positions on arbitrarily large dimensions.

3 A Bayesian model for SBM embeddings

For simplicity, the embeddings will be generically denoted as ${\mathbf{X}}=[\bm{x}_{1},\dots,\bm{x}_{n}]^{\top}\in\mathbb{R}^{n\times m},\ \bm{x}_{i}\in\mathbb{R}^{m}$ for some $m$, $d\leq m\leq n$. In this article, $m$ is always assumed to be fixed and obtained from a preprocessing step. Choosing an appropriate value of $m$ is arguably much easier than choosing the correct $d$, and, in the proposed model, the correct $d$ can be recovered for any choice of $m$, as long as $d\leq m$. Let $\mathbf{X}_{:d}$ denote the first $d$ columns of ${\mathbf{X}}$, and $\mathbf{X}_{d:}$ the $m-d$ remaining columns. The notation $\bm{x}_{i,:d}$ denotes the first $d$ elements $(x_{1},\dots,x_{d})$ of the vector $\bm{x}_{i}$, and similarly $\bm{x}_{i,d:}$ denotes the last $m-d$ elements $(x_{d+1},\dots,x_{m})$.

Suppose a latent dimension $d$, $K$ communities, and latent community assignments $\bm{z}=(z_{1},\dots,z_{n})$. The latent positions of nodes within each community are assumed to be generated from an $m$-dimensional community-specific Gaussian distribution, where the first $d$ components are modelled differently from the remaining $m-d$: the initial components $\bm{x}_{i,:d}$ are assumed to have unconstrained mean vector $\bm{\mu}_{k}\in\mathbb{R}^{d}$ and positive definite $d\times d$ covariance matrix $\bm{\Sigma}_{k}$; in contrast, for $\bm{x}_{i,d:}$, two constraints are imposed: the mean is an $(m-d)$-dimensional vector of zeros, and the covariance is a diagonal matrix $\bm{\sigma}^{2}_{k}\mathbf{I}_{m-d}$ with positive entries $\bm{\sigma}^{2}_{k}=(\sigma^{2}_{k,d+1},\dots,\sigma^{2}_{k,m})$. The validation of the model assumptions will be discussed in Section 7.1. For mathematical convenience, conjugate priors can be placed on the parameters as follows:

[TABLE]

where if, $n_{k}=\sum\nolimits_{i=1}^{n}\mathds{1}_{k}\{z_{i}\}$ is the size of community $k$, $K_{\varnothing}=\sum_{k=1}^{K}\mathds{1}_{\mathbb{N}_{+}}\{n_{k}\}$ is the number of non-empty communities. Note that the inverse Wishart has been partially re-parametrised using a parameter $\nu_{0}>0$. The dimension of the corresponding matrix is then added to obtain the required constraint $\nu_{0}+d-1>d-1$ for the generic parametrisation and interpretation of the degrees of freedom $\nu$ (in this case $\nu=\nu_{0}+d-1$) of the distribution. Also, note that $m$ can be generically chosen to be equal to $K$, when fixed, or equal to $n$ to have the maximum possible dimension of the embedding. Since one of the specific objectives of the analysis is to learn the number of components of the mixture, rather than density estimation, the widely used infinite Bayesian nonparametric mixtures are not appropriate in this context (Miller and Harrison, 2014).

Yang et al. (2019) have proposed a similar model in a frequentist framework. The conjecture on the distribution of $\mathbf{X}_{d:}$ is essentially the same, except for the choice of the diagonal elements of the cluster-specific covariance matrix: Yang et al. (2019) use a common variance parameter $\sigma^{2}_{k}$ for the last $m-d$ columns of the embedding, whereas a $(m-d)$-dimensional vector of variances $\bm{\sigma}^{2}_{k}$ is used in this paper. Additionally, as a second difference from Yang et al. (2019), the full model proposed here will also incorporate a second-level community cluster structure on these vectors of variances, which will be introduced in Section 5.

A cartoon representation of the model is given in Figure 1. Note that the condition $d\leq K$ is explicitly enforced in (3.7). More specifically, $d\leq K_{\varnothing}$, which avoids an artificial matching between $d$ and $K$ using empty clusters, which are given non-zero probability mass under the Dirichlet-Multinoulli prior on $(\bm{\theta},\bm{z})$. One can also model $d$ and $K$ separately in an analogous way, changing the prior $p(d)$ to, for example,

[TABLE]

independently of $K$ and $\bm{z}$; this will later be referred to as the unconstrained model. The alternative prior (3.8) is particularly useful in practical applications and provides a useful interpretation of $d$: when $d\leq K$, then $d=\mathrm{rank}(\mathbf{B})$, but when $d>K$, this implies that the observed data might deviate from the stochastic blockmodel assumption, and provides a useful diagnostic for model validation and goodness-of-fit.

The likelihood associated with the spectral embedding $\mathbf{X}\in\mathbb{R}^{n\times m}$ obtained from a stochastic blockmodel can be expressed as:

[TABLE]

where $\phi(\cdot)$ denotes the (possibly multivariate) Gaussian density function. Hence, the posterior, up to a normalising constant, has form

[TABLE]

The $\mathrm{NIW}_{d}(\bm{0},\kappa_{0},\nu_{0}+d-1,\bm{\Delta}_{d})$ prior on the pair $(\bm{\mu}_{k},\bm{\Sigma}_{k})$ is conjugate and yields a conditional posterior $(\bm{\mu}_{k},\bm{\Sigma}_{k})|\mathbf{X},\bm{z},d\overset{d}{\sim}\mathrm{NIW}_{d}(\bm{m}^{(k)}_{:d},\kappa_{n_{k}},\nu_{n_{k}}+d-1,\mathbf{D}^{(k)}_{:d})$. By standard methods for inference in a multivariate Gaussian mixture model with $\mathrm{NIW}$ prior, the covariance matrix $\bm{\Sigma}_{k}$ can be explicitly integrated out from the posterior to obtain

[TABLE]

the density of the multivariate Student $t$ distribution with $\nu_{n_{k}}$ degrees of freedom, where

[TABLE]

Henceforth, $\bm{\mu}_{k}$ can easily be resampled in a simple Gibbs sampling step, conditional on the actual value of $d$ and on the community allocations $\bm{z}$. In this work, the location vectors $\bm{\mu}_{k}$ are collapsed out too, but the distribution is instructive to present other distributional results below, and could be also used if the objective of the analysis is to also recover the explicit form of the latent positions.

In a multivariate Gaussian model with normal inverse Wishart prior, it is also possible to analytically express the marginal likelihood of the observed data. Here, conditioning on a community-specific Gaussian, on the assignments $\bm{z}$ and on the dimension of the latent space $d$:

[TABLE]

where $\mathbf{X}^{(k)}_{:d}$ is the subset of rows of $\mathbf{X}_{:d}$ such that $z_{i}=k$.

Given the $\mathrm{Inv}\text{-}\chi^{2}(\lambda_{0},\sigma^{2}_{0})$ prior, the posterior for $\sigma^{2}_{j,k}$ is $\mathrm{Inv}\text{-}\chi^{2}(\lambda_{n_{k}},s^{(k)}_{j})$, where

[TABLE]

Similar calculations give the full marginal likelihood for the remaining portion of the embedding $\mathbf{X}^{(k)}_{d:}$:

[TABLE]

Finally, if $d$ is considered as a nuisance parameter:

[TABLE]

which is easily computed using (3.16) and (3.18).

Also, note that the probabilities $\bm{\theta}$ associated to the community assignment can be easily integrated out, resulting in the following marginal likelihood, conditional on $K$:

[TABLE]

The distributional results presented in (3.15), (3.16), and (3.20) (for a proof, see, for example, Murphy, 2007) are the building blocks for the MCMC sampler which is used to make Bayesian inference on the model parameters of interest.

4 Inference via MCMC

Since the full posterior is not analytically tractable, inference is performed using MCMC sampling with trans-dimensional moves (Green, 1995). The main objective of the analysis is to cluster the nodes, and therefore the locations $\bm{\mu}_{k}$, the variance parameters $\bm{\Sigma}_{k}$ and $\bm{\sigma}^{2}_{k}$ and the community probabilities $\bm{\theta}$ are considered as nuisance parameters and integrated out. Essentially, in this type of collapsed Gibbs sampler (Liu, 1994), four moves are available (Richardson and Green, 1997; Zhang et al., 2004), described in the subsequent four subsections.

4.1 Change in the community allocations

A fully collapsed Gibbs update for each community assignment is available:

[TABLE]

In the special case where $d=K_{\varnothing}$ and $n_{z_{i}}=1$, the full conditional distribution for $z_{i}$ assigns probability one to retaining the same value since the model does not permit $d>K_{\varnothing}$. Otherwise, from (3.20):

[TABLE]

where $n_{k}^{-i}=n_{k}-\mathds{1}_{k}(z_{i})$. Similarly, the remaining term in (4.1), $p({\bm{x}}_{i}|\{{\bm{x}}_{j}\}_{j\neq i:z_{j}=k},d)$, can be obtained as the ratio of marginal likelihoods

[TABLE]

The ratio (4.3) can be decomposed as the product of two ratios of marginal likelihoods. Using (3.16), the first ratio is equivalent to the following multivariate Student $t$ distribution (Murphy, 2007):

[TABLE]

where the additional superscript $-i$ denotes a cluster quantity that is computed excluding the allocation $z_{i}$ of the $i$-th node. The second ratio, which accounts for the last $m-d$ dimensions, has the form

[TABLE]

4.2 Split or merge two communities

To vary the number of communities, move proposals inspired by Sequentially-Allocated Merge-Split sampling (Dahl, 2003) are used here, but other alternative choices are available (see, for example, Jain and Neal, 2004; Bouchard-Côté et al., 2017). Two indices $i$ and $j$ are sampled at random from the $n$ nodes, and without loss of generality assume $z_{i}\leq z_{j}$. If $z_{i}=z_{j}$, then the single cluster is split, whereas if $z_{i}>z_{j}$ the two clusters are merged. In both move types, node $i$ will remain in the same cluster, denoted $z^{\star}_{i}=z_{i}$. In the merge move, all elements of cluster $z_{j}$ are reassigned to cluster $z_{i}$ (with any higher indexed clusters subsequently decremented). For the split move, node $j$ is first reassigned to cluster $K^{\star}=K+1$ with new allocation $z^{\star}_{j}=K^{\star}$; then, in random order the remaining nodes currently allocated to cluster $z_{i}$ are randomly reassigned to clusters $z_{i}$ or $K^{\star}$ with probability proportional to their predictive distribution from the generative model (4.3). Denoting the resulting product of renormalised predictive densities from these reallocations by $q(K^{\star},\bm{z}^{\star}|K,\bm{z})$, the acceptance probability for a split move, for example, is

[TABLE]

Note that the ratio of densities for $\mathbf{X}$ in the acceptance ratio (4.6) will depend only upon the rows of the matrix corresponding to the cluster being split (or similarly, merged), and furthemore these expressions will decompose as a products of terms for the first $d$ and remaning $m-d$ components (cf. (4.4), (4.5)).

4.3 Create or remove an empty community

Adding or removing empty communities whilst fixing $\bm{z}$ corresponds to proposing $K^{\star}=K+1$ or $K^{\star}=K-1$ respectively, although the latter proposal is not possible if $K=K_{\varnothing}$, meaning there are currently no empty communities. The acceptance probability is simply

[TABLE]

where the proposal ratio $q_{\varnothing}=2$ if $K^{\star}=K_{\varnothing}$, $q_{\varnothing}=0.5$ if $K=K_{\varnothing}$ and $q_{\varnothing}=1$ otherwise.

4.4 Change in the latent dimension

This move is only required when the latent dimension is not marginalised out. Given a current value $d$, a new value $d^{\star}$ is proposed from a density $q(d^{\star}|d)\propto\xi^{|{d^{\star}-d}|}\mathds{1}_{\mathcal{D}}(d^{\star})$ on a neighbourhood

[TABLE]

typically with $l\leq 5$ and with $\xi\in(0,1)$. The acceptance ratio reduces to

[TABLE]

Notably, if $d^{\star}>d$, the ratio $p(\mathbf{X}|d^{\star},K,\bm{z})/p(\mathbf{X}|d,K,\bm{z})$ only depends on the first $d^{\star}$ components of the embedding, since the last $m-d^{\star}$ components remain independent by (3.7).

4.5 Inferring communities

Markov Chain Monte Carlo samplers for mixture models with varying number of clusters are well known to be affected by label switching (Jasra et al., 2005), since the likelihood is invariant to permutations of the cluster labels. However, the estimated posterior similarity between nodes $i$ and $j$, $\hat{\pi}_{ij}=\hat{\mathbb{P}}(z_{i}=z_{j}|\mathbf{X})=\sum_{s=1}^{M}\mathds{1}_{z_{i}^{(s)}}\{z_{j}^{(s)}\}/M$ is invariant to label switching. Communities can be estimated from the MCMC chains using the posterior similarity matrix and the PEAR method (maximisation of the posterior expected Rand adjusted index, Fritsch and Ickstadt, 2009, R package mcclust). Alternatively, if a configuration with a fixed number of clusters $K$ is required, the clusters can been estimated using hierarchical clustering with average linkage, using $1-\hat{\pi}_{ij}$ as distance measure (Medvedovic et al., 2004). Many alternatives to this method have been proposed in the literature (Rand, 1971; Binder, 1978; Dahl, 2006).

5 Second-level clustering of community variances

Empirical analysis of assuming the model (3.7) for simulations from the stochastic block model show that identifying and clearly separating the $K$ clusters in $\mathbf{X}_{d:}$ is particularly difficult for the sampler in settings when $d\ll m$. The problem is particularly evident when $m=n$ and $d$ is small. In this case, it has been assessed empirically that the within-cluster variance of the true communities in simulated datasets seems to converge to similar values, such that $\sigma^{2}_{k,j}\approx\sigma^{2}_{\ell,j}$ for $j\gg d$ and $k\neq\ell$. Therefore, when $m$ is large enough, the selected model tends to be under-specified: the correct dimension $d$ is identified, but the true number of communities $K$ is underestimated. This is also one of the main reasons why it is not advisable to directly fit a Gaussian mixture model on $\mathbf{X}\in\mathbb{R}^{n\times m}$ and allow $K$ to be random, ignoring the role of $d$.

The problem is illustrated in Figure 2, which shows the results from performing ASE for a simulation of $n=500$ nodes from a stochastic block model containing five equally probable communities with well-separated mean locations. The within-cluster variance of three of the five communities fluctuate around the same values on each dimension larger than $d$. For a dimension larger than approximately $150$, four of the five clusters have the same variance on the subsequent dimensions. Therefore, when $m=n$, the MCMC sampler selects the MAP estimate $\hat{K}=2$ for parsimony, and increases the variance of the Gaussian distributions on the first two dimensions, on which the clusters are well separated.

The solution proposed here is to assume shared variance parameters between some of the clusters for dimensions larger than $d$. Specifically, each community $k\in\{1,\dots,K\}$ is assigned a second-level cluster allocation $v_{k}\in\{1,\dots,H\}$, with $H\leq K$. If $v_{k}=v_{\ell}$, then for $j>d$, $\sigma^{2}_{k,j}=\sigma^{2}_{\ell,j}$. Formally,

[TABLE]

Essentially, the vector $\bm{v}=(v_{1},\dots,v_{K})$ defines a clustering of communities. Figure 3 depicts a cartoon example of this extended model. Note that if $H=1$ and all the communities are assigned to the same second-level cluster, the problem of selecting $d$ essentially reduces to an ordinal version of the feature selection problem in clustering (Raftery and Dean, 2006; Maugis et al., 2009).

Under this extended model, the posterior distribution for $\sigma^{2}_{j,k}$ changes due to the aggregation of communities in the second level. Under the $\mathrm{Inv}\text{-}\chi^{2}(\lambda_{0},\sigma^{2}_{0})$ prior, the posterior is $\mathrm{Inv}\text{-}\chi^{2}(\lambda_{n_{\bullet k}},s^{(\bullet k)}_{j})$ where

[TABLE]

Calculations similar to (3.18) give the correct form of the marginal likelihood for the right hand side of the matrix, restricted to a given value of $v_{k}$. Clearly, $\bm{\phi}$ can be again marginalised out, yielding the marginal likelihood

[TABLE]

The MCMC sampler described in Section 4 must be slightly adapted. For the Gibbs sampling move in Section 4.1, the product of univariate Student’s $t$ densities in (4.5) is modified using the appropriate $(\lambda_{n_{\bullet k}},s^{(\bullet k)}_{j})$ pair. For the change in dimension, the likelihood $p(\mathbf{X}|d,K,\bm{z},\bm{v})$ should be computed using the shared variances and the allocations $\bm{v}$. When an empty community is proposed, as in Section 4.3, the ratio $p(\bm{v}^{\star}|K)/p(\bm{v}|K)$ must be added, limited to the second level allocation of the additional community. The value $v_{k}$ for the proposed empty cluster can be simply chosen at random from $\{1,\dots,H\}$. Finally, for the split-merge move in Section 4.2, if $z_{i}=z_{j}$ for the two selected nodes, then $v_{z_{i}}=v_{z_{j}}$ after the split move. Alternatively, if $z_{i}\neq z_{j}$, then the new value of $v_{k}$ is sampled at random from $v_{z_{i}}$ and $v_{z_{j}}$.

Finally, three additional moves are required: resampling the second-level cluster allocations $\bm{v}$ using a Gibbs sampling step; proposing a second-level split-merge move; and adding or removing an empty second-level cluster. When $\bm{\phi}$ and the parameters of the Gaussian distributions are marginalised out, the second-level allocations are resampled according to the following equation:

[TABLE]

where the independence assumption between $\mathbf{X}^{(k)}_{:d}$ and $\mathbf{X}^{(k)}_{d:}$ is used. Similarly to (4.2):

[TABLE]

The calculations for the second term in (5.7) are similar to (4.3):

[TABLE]

which can be computed using (3.18). The second-level split-merge move and the proposal of an empty cluster follows the same guidelines in Sections 4.2 and 4.3.

Potentially, the model could be extended further using the same reasoning: from the plot in Figure 2, it is clear that the different clusters begin to share the same variance at different points in the plot. Empirically, all the variances approximately converge to the same values at large dimensions, and it is therefore possible to identify a $(K-1)$-vector of discrete points in $\{d,d+1,\dots,m\}$ at which different community variances coalesce. For the plot in Figure 3, such vector could be $(d,d,d,150,n)$, with $d=2$ and $n=m=500$.

6 Extension to directed and bipartite graphs

A directed graph $\mathbb{G}=(V,E)$ has the property that $(i,j)\in E\centernot\implies(j,i)\in E$, meaning the corresponding adjacency matrix $\mathbf{A}\in\{0,1\}^{n\times n}$ is not, in general, symmetric. Directed graphs are useful for representing directed interaction networks, such as email traffic patterns; knowing that individual $i$ broadcasts emails to individual $j$ does not immediately imply that $j$ also issues communications to $i$.

Definition 4** (Adjacency embedding of the directed graph).**

Given a directed graph with adjacency matrix $\mathbf{A}\in\{0,1\}^{n\times n}$, and a positive integer $d$, $1\leq d\leq n$, consider the singular value decomposition

[TABLE]

where $\hat{\mathbf{D}}\in\mathbb{R}_{+}^{d\times d}$ is diagonal matrix containing the top $d$ singular values in decreasing order, $\hat{\mathbf{U}}\in\mathbb{R}^{n\times d}$ and $\hat{\mathbf{V}}\in\mathbb{R}^{n\times d}$ contain the corresponding left and right singular vectors, and the matrices $\hat{\mathbf{D}}_{\perp}$, $\hat{\mathbf{U}}_{\perp}$, and $\hat{\mathbf{V}}_{\perp}$ contain the remaining $n-d$ singular values and vectors. The $d$-dimensional directed adjacency embedding of $\mathbf{A}$ in $\mathbb{R}^{d}$, is defined as the pair

[TABLE]

Writing $\mathbf{X}={\mathbf{U}}{\mathbf{D}}^{1/2}$ and $\mathbf{X}^{\prime}={\mathbf{V}}{\mathbf{D}}^{1/2}$, the rows of $\mathbf{X}$ characterise the activities of each node as a source, and the rows of $\mathbf{X}^{\prime}$ characterise the same nodes as destinations.

The model in (3.7) can be easily adapted to directed graphs. Treating the embeddings $\mathbf{X}$ and $\mathbf{X}^{\prime}$ as independent, each is modelled separately using the same Gaussian structure and prior distributions (3.7), except for three parameters which are initially assumed common to both embeddings: the latent dimension $d$, the number of communities $K$ and the vector of node assignments to those communities, $\bm{z}$.

In some contexts it will be more relevant to consider different community membership structures for the same set of nodes when considering them as source nodes or destination nodes. In this case, let $K$ denote the number of source communities and $K^{\prime}$ denote the number of destination communities; similarly let $\bm{z}$ denote the assignments of nodes to source communities, and $\bm{z}^{\prime}$ the allocations to destination communities. Jointly learning $\bm{z}$ and $\bm{z}^{\prime}$ (as well as $d$) is a problem commonly known as co-clustering, and the corresponding network model is known as the stochastic co-blockmodel (ScBM) (Rohe et al., 2016). Given an asymmetric matrix $\mathbf{B}\in[0,1]^{K\times K^{\prime}}$, then $\mathbb{P}(A_{ij}=1)=B_{z_{i}z^{\prime}_{j}}$. From a random dot product graph perspective, it is assumed that $B_{z_{i}z^{\prime}_{j}}=\bm{\mu}_{z_{i}}^{\top}\bm{\mu}^{\prime}_{z^{\prime}_{i}}$, for some latent positions $\bm{\mu}_{z_{i}},\bm{\mu}^{\prime}_{z^{\prime}_{i}}\in\mathbb{R}^{d}$ and $d=\operatorname{rank}(\mathbf{B})\leq\min(K,K^{\prime})$.

The Bayesian model for ScBMs can be easily represented as a separate model for $\mathbf{X}$ and $\mathbf{X}^{\prime}$, of the form given in (3.7), with the latent dimension of the embedding $d$ now the only common parameter. Inference via MCMC can be performed in an equivalent way to the method described in Section 4; the only difference is in the expression of the acceptance ratio for a change in the shared latent dimension $d$, but the procedure can exploit the results obtained in Section 4.4, using the fact that

[TABLE]

where all the marginal likelihoods can be equivalently obtained from (3.16). Furthermore, the model can be appropriately modified when $d\ll m$ to include the second-level cluster allocations proposed in Section 5.

Finally, in bipartite graphs, the observed nodes can be partitioned into two sets $V$ and $V^{\prime}$, with $V\cap V^{\prime}=\varnothing$ and $E\cap(V\times V\cup V^{\prime}\times V^{\prime})=\varnothing$. Assume that $V$ plays the role of the set of source nodes and $V^{\prime}$ of the set of destination nodes. Bipartite graphs are usually represented by a rectangular bi-adjacency matrix $\mathbf{A}\in\{0,1\}^{n\times n^{\prime}}$, with $n^{\prime}=|{V^{\prime}}|$. In this case, it is still possible to apply the methods described in this section to the SVD embedding obtained from the rectangular matrix $\mathbf{A}$. Note that the ScBM extends trivially to the bipartite graph case, which is essentially a special case of a directed graph, with the cluster configurations for source and destination nodes now inescapably unrelated and each node possessing only one latent representation in $\mathbb{R}^{d}$.

7 Applications and results

The Bayesian latent feature network models described in this article have been applied to a both simulated and real world network data from undirected, directed and bipartite graphs. The real network data analysed are from an undirected network obtained from the Santander bikes hires in London, and the Enron Email Dataset; details are given in the corresponding sections.

The model and MCMC sampler have been tested using different combinations of the hyperparameters, showing robustness to the prior choice. Inferential performance is sensitive to extreme values of the variance parameters, relative to their prior mean, but otherwise robust. So in practice, the expectation of the prior for the variance parameters should be chosen to be on the same scale as the observed data. The default settings for the MCMC sampler used in the next sections are the usual uniformative values $\kappa_{0}=\nu_{0}=\lambda_{0}=\alpha=\beta=1$, and $\omega=\delta=0.1$. For the proposal of change in dimension (cf. Section 4.4), $\xi=0.8$. The algorithms were run for a total of $M=\numprint{500000}$ samples with burn-in $\numprint{25000}$, for a number of different chains to check for convergence. The cluster configuration has been initialised using $K$-means for some pre-specified $K$, usually chosen according to the scree-plot criterion. The second-level clusters have been initialised setting $H=K$. In order to set the prior covariances to a realistic value, the correlations in the $\bm{\Delta}_{d}$ matrices are set to zero, and the elements on the diagonal of $\bm{\Delta}_{d}$ to the average within-cluster variance based on the $K$-means cluster configuration. Similarly, the prior $\sigma^{2}_{0j}$ values have been set to the total variance on the corresponding column of the embedding.

7.1 Synthetic data and model validation

A stochastic blockmodel can be simulated starting from a matrix $\mathbf{B}\in[0,1]^{K\times K^{\prime}}$ containing the probabilities of connection between communities, and a vector $\bm{\theta}$ of community allocation probabilities. Here, each element $B_{k\ell}$ of $\mathbf{B}$ was generated from a $\mathrm{Beta}(1.2,1.2)$ distribution, which produces communities with a moderate level of separation. For an undirected graph, $K=K^{\prime}$ and the constraint $B_{k\ell}=B_{\ell k}$ is imposed; similarly, for directed graphs with a shared cluster configuration (cf. Section 6), $K=K^{\prime}$.

The latent dimension $d$ can be interpreted as the rank of the matrix $\mathbf{B}$, and a random matrix $\mathbf{B}$ generated from independent beta draws has full rank with probability $1$. Therefore, to simulate $d<K$ a low-rank approximation of $\mathbf{B}$ must be used to generate the embedding. For undirected graphs, a truncated spectral decomposition can be used: $\tilde{\mathbf{B}}=\bm{\Gamma}_{d}\bm{\Lambda}_{d}\bm{\Gamma}_{d}^{\top}$ (recall Definition 1). Similarly, for the directed and bipartite graphs, the truncated SVD is an appropriate approximation: $\tilde{\mathbf{B}}=\mathbf{U}_{d}\mathbf{D}_{d}\mathbf{V}_{d}^{\top}$ (see Definition 4). Note that under this low-rank approximation, it must be checked that each element $\tilde{B}_{k\ell}\in[0,1]$.

Figure 4 shows results for a synthetic undirected graph with $d=2$ and $K=5$. The scatterplot of the first two columns of the adjacency embedding $\mathbf{X}$, plotted in Figure 4(a), show well-separated clusters, which can be suitably modelled using a Gaussian mixture. Figure 4(b), shows the next two dimensions. Clearly, the community mean locations are significantly different from zero in just the first two dimensions, and this is further illustrated in Figure 4(c). Similarly, in Figure 4(d), the difference between the within-cluster and overall variance is evident only in the first two dimensions, after which the quantities are of the same order of magnitude. The within-cluster variances differ across communities, suggesting that it is appropriate to have cluster-specific values of $\sigma^{2}_{j,k}$ for $j>d$; this phenomenon can also be witnessed in Figure 4(b). Nevertheless, it also seems appropriate to use a second-level clustering with $H=3$, since the variances of three of the five communities are approximately the same for dimensions larger than $d=2$. Also, in Figure 4(e), the within-cluster correlations between $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$ suggest dependence for at least one of the clusters. On the other hand, the sample within-cluster correlations for $\mathbf{X}_{d:}$ tend to be small and centred around 0, suggesting that the assumption of independence is appropriate in that part of the model (3.7). Finally, form Figure 4(f) the marginal likelihood strongly favours the true value $d=2$, resulting in a posterior distribution essentially consisting of a point mass at the true value.

Figure 5 shows results for a simulated bipartite graph with separate community structures for nodes as sources and destinations, with $d=2$, $K=5$ and $K^{\prime}=3$. Again, the clusters in Figure 5(a) and 5(b) are well-separated and can be easily estimated using the Gaussian mixture model, both for sources and receivers. From Figure 5(c), the zero-mean assumption for the columns with index larger than $d$ seems to hold even for a relatively small number of nodes per community. A similar plot can be produced for $\mathbf{X}^{\prime}$, showing a similar pattern. In Figure 5(d), the marginal likelihood strongly favours the true value $d=2$, which again results in a point mass posterior centred at the true value. Similar considerations hold for variances and correlations, with results which are similar to the plots in Figure 4(d) and 4(e) for the undirected graph.

Overall, for synthetic data, the model seems robust and able to detect the correct $d$ and $K$ in a variety of different settings.

7.2 Undirected graphs: Santander bikes

The Santander Cycle hire scheme is a bike sharing system implemented in central London. Transport for London periodically releases data on the bike hires111The data are available at the following URL: https://cycling.data.tfl.gov.uk/.. Considering this as a network, the nodes correspond to bike sharing stations, and a directed edge between stations $i$ and $j$ is drawn if at least one ride between station $i$ and $j$ is completed within the time period considered. In this example, one week of data were considered, from 5 September until 11 September, 2018. The total number of stations used, $n=783$; the total number of undirected edges $|E|=\numprint{96060}$, implying the adjacency matrix is fairly dense.

The results of the Bayesian inferential procedure, using the unconstrained prior (3.8) for $d$, applied to the adjacency and Laplacian embeddings for the Santander bike network are presented in Figure 6. The initial value of $K$ was set to $10$, with $m=25$, but similar estimates were obtained using different starting points for $K$ and different values of $m$. It is interesting to note the different shapes of the posterior barplots of $K_{\varnothing}$ and $H_{\varnothing}$, Figures 6(a) and 6(b), showing that the second-level clustering is crucial to obtain a more accurate model fit when the adjacency embedding is used. On the other hand, for the Laplacian embedding, the posteriors for $K_{\varnothing}$ and $H_{\varnothing}$ are extremely similar, suggesting that the second-level clustering is not required for $m=25$. The maximum a posteriori (MAP) values $d=11$ (adjacency) and $d=12$ (Laplacian) clearly correspond to the elbows in the corresponding scree-plots in Figures 6(c) and 6(d).

Note that, especially in the case of the adjacency embedding, $d$ and $K$ have similar values, showing that the two graphs might be well described by a stochastic blockmodel. Similarly, the constrained model with $d\overset{d}{\sim}\mathrm{Uniform}\{1,\dots,K_{\varnothing}\}$ (3.7) returns the same MAP estimates for $d$, but the constraint $d\leq K_{\varnothing}$ results in a larger number of small clusters; the posterior of $K_{\varnothing}$ essentially reduces to the rescaled probability mass function obtained from the unrestricted model, constrained such that $d\leq K_{\varnothing}$, since the posterior for $d$ is approximately a point mass.

The resulting estimated clustering for the unconstrained model (3.8) based on the adjacency embedding and $K=11$ (the MAP for $d$), plotted in Figure 7, shows a clear structure: the largest communities have approximately the same extension, with a few exceptions. This is expected, since the bikes are free for the first $30$ minutes and have limited speed, and are therefore used for small distance journeys. Two clusters are significantly smaller than the others, and correspond to touristic areas around Westminster, Covent Garden and Buckingham Palace. On the other hand, two clusters have a large geographical extension, and cover the East and West London areas. For the adjacency embedding, the MAP clustering obtained from the restricted model is almost identical. The PEAR method (Fritsch and Ickstadt, 2009) suggests $K=7$ communities instead. Similarly, if the Laplacian embedding is used, the MAP clustering structure suggested by PEAR has $K=7$ communities for the unconstrained model (3.8) for $d$ and $K=12$ for the constrained model (3.7).

7.3 Directed graphs: Enron Email Dataset

Next, the algorithm is applied to a directed network: the Enron Email Dataset222The entire version of the data is available at the following URL: https://www.cs.cmu.edu/~./enron/.. The Enron database is a corpus of emails sent by the employees of the Enron corporation. The version of the Enron data which has been analysed here is described in Priebe et al. (2005), and consists of $n=184$ nodes and $\numprint{3010}$ directed edges. A directed edge $i\to j$ is drawn if the employee $i$ sent an email to the employee $j$.

The results of analysing this network are presented in Figure 8. The initial value of $K$ was set to $10$, with $m=25$, but again similar results were obtained using different starting points for $K$ and different values of $m$. The plots in Figures 8(a) and 8(b) report the estimated posterior distributions of $K_{\varnothing}$ and $H_{\varnothing}$ for the constrained (3.7) and unconstrained (3.8) models. Interestingly, the MAP estimate for $d$ coincides with the MAP estimate for $K$ in the unconstrained model, which is promising. For the constrained model, the MAP for $K$ exceeds the MAP for $d$ by 1, allowing for $\mathrm{rank}(\mathbf{B})<K$. Overall, $d$ and $K$ have similar values, showing that the graph might be well described by a directed stochastic blockmodel.

The posterior distributions for $K_{\varnothing}$ and $H_{\varnothing}$ in Figures 8(a) and 8(b) are fairly different, showing that the second-level clustering might be relevant for this model. Inference on the model without second-level clustering confirms this impression: the posteriors for $K_{\varnothing}$, presented in Figures 8(c) and 8(d) have a more symmetric shape, and the MAP latent dimension is $d=6$. As before, the MAP for $K$ is $d+1=7$, providing some evidence for the possibility $\mathrm{rank}(\mathbf{B})<K$.

From Figure 8(e), the selected MAP values $d=6$ and $d=9$ for the models with and without second-level clustering seem to be a tradeoff between the two most popular criteria for selection of the appropriate latent dimension: the eigengap heuristic suggests $d=5$ if the second largest difference is considered, and the elbow in the scree-plot is approximately located around $d\approx 15$.

8 Conclusion

In this article, a novel Bayesian model has been proposed for automatic and simultaneous estimation of the number of communities and latent dimension of stochastic blockmodels, interpreted as special cases of generalised random dot product graphs. The Bayesian framework allows the number of communities $K$ and latent dimension $d$ to be treated as random variables, with associated posterior distributions. The postulated model is based on asymptotic results in the theory of network embeddings and random dot product graphs, and has been validated on synthetic datasets, showing good performance at recovering the latent parameters and communities. The model has been extended to directed and bipartite graphs, using SVD embeddings and allowing for co-clustering.

Overall, the main advantages of the proposed methodology is to allow for an arbitrarily large value of $m$, the number of columns (dimension) of the embedding at the first stage of the analysis, and then to treat $d$ and $K$ separately, allowing for the case $d=\mathrm{rank}(\mathbf{B})<K$, which is often overlooked in the literature. Problems arising from overspercification of $m$ are tackled using a second-level clustering procedure. Also, the model provides an automated procedure and criterion to select the dimension of the embedding and an appropriate number of communities. If $d$ is not constrained to be less than or equal to $K$, the model also provides empirical evidence of the goodness-of-fit of a stochastic blockmodel for the observed data. Results on real world network data sets show encouraging results in recovering the correct $d$, when compared to commonly used heuristic methods, and the community structure.

Code and datasets

The python code and datasets are available at https://www.github.com/fraspass/sbm.

Acknowledgements

The authors gratefully acknowledge funding from the EPSRC and the Heilbronn Institute for Mathematical Research.