Iterative Connecting Probability Estimation for Networks
Year: 2022

Abstract:
Estimating the probabilities of connections between vertices in a random network using an observed adjacency matrix is an important task for network data analysis. Many existing estimation methods are based on certain assumptions on network structure, which limit their applicability in practice. Without making strong assumptions, we develop an iterative connecting probability estimation method based on neighborhood averaging. Starting at a random initial point or an existing estimate, our method iteratively updates the pairwise vertex distances, the sets of similar vertices, and connecting probabilities to improve the precision of the estimate. We propose a two-stage neighborhood selection procedure to achieve the trade-off between smoothness of the estimate and the ability to discover local structure. The tuning parameters can be selected by cross-validation. We establish desirable theoretical properties for our method, and further justify its superior performance by comparing with existing methods in simulation and real data analysis.

heading: 1 Introduction
text: Network analysis has been a promising research area in data science and widely applied in many fields, including medicine, social sciences. Objects such as genes, proteins, or people are presented as vertices in the network, and relationships between objects (e.g., protein interactions or friendships) are presented as connections between pairs of vertices. In this paper, we focus on estimating the connecting probability between each pair of vertices based on the observed network.
In the literature, the most fundamental assumption is the Aldous-Hoover theorem (Aldous, 1981). Assume that the network is vertex-exchangeable, then each vertex i is associated with a latent variable ξi, and the connecting probability between any two vertices depends on a graphon function f(·) of the corresponding latent variables.
To estimate the connecting probabilities, some researchers focus on estimating f(·) under strong assumptions. Chan and Airoldi (2014) propose the sorting-and-smoothing (SAS) algorithm. It first sorts the adjacency matrix using the empirical degrees and smoothes the sorted matrix via total variation minimization. It requires strict monotonicity of the degrees, which is difficult to be satisfied in practice. Since the purpose of estimating f(·) is to estimate the connecting probabilities, another direction is to estimate the latter directly, such as stochastic block models approximation (SBA) (Airoldi et al., 2013; Holland et al., 1983; Wang and Wong, 1987). Under the SBM, vertices can be divided into several blocks, and the connecting probability between two vertices depends only on the blocks they belong to. Once the vertices are clustered correctly, the connecting probabilities can be estimated by a moment estimator. In recent years, methods based on SBM approximation and its variation, network histograms, have been widely studied (Choi and Wolfe, 2014; Olhede
∗Correspondence Author, Center for Applied Statistics and School of Statistics, Renmin University of China
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
and Wolfe, 2014). Gao and Ma (2021) establish the framework of minimax error rate analysis for these methods. Since the adjacency matrix can be viewed as a noisy version of the connecting probability matrix, general matrix denoising methods can be applied too. Chatterjee (2015) provides a universal singular value thresholding method (USVT), and Xu (2018) analyzes its error rate. Zhang et al. (2017) propose an efficient neighborhood smoothing method (NS). They select those with similar connecting patterns as the neighbors for each vertex, then the connecting probabilities can be estimated by neighborhood averaging. However, the estimate is not reused to improve neighborhood selection, which may lead to information loss.
Since the main purpose of estimating the graphon function is to obtain the connecting probabilities, we focus on estimating the connecting probability matrix and develop an iterative connecting probability estimation method (ICE) in this paper. We define a new vertex distance calculated directly on the probability matrix. In contrast to one-step neighborhood smoothing method, our method is iterative. In each iteration, with an estimate, it can update the pairwise vertex distances and the sets of similar vertices. Then, with better selected similar vertices, it can improve the performance on estimation, making it easier to discover complicated structure. The initial input of the iterations can be an estimate obtained by the existing methods as we introduce above. Even if it begins with a random neighborhood averaging estimate, it can converge quickly and perform well in many cases. We also propose a two-stage strategy for neighborhood selection, which uses different sizes of the similar vertex sets during and at the end of the iterations, aiming to cope with the trade-off between the smoothness of the estimate and the ability to discover local structure. The tuning parameters can be selected by network cross-validation. According to the results on simulated and real networks, our method is comparable to state-of-the-art methods in simple cases and outperforms them in the presence of complicated structure.

