Optimality of variational inference for stochastic block model with missing links
Year: 2022

Abstract:
Variational methods are extremely popular in the analysis of network data. Statistical guarantees obtained for these methods typically provide asymptotic normality for the problem of estimation of global model parameters under the stochastic block model. In the present work, we consider the case of networks with missing links that is important in application and show that the variational approximation to the maximum likelihood estimator converges at the minimax rate. This provides the first minimax optimal and tractable estimator for the problem of parameter estimation for the stochastic block model with missing links. We complement our results with numerical studies of simulated and real networks, which confirm the advantages of this estimator over current methods.

heading: 1 Introduction
text: The analysis of network data poses both computational and theoretical challenges. Most results obtained in the literature concentrate on the stochastic block model (SBM) which is known to be a good proxy for more general models, such as the inhomogeneous random graph model, [35]. Recently, variational methods ([28, 49]) have attracted considerable attention as they offer computationally tractable algorithms often combined with theoretical guarantees. Theoretical results that one can find for such variational methods provide asymptotic normality rates for parameter estimates of stochastic block data. For example, consistency has been shown for profile likelihood maximization [8] and variational approximation to the maximum likelihood estimator [13], [7]. These results have been extended to the case of dynamic stochastic block model [34] and sampled data [48]. These work focus on parameter estimation, as in [44] and [53], who establish the minimax optimality of variational methods in a large class of models (which does however not include the stochastic block model). Variational inference has also been successfully applied to the problem of community detection, see, e.g., [3, 54, 26, 45]. In particular, the authors of [54] show that an iterative Batch Coordinate Ascent Variational Inference algorithm designed for the two-parameters, assortative stochastic block model achieves statistical optimality for community detection problem. Note that this algorithm cannot be extended to the more general stochastic block model considered here.
In parallel with this line of work, the problem of statistical estimation of model parameters, in particular, the question of minimax optimal convergence rates, has been actively studied in the statistical community. In the case of dense graphs, a pioneering paper [17] shows that, for the problem of estimating the matrix of connection probabilities, the least square estimator is minimax optimal and [18] provides optimal rate for Bayes estimation. For the more challenging case of
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
sparse graphs, the minimax optimal rates have been first obtained in [29] building on the restricted least square estimator. In [16], the authors consider the least square estimator in the setting when observations about the presence or absence of an edge are missing independently at random with the same probability p. Unfortunately, least square estimation is too computationally expensive to be used in practice. Many other approaches have been proposed, for example, spectral clustering [39, 22, 46], modularity maximization [42, 8], belief propagation [14], neighbourhood smoothing [55], convex relaxation of k-means clustering [20] and of likelihood maximization [4], and universal singular value thresholding [11, 30, 51]. These approaches are computationally tractable but show sub-optimal statistical performances. So the question of possible computational gap when no polynomial time algorithm can achieve minimax optimal rate of convergence has been raised.
The present work goes in these two directions. We study the statistical properties of the mean field variational Bayes method and show that it achieves the optimal statistical accuracy. In particular, these results close the open question on the possible existence of a computational gap for the problem of global parameter estimation. We built our analysis on the approach developed in [13], [7] and [48] using the closeness of maximum likelihood and maximum variational likelihood and on the results that show the minimax optimality of the maximum likelihood estimator [19].
In the present paper, we deal with settings where the network is not fully observed, a common problem when studying real life networks. In many applications the network has missing data as detecting interactions can require significant experimental effort, see, [32, 52, 24, 21]. For example, in biology graphs are used to model interactions between proteins. Discovery of these interactions can be costly and time-consuming [9]. On the other hand, the size of some networks from social media or genome sequencing may be so large that only subsamples of the data are considered [6]. It has been observed that incomplete observation of the network structure may considerably affect the accuracy of inference methods [31] and missing data must be taken into account while analyzing networks data. A popular approach consists in considering the edges with uncertain status as non-existing. In the present paper, we use a different framework by considering such edges as missing and introducing a separated data missing mechanism. A natural application of our method is link prediction [36, 56], the task of predicting whether two nodes in a network are connected. Our approach allows to deduce the pairs of nodes that are most likely to interact based on the known interactions in the network. Behind inference of the networks structure, our algorithms can be used to predict the links that may appear in the future if we consider networks evolving over the time. For example, in a social network, two users that are not yet connected but are likely to be connected can be recommended as promising friends.