heading: 2 Framework
text: Suppose G = (V,E) is the observed undirected network where V = [n] denotes the vertex set and E denotes the edge set and E ⊆ {(i, j) : i, j ∈ V, i ̸= j}. The network can be represented by an adjacency matrix A = [Aij ]n×n ∈ {0, 1}n×n, where Aij = 1 if vertices i and j are connected and 0 otherwise. We only consider undirected networks without self-loops, so A is symmetric with Aii = 0 for i ∈ V . Assume that the network is vertex-exchangeable, according to the AldousHoover theorem (Aldous, 1981), we have the following representation: for each vertex i, there exists a latent random variable ξi
i.i.d.∼ Uniform(0, 1), then for each vertex pair (i, j), we have Aij ∼ Bernoulli(Pij), where Pij = f(ξi, ξj) and f(·) is called the graphon function. Let P = [Pij ]
n×n ∈ [0, 1]n×n denote the probability matrix. Let || · ||2 and || · ||F denote the ℓ2 and Frobenius norms respectively.
Given the observed adjacency matrix A, our goal is to estimate the connecting probability matrix P. Since for each pair (i, j), Aij is the only realization of Pij , it is hard to get a sample-based estimation like maximum likelihood estimation. However, if there exists another vertex i′ satisfying f(ξi′ , ·) ≈ f(ξi, ·), then Ai′j can be roughly viewed as another realization of Pij for the reason that Pij = f(ξi, ξj) ≈ f(ξi′ , ξj) = Pi′j . Assume that there is such a group of vertices who share similar connecting patterns with vertex i; we refer to them in the rest of this paper simply as its “similar vertices”. Let S∗i denote the set of similar vertices for vertex i, that is, S ∗ i = {i′ : Pi· ≈ Pi′·}. Then we could estimate Pij by simply averaging the entries Ai′j′ over i′ ∈ S∗i and j′ ∈ S∗j , that is,
P̃ij = P̃ji = (s ∗)−2 ∑ i′∈S∗i ∑ j′∈S∗j Ai′j′ , (1)
where s∗ = |S∗i | is the number of similar vertices. Ideally, with known P, it is easy to define the set of similar vertices for vertex i as S∗i = {i′ : 0 < d∗ii′ ≤ ds∗}, where d∗ii′ = ||Pi·−Pi′·||22/n denotes the distance between vertex i and vertex i′ and ds∗ is the s∗-th smallest elements of {d∗ii′ : i′ ̸= i}.

heading: 3 Methodology

heading: 3.1 Iterative Connecting Probability Estimation
text: The neighborhood averaging method in Section 2 is only conceptually feasible. In practice, without knowledge of P, it is difficult to obtain the pairwise vertex distances and the sets of similar ver-
tices. We propose an iterative connecting probability estimation method (ICE) to solve the problem. Assume that we have an initial estimate P̂, the estimated vertex distance dii′ can be formulated as dii′ = ||P̂i· − P̂i′·||22/n. Then we can estimate S∗i with Si = {i′ : 0 < dii′ ≤ ds}, where s = |Si| is the estimated size of the similar vertex set and ds is the s-th smallest elements of {dii′ : i′ ̸= i}. Then a new estimate can be obtained by neighborhood averaging with
P̂ij = P̂ji = s −2 ∑ i′∈Si ∑ j′∈Sj Ai′j′ . (2)
Figure 1: The relationship among D, S and P̂.
Using such an estimate, we can update the pairwise vertex distances and the sets of similar vertices again, and obtain a new P̂. We repeat the process above until P̂ converges. Let D = [dij ]n×n and S = {S1, . . . , Sn} denote the pairwise vertex distances and the sets of similar vertices respectively. The iterative procedure is inspired by the interdependence among D, S and P̂, as is shown in Figure 1. We outline the details of our method in Algorithm 1.
Algorithm 1 Iterative connecting probability estimation method
Input: observed adjacency matrix A; initial connecting probability estimate P̂(0); neighborhood size s; threshold δ0 > 0. Output: connecting probability estimate P̂. 1: Let δP = +∞ and m = 0. 2: while δP > δ0 do 3: For each vertex pair i, i′ ∈ V , obtain their distance dii′ = ||P̂(m)i· − P̂ (m) i′· ||22/n.
4: For each vertex i ∈ V , obtain its set of similar vertices Si = {i′ : 0 < dii′ ≤ ds}. 5: For each vertex pair i, j ∈ V , update their connecting probability estimate P̂ (m+1)ij =
s−2 ∑ i′∈Si ∑ j′∈Sj Ai′j′ .
6: Let δP = ||P̂(m+1) − P̂(m)||F /||P̂(m)||F . 7: Let m = m+ 1. 8: end while 9: return P̂ = P̂(m).
In Algorithm 1, an initial estimate P̂(0) is required. Any P̂ via existing methods can be used as P̂(0). Alternatively, a random initial value would also work as discussed in Section 3.2, since our method is not sensitive to the initialization. For each vertex i, we get Si by randomly sampling s vertices from V \{i} without replacement, then we can obtain such a random P̂(0) by neighborhood averaging with these random selected vertices.
As to the size of the similar vertex set s, Zhang et al. (2017) set s = C(n log n)1/2 for each vertex i, where C is recommended set as 1. In Section 3.3, we develop a two-stage procedure using different levels of C in different stages, which improves the performance of our method significantly.
It should be noted that although we discuss our estimate as defined in (2) for convenience, in practice, in case that Si and Sj may overlap, we modify it as
P̂ij = P̂ji = (s 2 − |Si ∩ Sj |)−1 ∑ i′∈Si,j′∈Sj ,i′ ̸=j′ Ai′j′ .

heading: 3.2 Iterative Estimation versus One-step Estimation
text: A main contribution of our work is to estimate P iteratively rather than output a one-step estimate. As for those one-step methods (Zhang et al., 2017), they obtain the set of similar vertices based on the distances calculated from A, then obtain P̂ by neighborhood averaging. However, if we update the set of similar vertices based on P̂, the previous set of similar vertices may be changed. The
contradiction may cause inefficiency on making full use of the observed network. Indeed, if P̂ is reliable, we should update the set of similar vertices to get a better estimate until the sets remain unchanged.
We compare the performances of our iterative estimation and one-step estimation with a simulated network which consists of 1000 vertices as an example. As shown in Figure 2(a), the 1000 × 1000 connecting probability matrix is full rank with complicated structure. We use root mean squared errors (RMSE) to evaluate the performance. We present the results of the 1st and the 10th iterations with two kinds of initial inputs and varying value of tuning parameter C in Figure 2(b). It is clear that the iterations improve the precision of the estimates with both initializations.
The results also illustrate that our method is not sensitive to the initialization. Without iterations, the one-step estimate obtained by the NS method proposed by Zhang et al. (2017) outperforms random initialization significantly. However, after sufficient number of iterations, the performance of our method with different initial values becomes comparable. In practice, as the true network structure is unknown, we recommend using the random initial input.