heading: 1.1 Contribution and outline
text: The paper is organized as follows. After summarizing notations, we introduce our model and the maximum likelihood estimator for the stochastic block model with missing observations in Section 2. In Section 3, we introduce the mean field variational Bayes method and present a new estimator which combines the labels obtained using the variational method and the empirical mean for estimation of connection probabilities. In Section 3.2, we show that our estimator is minimax optimal for dense stochastic block models with missing observations as well as for sparse stochastic block models. Finally, in Section 4 we provide an extensive numerical study both on synthetic and real-life data which shows clear advantages of our estimator over current methods.

heading: 1.2 Notations
text: We provide here a summary of the notations used throughout the paper. For all d 2 N⇤, we denote by [d] the set {1, ..., d}. For z : [k] ! [n] and all (a, b) 2 [d]⇥ [d], we abuse notations and denote z 1(a, b) = {(i, j) : z(i) = a, z(j) = b, i 6= j}. For any two label functions z, z0, we write z ⇠ z0 if there exists a permutation of {1, ..., k} such that (z( (a)))ak = (z(a))ak. For any set S, we denote by |S| its cardinality. For any matrix A, we denote by Aij its entry on row i and column j. If A 2 [0, 1]n⇥n and A is symmetric, we write A 2 [0, 1]n⇥nsym . We denote by A B the Hadamard product of two matrices A and B. The Frobenius norm of a matrix A is denoted by kAk2 = qP i,j A 2 ij . We denote by C and C
0 positive constants that can vary from line to line. These are absolute constants unless otherwise mentioned. For any two positive sequences (an)n2N, (bn)n2N, we write an = !(bn) if an/bn ! 1.

heading: 2 Maximum likelihood estimation in the stochastic block model with missing links