heading: 3.3 Two-stage Neighborhood Selection
text: Setting an appropriate C is an important task. A small C is unable to make full use of the similar vertices while a large C may lead to over-smoothing. To deal with the dilemma, we propose a twostage strategy. In stage one, we use a small C = Cit during iterations to avoid over-smoothing and to discover local structure. In stage two, after the estimate converges, a large C = Cest is used to include more similar vertices to obtain a smooth estimate. The reason is that, as the iterations progress, the precision on neighborhood selection will be improved, so we should enlarge the size of the similar vertex set gradually. However, it is impracticable to select the upper bound, the lower bound and the step-size of C simultaneously. So, we simplify the problem by setting two tuning parameters. The first tuning parameter helps us better order the vertices according to their latent positions, whereas the second tuning parameter helps us obtain better estimates using the correct number of vertices.
We compare our two-stage strategy and the strategy using a fixed C (Cit = Cest) on the above simulated network. We consider the estimate of Pi·, where the index i is set as i = 350 to present the local structure. As is shown in Figure 3(a), using a small fixed C helps to obtain an estimate with small bias but introduces large variance. In contrast, according to Figure 3(b), by using a large fixed C, we obtain a smooth estimate but fail to capture the local structure. With the well-selected tuning parameters, the two-stage strategy outputs a smooth estimate of Pi· while capturing the local structure well.
We present the RMSE of the estimates with different combinations of (Cit, Cest) on the whole P over 100 repetitions in Figure 3(c). The horizontal axis stands for the value of Cest used while different curves present the results with different Cit. The curve of “Oracle” displays the RMSE when estimating P with the true similar vertices, which provides a golden standard. For all the curves, as Cest increases, RMSE decreases sharply first, reaches the bottom near Cest = 1 and increases gradually, which illustrates that using an appropriate Cest is important. Comparing the different curves, it is easily seen that with a fixed Cest, a smaller Cit will lead to smaller RMSE since it helps to discover local structure. Another idea to achieve the trade-off is to use a fixed medium
size of the similar vertex set. However, according to the orange curve, if we force Cit = Cest, no matter what the value is, this strategy is always defeated by our two-stage strategy. The performance of the combination (Cit = 0.2, Cest = 1) is the most competitive and even comparable to that of “Oracle.”
The two-stage strategy makes our proposed method more data-driven, thus estimating those connecting probability matrices with complicated structure more powerfully. To demonstrate this advantage, we display the estimates with different methods on the network above in Figure 4. The estimate of USVT is not smooth overall as the low rank assumption is not satisfied. SAS also performs poorly, as the degree is non-monotone. NS recovers the smooth areas successfully but tends to over-smooth the local structure, as Zhang et al. (2017) have claimed. Our method combines the advantages of USVT and NS. It discovers the local structure while maintaining the smoothness of the whole matrix.
For our two-stage neighborhood selection strategy, it is important to select appropriate tuning parameters. We apply the grid search method based on the network cross-validation proposed by Li et al. (2020). As Algorithm 2 displays, the edges are split into the training set Etrain and the validation set Eval. With estimate P̂ obtained from the training network G(V,Etrain), the error on the validation set can be evaluated with negative log-likelihood loss, which is defined as L = − ∑ (i,j)∈Eval [ Aij log P̂ij + (1−Aij) log ( 1− P̂ij )] . The black horizontal line in Figure 3(c) presents the RMSE with selected parameters, which shows the reliability of the crossvalidation.

heading: 4 Theoretical Properties
text: We study the theoretical properties of our iterative method. All of the results below are based on the following assumptions about the network structure.
Assumption 1 Let 0 = z0 < z1 < · · · < zK = 1, Ik = [zk−1, zk) for 1 ≤ k ≤ K − 1 and IK = [zK−1, zK ]. Assume that the graphon function f : [0, 1]2 → [0, 1] is bi-Lipschitz on each
Algorithm 2 Tuning parameters selection of ICE via edge cross-validation Input: observed adjacency matrix A; the training proportion p; the candidate set T of the tuning
parameters (Cit, Cest); the number of replications M . Output: a combination of the tuning parameters (Ĉit, Ĉest).
1: for m = 1, . . . ,M do 2: Randomly sample a subset of edges from E with probability p to obtain the training set of
the edges Etrain. Let Eval = E − Etrain denote the validation set. 3: Apply matrix completion method to the adjacency matrix corresponding to the leftover train-
ing network G(V,Etrain) to obtain Atrain. 4: Apply ICE method on Atrain to estimate P with each (Cit, Cest), and calculate the negative
log-likelihood loss Lm(Cit, Cest) of the estimate on Eval. 5: end for 6: Let L(Cit, Cest) = ∑M m=1 Lm(Cit, Cest)/M . 7: return (Ĉit, Ĉest) = argmin(Cit,Cest)∈T L(Cit, Cest).
Ik×Il for 1 ≤ k, l ≤ K. That is, both |f(x1, y)−f(x2, y)| ≤ L|x1−x2| and |f(x, y1)−f(x, y2)| ≤ L|y1 − y2| hold for all x, x1, x2 ∈ Ik and y, y1, y2 ∈ Il, where L is a global constant.
Assumption 2 The number of K grows with n as mink |Ik|/(n−1 log n)1/2 → ∞, where |Ik| denotes the length of Ik.
For any ξi ∈ [0, 1], let I(ξi) denote the interval that includes ξi, let Ni(∆n) = [ξi −∆, ξi +∆n] ∩ I(ξi) denote the neighborhood of ξi in which f(x, y) is Lipschitz in x ∈ Ni(∆n) for any fixed y. Then according to the Lemma 1 in Zhang et al. (2017), for arbitrary global constant C,C1 > 0, define ∆n = [ C + (C1 + 4) 1/2 ] (n−1 log n)1/2, then with probability 1 − 2n−C1/4, we have
mini∈V |{i′ : ξi′ ∈ Ni(∆n)}| ≥ C(n log n)1/2. That is, if n is large enough, the size of the similar vertex set of each vertex is larger than C(n log n)1/2 with high probability. For simplicity, we let C = 1 as Zhang et al. (2017) recommend.
Based on the assumptions, our goal is to obtain P̂ which minimizes ||P − P̂||F by neighborhood averaging. If we can get the true set of similar vertices S∗, then based on A, the corresponding estimate is P̃ with P̃ij defined as (1), which serves as the “Oracle” estimate for neighborhood averaging methods. We have the following results.
Theorem 1 Let C2, C3, C4, C5, C6 > 0 be arbitrary global constants and assume n is large enough so that n−2 + (C4 + 8)1/2n−1 ≤ [ C5 − (C2 + 3)1/2 − (C3 + 3)1/2 ] ( n−3 log n )1/4 and 8L2 [ 1 + (C1 + 4) 1/2 ] (n−1 log n) ≤ (C6 − 2C5)(n−3 log n)1/4, let C∗ =
max{C1/4, 2C2/3, 2C3/3, C4/2}, then with probability 1− 8n−C∗ , we have ||P̃−P||2Fn−2 ≤ C6(n−3 log n)1/4.
Theorem 1 shows the best error rate we can achieve once we succeed to find all the true similar vertices, which seems unrealistic in practice. Luckily, we can reach the same error rate even if the set of similar vertices includes a small number of non-similar vertices. We formulate this property in Theorem 2.
Theorem 2 With arbitrary global constants C7, C8 satisfying C8 > 2C5 + 4C7 > 0 for the C5 from Theorem 1, assume that Pij ∈ [a, b] for all (i, j) ∈ V × V where 0 < a < b < 1, for any solution (S, P̂) where S satisfies
max (i,j) ∑ i′∈Si ∑ j′∈Sj I((i′, j′) /∈ S∗i × S∗j ) ≤ √ C7n 5/8(log n)9/8 b− a = e(n), (3)
if n is large enough so that 16L2 [ 1 + (C1 + 4) 1/2 ] (n−1 log n) ≤ (C8−4C7−2C5)(n−3 log n)1/4,
then with probability 1− 8n−C∗ , we have ||P̂−P||2Fn−2 ≤ C8(n−3 log n)1/4.
According to Theorem 2, for each pair (i, j), since |Si × Sj | = s2, the maximum error rate allowed for neighborhood selection is e(n)/s2 ≈ e(n)/n log n = √ C7(n
−3 log n)1/8/(b − a), which is greater then C8(n−3 log n)1/4. That is, we can obtain an estimate with low error rate even with relatively high error rate on neighborhood selection. The smaller b− a is, the larger the error rate is allowed.
If we use the estimates with existing methods as the initial values, we can achieve a small error rate on neighborhood selection with a high probability, thus improving the estimation performance. For
any estimate P̂(0) which satisfies that maxi∈V ∑n
j=1
( P̂
(0) ij − Pij )2 n−1 ≤ C9E(n) with probabil-
ity 1−n−C10 , where the error rate E(n) satisfies limn→∞ E(n)/(n−3 log n)1/4 → ∞, we have the following theorem for the estimate in the next iteration.
Theorem 3 Define the lower bound of the distance between each vertex i and any vertex i′′ /∈ S∗i as
C(n) = min i∈V,i′′ /∈S∗i
|Pij − Pi′′j |,
assume that C2(n) ≥ 8L2 [ 1 + (C1 + 4) 1/2 ]2
(n−1 log n) + 20C9E(n), then with the global constants C2, C3, C4, C6 from Theorem 1, when we use P̂(0) as the initial value to obtain a new estimate P̂new, with probability 1− 8n−C∗ − n−C10 , we have
||P̂new −P||2Fn−2 ≤ C6(n−3 log n)1/4.
The results above show that if our proposed method start from a reasonably good estimate obtained by an existing method, then it will improve upon the initial estimate iteratively, making the final estimate comparable to that with all the true similar vertices. The next question is whether we can stop the iterations when the current estimate is near P. Assume that the current estimate P̂(m) is smooth enough to satisfy Assumptions 1 and 2, we have the following theorem.
Theorem 4 With arbitrary global constant C11 and C12 satisfying C12 > 4C5 + 2C11 > 0 for C5 from Theorem 1, and let n be large enough so that 32L2 [ 1 + (C1 + 4) 1/2 ]2 (n−1 log n) ≤
(C12 − 4C5 − 2C11) ( n−3 log n )1/4 , if P̂(m) satisfies ||P̂(m) − P||2F /n2 ≤ C11(n−3 log n)1/4, then with probability 1− 8n−C∗ , we have
||P̂(m) − P̂(m+1)||2Fn−2 ≤ C12(n−3 log n)1/4.
For our learning task, the adjacency matrix A is observed while the pairwise vertex distances and the sets of similar vertices can be viewed as unobserved data, and our goal is to estimate P with incomplete data, which is a typical task for an EM type algorithm and motivates our iterative procedure. However, although updating P̂ with S is the traditional M step, updating S using the current estimate P̂ is not equivalent to the E step. Therefore, we cannot directly use the existing theory.