heading: 2.1 Network model and missing data scheme
text: In the simplest situation, a network can be represented as undirected, unweighted graph with n nodes indexed from 1 to n. Then, the network can be encoded by its adjacency matrix A = (Aij). The adjacency matrix is a n⇥n symmetric matrix such that for any i < j, Aij = 1 if there exists an edge between node i and node j, Aij = 0 otherwise. We consider that there is no edge linking a node to itself, so Aii = 0 for any i. A common approach in network data analysis is to assume that the observations are random variables drawn from a probability distribution over the space of adjacency matrices. More precisely, for i < j the variables Aij are assumed to be independent Bernoulli random variables of parameter ⇥⇤ij , where ⇥
⇤ = (⇥⇤ij)1i<jn is a n⇥ n symmetric matrix with zero diagonal entries. The matrix ⇥⇤ corresponds to the matrix of probabilities of observing an edge between nodes i and j. This model is known as the inhomogeneous random graph model:
81  i < j  n, Aij |⇥⇤ij ind.⇠ Bernoulli ⇥⇤ij . (1)
Our focus is on the problem of estimation of the generative matrix ⇥⇤ which determines the overall structure of the network. This question is of particular interest for the task of link prediction.
Many of real-life networks are characterized by block structure. Loosely speaking, the block structure means that the nodes of the network are partitioned into groups called blocks, and that the distribution of the connections between nodes depends on the blocks to which the nodes belong. For example, when considering citation networks, where two articles are linked if one is cited by the other, it amounts to saying that the probability that two articles are linked only depends on their topic. Similarly, if one considers students of a school in a social network, it is a reasonable assumption to say that the probability that two students are linked only depends on their cohorts.
A very popular model that formalizes this idea is the stochastic block model (see, e.g., [27]). In this model, nodes are classified into k communities: each node i is associated with a community z ⇤(i), where z⇤ : [n] ! [k] is called the label function. This label function can either be treated as a parameter to estimate, or as a latent variable. In this last case, it is assumed that the indexes follow a multinomial distribution: 8i, z⇤(i) i.i.d⇠ Multinomial(1;↵⇤) where 8a 2 [k], ↵a is the probability that node i belongs to the community a. Given this label function, the probability that there exists an edge between nodes i and j depends only on the communities of i and j. Thus, the matrix of connection probabilities ⇥⇤ can be factorized as follows: ⇥⇤ij = Q ⇤ z⇤(i)z⇤(j), with Q
⇤ a k⇥ k symmetric matrix such that Q⇤ab is the probability that there exists an edge between a given member of the community a and a given member of the community b. The conditional stochastic block model can be written as:
9Q⇤ 2 [0, 1]k⇥ksym , 9z⇤ : [n] ! [k]
81  i < j  n, Aij | (Q⇤, z⇤) ind.⇠ Bernoulli ⇣ Q ⇤ z⇤(i)z⇤(j) ⌘ , Aii = 0.
(2)
Assuming that the network follows the stochastic block model, the problem of estimating the matrix of connection probabilities reduces to estimating the label function z⇤ and the matrix of probabilities of connections between communities Q⇤. Note that the conditional stochastic block model is at best identifiable up to a simultaneous permutation of the communities and of the rows and columns of the parameters Q⇤.
The stochastic block model has attracted considerable interest from the learning community. An important line of work has focused on the problem of estimation of the latent variables z⇤, see, for example, [38, 10, 1, 41]. The best understood framework is the binary, balanced, symmetric, assortative block model. In this simpler model, the two communities have the same size, the same probability of intra-community connection (Q⇤11 = Q ⇤ 22 = p), and nodes are assumed to be more connected with nodes of the same community (p > q = Q⇤12). Much work has been done on the precise characterisation of the conditions on p, q that allow for strong recovery of z⇤, i.e. to estimate z⇤ exactly with high probability. Closest to model (2) is perhaps the setting considered in [15]. In this work, the authors consider the related problem of community recovery in the binary block model [23],[2], and provide tight bounds on the recovery threshold for the balanced, two communities stochastic block model with missing observations. They propose a computationally
efficient algorithm for estimating z⇤ in regime where strong recovery is possible; this, however, requires prior knowledge of the parameter Q⇤.
Missing observations scheme Usually, when working with network data, not all the edges are observed. To account for this situation we introduce X 2 {0, 1}n⇥nsym the known sampling matrix where Xij = 1 if Aij is observed and Xij = 0 otherwise. We assume that X is random and independent from the adjacency matrix A and its expectation ⇥⇤. For any 1  i < j  n, its entries Xij are mutually independent and Xij
ind.⇠ Bernoulli(p) for some sampling rate p ! 0 such that p = ! (log(n)/n) when n ! 1.

heading: 2.2 Conditional maximum likelihood estimator
text: The log-likelihood of the parameters (z,Q) with respect to the adjacency matrix A and the sampling matrix X is given by LX(A; z,Q) = X
1i<jn Xij
⇣ Aij log(Qz(i)z(j)) + (1 Aij) log(1 Qz(i)z(j)) ⌘
= X
ab log(Qab)
X
(i,j)2z 1(a,b)
XijAij + X
ab log(1 Qab)
X
(i,j)2z 1(a,b)
Xij(1 Aij).
Let us denote by Zn,k the set of all label functions z : [n] ! [k]. For a given label function z 2 Zn,k, the log-likelihood is maximized by taking
Qab =
P (i,j)2z 1(a,b) XijAijP
(i,j)2z 1(a,b) Xij .
It is interesting to note that, for a fixed label function z, maximizing the likelihood or minimizing the least square criterion defined as CX (A; z,Q) = P
i<j Xij ⇣ Aij Qz(i),z(j) ⌘2 yields the same
estimator for the matrix Q. The main difference between these two methods is rooted in the label functions selected by the two criteria, see, e.g. [19].
To bound the risk of the maximum likelihood estimator, it is usual to assume that there exists sequences ⇢n and n such that 8i < j, 0 < n  ⇥⇤ij  ⇢n < 1. (3) This assumption ensures that the loss associated to the maximum likelihood estimator is Lipschitz continuous. See, for example, [7] and [50], where the authors assume that the adjacency matrix is generated by an homogeneous stochastic block model for which the matrix Q⇤/⇢n has entries bounded away from 0.
The restricted maximum likelihood estimator, b⇥, is based on the maximization of the likelihood among block constant matrices with entries in [ n, ⇢n]:
b⇥i<j = bQbz(i)bz(j), b⇥ii = 0
(bQ, bz) 2 argmax Q2[ n,⇢n]k⇥ksym ,z2Zn,k LX(A; z,Q). (4)
In (4), n and ⇢n are assumed to be known (see [19] for a discussion on how to estimate these parameters). Note that the Expectation-Maximization algorithm used in practice to obtain the variational approximation to the maximum likelihood estimator does not require the knowledge of these parameters. We also assume that k is known and that it can depend on the number of nodes n; it can be chosen using a network cross-validation method [12], a sequential goodness-of-fit testing procedure [33] or a likelihood-based model selection method [50]. The following result provides the upper bound on the estimation risk of the maximum likelihood estimator: Theorem 1 (Corollary 2 in [19]). Assume that A is drawn according to the conditional stochastic block model and ⇢n = !(n 1). Then, there exists absolute constants C,C 0 > 0 such that, with probability at least 1 9 exp
C⇢n
k 2 + n log(k) ,
k⇥⇤ b⇥k22  C 0 ⇣
⇢2n ((1 ⇢n)2^ 2n)
⌘ ⇢n k 2 + n log(k)
p . (5)
When all network entries are observed, we have p = 1. Note that this results implies that, when ⇢n = O( n), the maximum likelihood estimator is minimax optimal (see, [29, 16] for a statement of the lower bound).

heading: 3 Variational approximation to the maximum likelihood estimator