heading: 5 Experiments

heading: 5.1 Simulated Networks
text: To evaluate the effectiveness of our proposed method, we compare its performance with several popular estimation methods using simulated networks with different features, including low-rank, degree monotonicity and local structure. Each network contains 1000 vertices. The connecting probability matrices are displayed in Figure 5, where the rows and the columns are ordered by ξ. Graphon 1 is low rank and has vertex degree monotonicity. Graphon 2 is roughly monotone with local structure. Graphon 3 is a periodic function. Other graphons are all full-rank with extreme values in local area like the diagonal and the corner.
We compare the following methods: universal singular value thresholding algorithm (USVT) (Chatterjee, 2015), stochastic blockmodel approximation algorithm (SBA) (Airoldi et al., 2013), sorting and smoothing method (SAS) (Chan and Airoldi, 2014), neighborhood smoothing method (NS)
(Zhang et al., 2017), our proposed method (ICE). For fair comparison, we use random neighborhood averaging estimate as the initial value. We also present the result of neighborhood averaging with true similar vertices (Oracle) as a golden standard.
We show the results on average of 100 repetitions in Table 1. For each method, we report the mean of RMSE with standard deviation in the bracket. For probability matrices with low rank (Graphons 1 and 3), USVT does a good job. However, if the graphon is full rank, the low-rank approximation by USVT becomes insufficient. Although SAS method performs well on Graphons 1 and 2 since the expected degrees is monotone, it does not work on other graphon functions. NS method performs well in many cases, which coincides with the numerical results of Zhang et al. (2017). Our proposed method is comparable to NS on simple networks generated by Graphons 1 and 3. Notably, on networks with complicated structure (Graphon 4-8), our method shows significant advantage for its ability to discover local structure.
We also display the tuning parameters (Cit, Cest) selected via cross-validation and the best Cest if we use true similar vertices (Oracle). The Cest we select is comparable to that with true similar vertices while the selected Cit is always smaller, which again shows the advantage of our method on the trade-off between the smoothness of the estimate and the ability to discover local structure.