heading: 3.1 Definition of the estimator
text: The optimization of the log-likelihood function LX requires a search over the set of kn labels. As a consequence, the maximum likelihood estimator defined in (4) is computationally intractable. Celisse et al. [13] and Bickel et al. [7] are the first to study a variational approximation to this estimator. More recently, the authors of [48] used variational methods to approximate the maximum likelihood estimator in networks with missing observations. We start by formally introducing the variational approximation to the maximum likelihood estimator. We consider a stochastic block model with random labels with parameters (↵,Q). For this model, the likelihood of the observed adjacency matrix A and sampling matrix X is given by
lX(A;↵,Q) = X
z2Zn,k
0 @ Y
in ↵z(i)
1
A exp (LX(A; z,Q)) .
Note that the maximization of lX still requires to evaluate the expectation of the label function z for given parameters (↵,Q) by summing over kn possible labels. To circumvent this problem, one can use the mean-field approximation, which amounts to approximating the posterior distribution P (·|X A,↵,Q) by a product distribution. To ensure that this product distribution remains close to the posterior distribution, the objective function is penalized by the Kullback-Leibler divergence of the two distributions. More precisely, the posterior distribution P (·|X A,↵,Q) is approximated by a multinomial distribution denoted P⌧ , such that P⌧ (z) = Q 1in m(z|⌧ i), where m(·|⌧ i) is the density of the multinomial distribution with parameter ⌧ i = ⌧ i 1, ..., ⌧ i k , and ⌧ = ⌧ 1 , ..., ⌧ n . Then, the variational estimator is defined as ⇣ b↵V AR, bQ V AR , b⌧V AR ⌘ = argmax
↵2A,Q2Q,⌧2T JX(A; ⌧,↵,Q) (6)
for JX(A; ⌧,↵,Q) = log (lX(A;↵,Q)) KL (P⌧ (·)||P (·|X A,↵,Q))
where A, Q and T are the respective parameter spaces for the parameters ↵, Q and ⌧ , KL denotes the Kullback-Leibler divergence between two distributions, and X A denotes the observed entries of A. Since for any parameter (↵,Q), KL (P⌧ (·)||P (·|X A,↵,Q)) 0, we see that exp (JX(A; ⌧,↵,Q)) provides a lower bound on lX(A;↵,Q). The expectation - maximization (EM) algorithm derived in [48] can be used to iteratively compute the variational estimator. For details, see Appendix C. Since this algorithm is not guaranteed to converge to a global maximum, it should be initialized with care, by using, for example, a first clustering step. This solution is implemented in the package missSBM.
Statistical guarantees for the variational estimator obtained in [13, 7, 37] establish that maximizing max⌧2T JX(A; ⌧,↵,Q) is equivalent to maximizing lX(A;↵,Q), and that the estimator obtained by maximizing lX(A;↵,Q) converges to the true parameters (↵⇤,Q⇤). This in turn implies that (b↵V AR, bQ V AR ) also converges to (↵⇤,Q⇤). Note that these results do not provide guarantees on the recovery of the true labels z⇤ or on the matrix of connection probabilities ⇥⇤. In order to estimate ⇥⇤, we first define the label estimator bzV AR using the minimizer of the objective function (6):
8 i  n, bzV AR(i) , argmax ak b⌧V AR i a . (7)
Once we have estimated the community labels using (7), we replace the estimator bQ V AR
of the matrix of connection probabilities by the empirical mean estimator:
8a < k and b < k, bQ ML V AR ab ,
P (i,j)2(bzV AR) 1(a,b) XijAijP
(i,j)2(bzV AR) 1(a,b) Xij
and define b⇥ V AR as b⇥ V AR i 6=j = bQ ML V AR bzV AR(i),bzV AR(j), b⇥ V AR ii = 0. (8)
We will show respectively in Theorems 2 and 3 that this new estimator ⇣ bzV AR, bQ ML V AR⌘ is
minimax optimal for dense networks with missing observations as well as for sparse networks. The simulation study provided in Section 4 reveals that this estimator also has good performances in practice.

heading: 3.2 Convergence rates of variational approximation to the maximum likelihood estimator
text: In this section, we show the asymptotic equivalence of bzV AR and bz, where
(bQ, bz) 2 argmax Q2Q,z2Zn,k LX(A; z,Q) (9)
is the maximum likelihood estimator. More precisely, we show that, with large probability, there exists a permutation of {1, ..., k} such that bzV AR( (a))
ak = (bz(a))ak and
⇣ bQ ML V AR (a), (b) ⌘
a,bk =
⇣ bQa,b ⌘
a,bk . When this hold, the error of the tractable estimator
⇣ bzV AR, bQ ML V AR⌘ matches the
minimax optimal error rate. These results are established under the following assumptions:
A.1 There exists c > 0 and a compact interval CQ ⇢ (0, 1) such that A ⇢ [c, 1 c]k and Q ⇢ Ck⇥kQ ; A.2 The true parameters ↵⇤ and Q⇤ lie respectively in the interior of A and Q; A.3 The coordinates of ↵⇤Q⇤ are pairwise distinct.
Note that Assumption A.2 and A.3 are standard. Assumption A.2 requires that the true parameters lie in the interior of the parameter space, which is classical in parametric estimation. In the most simple case, the parameters ↵⇤ and Q⇤ lie respectively in the interior of sets A and Q of the form A = [c, 1 c]k, Q = [c0, 1 c0]k⇥ksym, for some c, c0 2 (0, 1/2). Assumption A.3 ensures the identifiability of stochastic block model parameters. Then, under the assumption that p = ! (n/ log(n)), strong recovery of the labels is possible. Assumption A.1 is more restrictive, as it implies that the network is dense. This assumption will be relaxed in Theorem 3, where we consider sparse stochastic block models such that Q⇤ = ⇢nQ0 for some fixed Q0 and some decreasing, sparsity inducing sequence ⇢n.
The following Theorem shows the minimax optimality of the tractable estimator b⇥ V AR
under assumptions A.1 - A.3. Theorem 2. Assume that A is generated from a stochastic block model with parameters (↵⇤,Q⇤) satisfying assumptions A.1 - A.3. Then, P bzV AR ⇠ bz ! 1 when n ! 1. Moreover, there exists a constant CQ⇤ > 0 depending on Q⇤ such that
P ⇥⇤ b⇥
V AR 2
2 
CQ⇤ k 2 + n log(k)
p
! !
n!1 1.
Let us now discuss the extension of Theorem 2 to the case of sparse networks. To avoid technicalities, we will consider the case when the network is fully observed. We will also assume that the proportions of different communities are held constant, while the probabilities of connections between communities may decreases at rate ⇢n. That is, the parameters (↵⇤,Q⇤) verify the following assumptions:
A.4 ↵⇤ = ↵0 for some fixed ↵0 such that ↵0a > 0 for any a 2 {1, ..., k}
A.5 Q⇤ = ⇢nQ0 for some fixed Q0 2 (0, 1)k⇥k such that kP
a,b=1 ↵ 0 a↵ 0 bQ 0 ab = 1
Assumption A.5 relaxes Assumption A.1 and allows us consider sparse networks. The normalization constraint P 1a,bk ↵ 0 a↵ 0 bQ 0 ab = 1 ensure the identifiability of the parameters (Q 0 , ⇢n) (see [7]). In the following, we denote by Q the set of parameters (↵,Q) verifying Assumptions A.4 and A.5.
The following theorem provides the analogous of Theorem 2 in the case of fully observed sparse networks. It is obtained by combining Propositions 2 and 3 in [19]: Theorem 3. Assume that A is fully observed, and is generated from a stochastic block model with parameters (↵⇤,Q⇤) satisfying Assumptions A.4 and A.5, such that Q0 has no identical columns and the sparsity inducing sequence ⇢n satisfies ⇢n log(n)/n. Then, P bzV AR ⇠ bz ! 1 when n ! 1. Moreover, there exists a constant CQ0 > 0 depending on Q0 such that
P ✓ ⇥⇤ b⇥ V AR 2
2  CQ0⇢n k 2 + n log(k)
◆ !
n!1 1. (10)
Theorems 2 and 3 establish that the variational estimator b⇥ V AR
is minimax optimal for both the estimation of dense networks with observations missing uniformly at random, and sparse networks. For proofs and discussion see Appendix A.

heading: 4 Numerical Results