heading: 5.2 Real Networks
text: We analyze a human brain projectome dataset from an experiment of Beijing Normal University in China (Yan et al., 2009)2. The dataset is available on https://NeuroData.io/, a platform that enables large-scale neurodata storing, analyzing, and modeling. A “projectome” is a large-scale mapping between regions of the brain. This data include 48 anonymous healthy students as controls and each of them has 3 scans in different statuses. A valuable task is to the study the relationship between the regions, which may reveal the way they work together.
We have evaluated the performance of our methods on several subjects. For the sake of brevity, we only display the result on estimating the projectome of the subject with ID 0027055 in the first status3. The network consists of 349 vertices and 3772 edges. We present the adjacency matrix in Figure 6(a). The rows and columns are ordered according to the function of the regions.
As is shown in Figure 6, the estimates of these methods are quite different. USVT tends to estimate the probability matrix as a SBM model with many blocks. It fails to recover the local structure in the top left. SBA and SAS both perform poorly, because the degree disruption does not coincide with the connecting behaviors. The result of NS is more similar to that of ICE. They both obtain a
2http://dx.doi.org/10.15387/fcp_indi.corr.bnu1 3http://mrneurodata.s3-website-us-east-1.amazonaws.com/BNU3/ndmg_0-0-48/graphs/DS00350/sub-
0027055_ses-1_dwi_DS00350.gpickle
smooth estimate with local structure. The difference is that, since the whole network is sparse, NS is easier to underestimate the connecting probability of moderately connected pairs. This shortcoming derives from its assumption about the smoothness, which may lead to over-smoothing.
As the true probability matrix is unknown, it is hard to examine and compare the results of different methods on real networks. Zhang et al. (2017) evaluate their method by applying it to the task of link prediction. Assume that A is the true adjacency matrix while Amiss is the observed version. Assume that each edge in the true network will be removed in the observed network with probability p, which can be formulated as Amiss = M ∗A, where Mij
i.i.d.∼ Bernoulli(1− p). Based on Amiss, the estimate P̂ gives the scores proportional to the connecting probabilities (Gao et al., 2016). With threshold t > 0, Zhang et al. (2017) define the false positive rate and the true positive rate by
FP(t) = ∑ i,j I ( P̂ij > t,Aij = 0,Mij = 0 ) / ∑ i,j I (Aij = 0,Mij = 0) ,
TP(t) = ∑ i,j I ( P̂ij > t,Aij = 1,Mij = 0 ) / ∑ i,j I (Aij = 1,Mij = 0) .
By varying t, we can obtain the ROC curve and compare the precision of different methods with AUC. The missing rate p is set as 0.1.
Here we develop a new approach for evaluation of connecting probability estimation methods in practice. Although P is unknown, we can estimate it by one of the existing methods and apply the result Phost as the baseline. Then we generate the adjacency matrix Agen with Phost and apply different methods to estimate it. We refer to the method used to obtain Phost as the “host” method. The “host” is always the best due to the home-court advantage. We focus on those methods whose performances are competitive and even comparable to the “host”, which demonstrates the capability of generalization.
To evaluate the precision on estimation, we display the results of the link prediction task in Figure 7(a). It is obvious that SBA and SAS are outperformed by others significantly. Our method is slightly better than NS and USVT due to its ability to recover the local structure with higher precision.
To evaluate the the capability of generalization, we use the approach we propose. For each “host”, we generate the adjacency matrix for 100 repetitions and present the RMSE on average in Figure 7(b). When Phost is based on USVT, it outperforms others significantly because Phost is low rank. However, if the “host” is our proposed method, USVT performs poorly as it is hard to recover local structure. Under this setting, NS is also defeated by our method due to over-smoothing. When the “host” are SBA and SAS, NS outperforms others slightly. The reason is that both SBA and SAS output smooth estimated probability matrices and satisfy the assumption about smoothness of NS. No matter what the “host” is, our method is always competitive, which manifests its ability of generalization.