heading: 4.1 Synthetic data
text: In this section we provide a simulation study of the performances of the maximum likelihood estimator defined in (8), and compare it to the variational estimator defined in [48] and implemented in the package missSBM, as well as to the Universal Singular Value Thresholding estimator introduced in [25] and implemented in the package softImpute. The results are reported in Figure 1. Thorough descriptions of the simulation protocols are provided in the Appendix.
Dense stochastic block model First, we evaluate the empirical performances of the variational approximation of the maximum likelihood estimator defined in (8) on dense stochastic block models. We estimate the matrix of probabilities of connections, and we compare our estimator with the estimator given by the methods missSBM and softImpute. The quality of the inference is assessed by computing the squared Frobenius distance between the estimators and the true matrix of connection probabilities ⇥⇤.
We consider three types of three-communities stochastic block model. The first model, given by (↵assort.,Qassort.), provides a simple assortative network, where individuals are more connected with people from their communities than with other individuals. On the contrary, the second model, given by (↵disassort.,Qdisassort.), is disassortative: individuals are more connected with individuals from outside of their communities. Both the assortative and disassortative models have balanced communities. The third model considered, given by (↵mix.,Qmix.), exhibits neither assortativity nor disassortativity, and the communities are unbalanced. We introduce missing data by observing each entry of the adjacency matrix independently with probability 0.5.
The variational approximation to the maximum likelihood estimator defined in (8) outperforms the softImpute method across all models and all number of nodes. Its error is equivalent to that of the oracle estimator with hindsight knowledge of the true label function z⇤ when the network is a few hundred nodes large. Interestingly, our estimator also outperforms the variational estimator implement in the package missSBM. We underline however that the primary focus of the missSBM method is to infer the parameters (↵⇤,Q⇤).
Additional experiments illustrating the strong consistency of the variational estimator canbe found in Appendix B.2.
Sparse stochastic block model Next, we investigate the behaviour of our estimator on increasingly sparse networks. We consider a three-communities assortative stochastic block model of 500 nodes with balanced communities, and 50% missing values. The probabilities of connections are given by Q⇤ = ⇢Q0, where ⇢ is a parameter controlling the sparsity, which ranges from 0.05 to 1. We compare the performance of the variational approximation to the maximum likelihood estimator to that of the methods softImpute and missSBM. We also compare these estimators to the trivial estimator with all entries equal to the average degree divided by the number of nodes. The error is measured as the squared Frobenius distance between the estimator and the matrix ⇥⇤ divided by ⇢2.
As the network sparsity increases, the clustering of the nodes becomes more difficult. The normalized error of the estimator b⇥ V AR increases up to a threshold corresponding to the normalized error of the trivial estimator with all entries equal to the empirical degree, divided by the number of nodes. Note that when considering very sparse networks, with ⇢ ⌧ log(n)/n, it is known that the trivial estimator with entries equal to the empirical mean degree is minimax optimal (see, eg, [29])). Thus, the estimator enjoys relatively low error rates in both high and low signal regime. By contrast, the normalized error of the softImpute method diverges as the network becomes increasingly sparse.
Stochastic block model with missing observations To conclude our simulation study, we evaluate the robustness of the methods against missing observations. We consider a three-communities assortative stochastic block model with balanced communities and 500 nodes. We increase the proportion of missing observations, and we compare the performance of the variational approximation to the maximum likelihood estimator to that of the methods softImpute and missSBM. The error is measured as the squared Frobenius distance between the estimator and the matrix ⇥⇤.
As the sampling rate p decreases, the clustering becomes impossible and the error rate of the estimator b⇥ V AR increases up to that of the trivial estimator obtained by averaging the observed entries of the
adjacency matrix. By contrast, the methods softImpute and missSBM lack robustness against missing observations, and their error diverges as the number of missing observations increases.

heading: 4.2 Analysis of real networks

heading: 4.2.1 Prediction of interactions within a elementary school
text: We apply our algorithm to analyze a network of interactions within a French elementary school collected by the authors of [47]. The network records durations of physical interactions occurring within a primary school between 222 children divided into 10 classes and their 10 teachers over the course of two consecutive days; this dataset was collected using a system of sensors worn by the participants. We consider that an interaction has occurred if the corresponding duration is greater than one minute. If an interaction of less than one minute is observed, we consider that this observation may be erroneous, and treat the corresponding data as missing. By doing so, we remove respectively 11 and 13% of the observations on Day 1 and Day 2.
The graphs of interactions recorded during Day 1 and Day 2 can be considered as two outcomes of the same random network model characterized by the matrix of connection probabilities ⇥⇤. In this spirit, we use the observations collected on Day 1 estimate the matrix ⇥⇤, and evaluate those estimators on the network of interactions corresponding to Day 2. We note that the network of interactions for Day 1 has rather homogeneous degrees (the maximum degree is 41 and the minimum degree is 5, while the mean degree is 20). Moreover, it exhibits a strong community structure. Therefore, we expect the networks of interactions to be well approximated by a stochastic block model.
We compare the performance in terms of link prediction of the estimator b⇥ V AR
defined in (8) to that of the method missSBM, and that of the method softImpute. In this last method, we set the penalty to 0, and we choose the rank of the estimator to be equal to the number of communities, which is estimated according to the Integrated Likelihood Criterion. We also compare these methods to the naive persistent estimator b⇥ naive given by b⇥ naive
ij = 1 if an interaction between i and j has been
recorded on Day 1, b⇥ naive = 0 if no such interaction has been recorded, and b⇥ naive
ij = d/n if the information is missing, where d is the average degree of the graph for Day 1. Table 1 present the error of the different estimators, measured as the squared Frobenius distance between the adjacency matrix of Day 2 and its predicted value, divided by the squared Frobenius norm of the adjacency matrix of Day 2 (i.e, the error of the trivial null estimator).
The variational method predicts most accurately the interactions on Day 2. It is closely followed by the estimator provided by the package missSBM. By contrast to the simulation study, the reduction in error when using the new estimator is moderate : the error of b⇥ V AR is respectively 1.4% and 12.4% smaller than that of b⇥ missSBM and b⇥ softImpute
. In addition, the precision-recall curve presented in the Appendix indicates that no estimator is better across all sensitivity levels. Interestingly, the naive estimator obtains a high error, which suggests a certain versatility in the children’s behaviour.