heading: 6 Conclusions
text: We propose an iterative connecting probability matrix estimation method which updates the set of similar vertices and the estimate progressively. It is intuitive, easy to implement and computationally
feasible for its fast convergence. Even starting with an estimate by random neighborhood averaging performs well. The main strength of our method is the ability to discover local and complicated structure for networks. In our algorithm, deciding the size of the similar vertex set is crucial. We propose a two-stage strategy using different sizes to balance its ability to discover local structure and output a smooth estimate. In practice, the appropriate size may vary from vertex to vertex. A possible variation of our method is to set specific size of the similar vertex set for each vertex. Although it may introduce more turning parameters, we believe that this issue is worth of further investigation.
Acknowledgments and Disclosure of Funding
We thank the program chair and reviewers for their insightful comments that greatly improve the article. Funding in direct support of this work: None. Additional revenues related to this work: None.

Review summaries:
The authors propose an iterative approach to estimate the probability of forming an edge between any two vertices in a graph, which they denote as the *connecting probability*. With only a single realization of the graph, it is a difficult task to estimate the connecting probabilities. The proposed iterative connecting probability estimation (ICE) method uses neighborhood averaging, similar to Zhang et al. (2017), but with iterative refinement of the neighborhood. They provide a theoretical analysis demonstrating that the error rate of their proposed estimation procedure is asymptotically as good as that of an "oracle" estimate that has access to the true connecting probabilities to choose the neighborhood. They also provide extensive simulations and a real data experiment demonstrating good empirical performance compared to prior work. 

This paper deals with the task of estimating the connection probabilities in a random graph, assuming that it is drawn from a graphon model; i.e. $P_{ij} = f(\xi_i, \xi_j)$, where $f$ is a smooth function and the $\xi_i$ are uniform random variables. The authors do so by exploiting the same ansatz as in Zhang et al: for $i \in [n]$, find several $i'$ with connection probabilities "close" to those of $i$, and average the $P_{i'j}$ to find an estimate for $P_{ij}$.  The main difference with the existing litterature is the choice of the $i'$: the authors choose an iterative procedure, building ever improving estimators of the $P_{ij}$ and using the current estimate to find the $i'$ and build the next estimator. The paper also presents a two-step procedure to choose the right number of "close" vertices.  The authors then show several theoretical properties of their algorithm, in a model including (but not limited to) the classical stochastic block model. Their bounds depend on the number of vertices correctly classified as "close" to each vertex in the graph.  Finally, they present experiments on both synthetic and real-world data; these seem to indicate that the ICE algorithm outperforms its counterparts when the graphon function $f$ becomes more complex. 

The paper studies an interactive algorithm for connecting probability estimation. The proposed algorithm updates the pairwise vertex distance, neighborhood set, and the connecting probabilities iteratively to improve the estimation until the result converges. The authors have shown both theoretical and empirical results. Some practical techniques to improve the performance of the algorithm are also discussed in the paper.