heading: 4.2.2 Network of co-authorship
text: Finally, we use variational approximation to predict unobserved links in a network of co-authorship between scientists working on network analysis, first analysed in [43]. We discard the smallest connected components (with less than 5 nodes), and we obtain a network of 892 nodes. By contrast to the network of interaction in an elementary school, the network of co-authorship is quite sparse, and presents heterogeneous degrees: the average number of collaborators is 5, while the maximum and minimum number of collaborators are respectively 37 and 1.
In order to obtain unbiased estimates of the error of the estimators b⇥ V AR
, softImpute, and missSBM, we introduce 50% of missing values in the dataset. We train the three estimators on the observed
entries of the adjacency matrix, and we use the unobserved entries to evaluate their imputation error. Table 2 present the mean imputation error of the different estimators over 100 random samplings, measured in term of squared Frobenius error and normalized by the squared Frobenius norm of the adjacency matrix of the remaining entries (i.e, the error of the null estimator). Here again, the
variational approximation to the maximum likelihood estimator obtains the best performance. The precision-recall curves of these methods, included in the Appendix, indicates that this new estimator is preferable across almost all sensitivity levels. We underline however that the errors in term of Frobenius norm of the three estimators are close, and relatively high. This comes as no surprise, as the high sparsity of the network causes the link prediction problem to be difficult.

heading: 5 Conclusion
text: In this work, we have introduced a new tractable estimator based on variational approximation of the maximum likelihood estimator. We show that it enjoys the same convergence rates as the maximum likelihood estimator, and that it is therefore minimax optimal. Our simulation studies reveal the advantages of our estimator over current methods. In particular, they highlight its robustness against network sparsity and missing observations. Our results pave the way for analysing variational approximations of more general structured network models such as the latent block model.
Broader Impact
This paper proposes a new method for network analysis which is theoretically optimal, and has empirical advantages compared to current algorithms. It can be useful to researchers and practitioners with applications, e.g., in biology [40], ethnology [48] and ecology [5]. We acknowledge that improvements in network analysis and link prediction can have negative downstream societal impacts such as contributing to community-based polarization of opinion, or improving online surveillance.

heading: Acknowledgments
text: We thank the anonymous referees for their helpfull comments.

Review summaries:
The main contribution of this paper is to propose and analyze a variational approximation (VA) method for the inference of the SBM. The main result is to establish that VA leads to a tractable and minimax optimal algorithm, specifically the latter is based on the EM principle with a variational approximation scheme. Numerical experiments are also performed to verify the theory. 

This paper considers estimation of the parameters of a stochastic block model (community membership and probability parameters) using a variational approach. Both the censored, dense model and the sparse model are considered. Experiments on real and synthetic data compare the estimator to another variational estimator, as well as another estimator based on Universal Singular Value Thresholding.

This paper proves that mean-field variational inference is minimax optimal for parameter estimation of the stochastic block model (SBM) with missing data. The paper establishes this result for both dense networks as well as for (more realistic) sparse networks. The paper further reports the results of experiments on both real and synthetic data that clearly corroborate the paper's main theoretical results.