{'fraction_non_alphanumeric': 0.039310018709266224, 'fraction_numerical': 0.017023376302966756, 'mean_word_length': 4.534820209376423, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Network Estimation by Mixing: Adaptivity and More

Networks analysis has been commonly used to study the interactions between units of complex systems. One problem of particular interest is learning the network's underlying connection pattern given a single and noisy instantiation. While many methods have been proposed to address this problem in recent years, they usually assume that the true model belongs to a known class, which is not verifiable in most real-world applications. Consequently, network modeling based on these methods either suffers from model misspecification or relies on additional model selection procedures that are not well understood in theory and can potentially be unstable in practice. To address this difficulty, we propose a mixing strategy that leverages available arbitrary models to improve their individual performances. The proposed method is computationally efficient and almost tuning-free; thus, it can be used as an off-the-shelf method for network modeling. We show that the proposed method performs equally well as the oracle estimate when the true model is included as individual candidates. More importantly, the method remains robust and outperforms all current estimates even when the models are misspecified. Extensive simulation examples are used to verify the advantage of the proposed mixing method. Evaluation of link prediction performance on 385 real-world networks from six domains also demonstrates the universal competitiveness of the mixing method across multiple domains.


Introduction
Networks are widely used to represent complex interactions between entities, and network analysis has been an intensively studied topic in statistics, computer science, physics, social science, and many areas in recent years. Analyzing the structures of network data can render salient insights about relation formulation or interaction mechanisms between individuals. Network analysis has been used in various applications in science, engineering, and social studies, revealing many new aspects in these studies [Barabási and Albert, 1999, Steglich et al., 2006, Newman, 2010, Fortunato, 2010, Lü and Zhou, 2011.
Due to the noisy nature of network data, statistical modeling has been a common approach network analysis. Multiple probabilistic frameworks have been proposed for network data (e.g., Barabási and Albert [1999], Bollobas et al. [2007], Crane [2018], Ghafouri and Khasteh [2020]). These frameworks are understood as approximations of real-world network formulation from different aspects with pros and cons in various situations. Our discussion in this paper is embedded in arguably one of the most popular frameworks for statistical network modeling -the so-called "inhomogeneous Erdős-Rényi model" [Bollobas et al., 2007, Goldenberg et al., 2010, Newman, 2010. This framework includes as special cases several important models in the network literature, such as the stochastic block model (SBM) [Holland et al., 1983] and its variants [Airoldi et al., 2008, Karrer and Newman, 2011, Jin et al., 2017, Sengupta and Chen, 2018, Noroozi et al., 2019 for community structures, the random dot product model [Young and Scheinerman, 2007] and latent space model [Hoff et al., 2002, Hoff, 2008 for latent vector spaces, and the graphon model node-exchangeable graphs [Aldous, 1981, Lovász and Szegedy, 2006, Diaconis and Janson, 2008.
While effective modeling strategies have been proposed and studied in many settings, it is observed that real-world networks may exhibit a much wider range of patterns [Ugander et al., 2013, Ghasemian et al., 2020, Miao and Li, 2021 and it is generally difficult to know which model is the proper one, or whether there is a proper choice. Therefore, an indispensable task in network modeling is to select a proper model or algorithm to use. Depending on the scope of generality, model selection procedures have been proposed in different categories. The model-based methods are designed under a family of models, such as testing procedures [Bickel and Sarkar, 2016, Gao and Lafferty, 2017, Lei, 2016, Mukherjee and Sen, 2017, Banerjee and Ma, 2017, Jin et al., 2019, Zhang and Amini, 2020 and criterion-based methods [Saldana et al., 2017, Wang and Bickel, 2017, Le and Levina, 2015, Yan et al., 2017. The more general model selection approach is cross-validation [Chen and Lei, 2018, Li et al., 2020b, Chang et al., 2020, which can be used to compare models across different families. Cross-validation intuitively seeks the most proper model for the data even if none of the models under consideration are true. Despite the many good properties derived from these methods, model selection can be unstable [Breiman, 1996a] in the sense that a small perturbation of the data can result in selecting a different model. Moreover, even if one does know the correct class of candidate models, it is not always possible to effectively fit such models due to various constraints. For example, when the network is sparse, it is known that one cannot accurately estimate an SBM with a large number of communities [Rohe et al., 2011, Lei and Rinaldo, 2014, Li et al., 2020a or a general graphon structure [Chatterjee, 2015, Gao et al., 2015.
To address the aforementioned shortcomings, this paper aims for an "off-the-shelf" and principled network modeling strategy. From a practical perspective, it is computationally efficient for largescale networks and does not require careful tuning or model specification. From the theoretical perspective, our strategy can be general enough to incorporate the well-studied models in the literature with a theoretical guarantee on performance. These properties are achieved by a simple but powerful idea: instead of selecting a single model, we take a constrained linear combination of all candidate models in a data-driven way, such that the aggregated estimate is provably more stable and adaptive to different underlying models.
Model aggregation is not unique to network data. The idea has been studied in statistical learning on metric space data, such as density estimation and regression problems, often known as the "mixing" or "stacking" procedure [Stone, 1974, Wolpert, 1992, Breiman, 1996b, LeBlanc and Tibshirani, 1996, Catoni, 1997, Haussler et al., 1998, Yang, 2000, 2001, Juditsky and Nemirovski, 2000, Tsybakov, 2003, Bunea et al., 2007. However, the mixing strategy in network settings has not been well studied, except for the recent empirical work of Ghasemian et al. [2020]. Designing a network version of mixing estimation as an "off-the-shelf" method is non-trivial. The previous mixing strategies either fail to produce good theoretical guarantees in network settings [Tsybakov, 2003] or are computationally prohibitive [Yang, 2000, 2001, Tsybakov, 2003, Bunea et al., 2007, van der Laan et al., 2007 and require additional tuning [Bunea et al., 2007]. The unique nature of network data (e.g., no i.i.d sample, discreteness, sparsity) also require special methodological designs and a different approach to theoretical analysis.
Our mixing methods rely on splitting the network node-pairs into a training and a validation data set, as in Li et al. [2020b], although other methods may be considered [Spielman and Srivastava, 2011, Rohe and Qin, 2013, Le, 2021. Multiple individual models are fit using the training data and then aggregated according to the validation data. A simple aggregation approach with exponential weights is introduced to achieve the basic model adaptivity, an ideal property one could hope for the cross-validation procedure. Then, based on special geometrical properties of network data, we design a sign-constrained least square procedure to find the aggregation, which significantly improves the strategy. Our method is much simpler to implement and possesses strong theoretical guarantees. Our numerical and theoretical analyses show that the proposed mixing methods nearly achieve the oracle estimators' accuracy regardless of whether the true model is contained in the candidate model set. Furthermore, an evaluation of 385 real-world networks shows that the proposed methods are more accurate and much more computationally efficient than the stacking method of Ghasemian et al. [2020], which still lacks theoretical support.
The rest of the paper is organized as follows. Section 2 introduces the proposed methods and their theoretical properties. We start from the basic setup of the mixing strategy. Exponential mixing is first introduced as a warm-up, which is shown to achieve single-model adaptivity. Then we introduce the more general non-negative linear mixing method and show that the non-negative property is crucial in network settings. Section 3 includes extensive simulation studies of the (almost) tuningfree properties and compares the mixing estimators with several benchmark methods under both graphon (nonparametric) models and parametric models. In Section 4, we consider link prediction tasks on 385 real-world networks from 6 domains, introduced in Ghasemian et al. [2020]. The evaluation demonstrates the advantage of the proposed method as an "off-the-shelf" method for link prediction. Section 5 concludes the paper.
2 The network mixing method and its properties Notations. Given a matrix M , we use M and M F to denote its spectral norm and Frobenius norm, respectively, and let M ∞ = max ij |M ij |. In addition, for any two matrices P 1 , P 2 of the same dimensions, we define P 1 , P 2 = trace(P T 1 P 2 ). Given a matrix M ∈ R n×n and an index set Ω ⊂ [n] × [n], denote by M (Ω) the matrix obtained from M by setting all entries of M with indices in the complement Ω c of Ω to zero. For two sequences {a n } and {b n }, we write a n b n if a n = O(b n ), in which case we can also write it as b n a n . If a n b n and a n b n , we write a n ≈ b n . We say that an event E occurs with high probability if P(E) ≥ 1 − n −δ for some constant δ > 0. We use C > 0 to denote an absolute constant whose value may change from line to line.
We focus on unweighted undirected networks in this paper. Let n be the network size. A network can be represented by a binary adjacency matrix A ∈ {0, 1} n×n , where A ij = 1 if and only if node i and node j are connected. Since the network is undirected, we have A = A T . Furthermore, assume A ii = 0 for all i (we do not consider self-loops). Our discussion of networks will be embedded in the inhomogeneous Erdős-Rényi network model [Bollobas et al., 2007], which requires that the upper diagonal entries of A are independent Bernoulli random variables. Denote P = EA and assume that there exists a set M of m available methods for estimating P from A. These may be parametric estimation procedures, such as those based on the stochastic block model Lei, 2018, Li et al., 2020b] or the latent space model , or nonparametric algorithms [Chatterjee, 2015, Airoldi et al., 2013, Chan and Airoldi, 2014, Gao et al., 2015. Our goal is to leverage these methods to derive a new estimatorP that can achieve a better approximation of P in terms of the mean squared error. A natural baseline target is, not surprisingly, to ensure thatP is as good as the unknown optimal estimator in M. Such a property is often referred to as adaptivity [Catoni, 1997, Yang, 2000, 2001. Ultimately, our goal will be more ambitious than adaptivity: we would like to achieve optimality within a broader class of estimators than the m estimators of M. We aim for an off-the-shelf option for modeling networks. That is, in addition to strong theoretical properties, empirical applicability is equally important. In a linear regression setting, procedures with provable adaptivity may not be computationally feasible or may require tuning parameters [Catoni, 1997, Yang, 2000, 2001, Tsybakov, 2003, Bunea et al., 2007, van der Laan et al., 2007, but we will introduce a computationally efficient and tuning-free method that can be scalable even to large networks.
Motivated by the adaptive mixing method in the classical regression setting, we resort to a datasplitting strategy to design our method. The strategy is based on the dyad-splitting of the network cross-validation procedure of Li et al. [2020b]. Specifically, given a fixed p ∈ (0, 1), we randomly choose a subset of node pairs Ω ⊂ [n] × [n] by independently selecting each node pair with probability p; in our setting, p is usually a constant between 0 and 0.5. Then A(Ω c ) can be viewed as an adjacency matrix generated from the inhomogeneous Erdős-Renyi model with probability matrix (1 − p)P . Therefore, any model fitting methods for inhomogeneous Erdős-Renyi networks can be used with the same theoretical guarantees for consistency as their original version as discussed in Gao and Ma [2020] and Li et al. [2020c]. We apply the methods from M to A(Ω c ) to get estimators of (1 − p)P , and multiply them by (1 − p) −1 to obtain estimatorsP (1) , ...,P (m) of P . We assume that the entries of these estimators always fall between 0 and 1. If this is not the case, we can always threshold the entries without increasing the error. The mixing strategy is essentially an aggregated version of the m estimators for P .
We use the hold-out entries A(Ω) to determine a reasonable way to combine the available estimators. A potential approach for this purpose is the cross-validation method of Li et al. [2020b], which calculates A(Ω) −P (r) (Ω) 2 F , 1 ≤ r ≤ m, and picks the one with the smallest error: (In practice, one often repeats this procedure multiple times and uses the average validation error as the criterion. Here we focus on the single split for conceptual simplicity.) This is one of the most common model selection strategies in statistical modeling [Hastie et al., 2009]. It has been verified, for example, in the context of multivariate outcome prediction and density estimation [van der Laan et al., 2006, Feng and Yu, 2019, Lei, 2020. But so far it is unclear whether such an approach is valid for network cross-validation. For example, although the estimator in (1) works well for moderately dense networks when M contains the true model, its performance deteriorates quickly in a sparse regime or under model misspecification (as we show in Section 3). We now introduce a soft-selection version of the cross-validation method to address this shortcoming.

Network mixing with exponential weights and estimation adaptivity
We view the validation error A(Ω) −P (r) (Ω) 2 F as a goodness-of-fit metric for the rth model. To match the performance of the optimal model in M, we focus on the models with small valida-tion errors. In particular, we propose the following simple rule to combine the models based on exponential weights: Compared with (1), this estimator is a soft-selection version of the cross-validation procedure. Despite its simplicity, we now show that it achieves model adaptivity, matching the performance of the unknown optimal estimator produced by M.
Since Ω is independent of the data, all theoretical analysis is conditioned on Ω. The data-splitting proportion p is a fixed constant, independent of the network size n.
Theorem 1 (Mixing estimator with exponential weights). Let A be the adjacency matrix of a random network drawn from the inhomogeneous Erdős-Rényi model with probability matrix P = EA. Denote byP (r) , 1 ≤ r ≤ m, the estimators of P obtained by applying the methods from M to A(Ω), and letP be the convex combination of these estimates defined by (2). Assume m = o(n 2 ). Then with high probability, where C > 0 is a constant and Remark 1. Although Theorem 1 describes the error restricted to Ω, extending the error bound to the entire matrix P is trivial.
Replacing Ω by Ω c = [n] 2 \ Ω in the estimating procedure and applying Theorem 1 would generate another estimatorP that admits the same type of error bound. Then combiningP (Ω) andP (Ω c ) as a full estimator would produce the same error bound for the entire matrix P . For simplicity, we will state our results only in terms of P (Ω) in all theoretical discussions to follow.
Theorem 1 states that the estimator given in (2) is nearly as good as the best estimate produced by a single method from M. To better understand the additional error ε, assume that the network is generated from a stochastic block model, under which nodes are partitioned into k groups according to a label vector c ∈ [n] k and edges are formed independently between nodes with probabilities P ij = B c i c j for some fixed matrix B ∈ R k×k . According to Gao et al. [2015], while the square of the additional error is It is easy to see that ε 2 is smaller than the rate-optimal error k 2 + n log k. Therefore, our estimator can always match the optimal estimator from M.
As a connection to similar properties of statistical estimation problems, Catoni [1997] and Yang [2000Yang [ , 2001 introduce mixing methods for regression and density estimation that can achieve estimation adaptivity. Unfortunately, their strategies are computationally expensive even in regression settings, let alone in network problems, the sample size of which scales in O(n 2 ). In contrast, our estimator (2) achieves the same type of adaptivity with negligible computational cost in addition to the m model-fitting procedures.

Beyond adaptivity: Non-negative linear network mixing
The exponential mixing estimator (2) performs well when at least one of the m individual estimators is accurate. In practice, however, all the available estimators may be inaccurate, for example, when the network is very sparse. Given the m estimators at hand, it is natural to be more ambitious: to seek a better candidate than any single model estimate.
Consider the class of linear combinations of the m estimatorŝ The most natural option is the linear mixing estimator, which uses the weights provided by solving the ordinary least squares (OLS) problem: Although linear mixing performs very well in many settings, it is a generic method and does not leverage many features of network estimation. For example, all the entries of P are non-negative, and therefore it is expected that P (r) (Ω), P (Ω) > 0 for any reasonable estimateP (r) (Ω). Similarly, P (r) (Ω),P (s) (Ω) > 0 for all 1 ≤ r, s ≤ m (unless their supports are disjoint), and the angles between reasonably good estimators tend to be small. Moreover, in sparse networks the noises (entries of A − P ) are heavy-tailed random variables, so the available estimators of P can be very noisy. These observations motivate us to further improve the linear mixing method by introducing a natural non-negative constraint to the OLS weights. In particular, we consider the non-negative linear (NNL) mixing estimator based on the following weights: where π 0 denotes the constraint that all the entries of the weight vector π are non-negative. This is a simple convex optimization problem and can be solved efficiently by either a projected quasi-Newton algorithm or a sequential coordinate descent algorithm [Kim et al., 2006, Chen andPlemmons, 2010]. The NNL mixing estimator is defined to bê We can see thatP (nnl) (Ω) is the projection of A(Ω) on the convex cone formed by the conical combination of m individual estimatorsP (r) (Ω), 1 ≤ r ≤ m. The non-negative sign constraint directly imposes a regularization effect that significantly helps the method handle the potentially large number of individual estimators. This is crucial for a tuning-free procedure and matches our aim for an "off-the-shelf" method.
The non-negative constraint has proved effective in high-dimensional linear regression problems [Slawski and Hein, 2011, Slawski et al., 2013, Meinshausen et al., 2013, with properties similar to the LASSO estimator. However, in the current context, NNL estimation is not motivated by the curse of dimensionality. In network mixing problems, the sample size for (5) is of order n 2 , and m is usually much smaller than n 2 . Therefore, though m can be large if one wants to ensure the expressiveness of the mixing estimator, (5) is usually not an ultra-high-dimensional problem. Instead, the strong correlation between theP (r) s and the low concentration of the adjacency matrix (due to the network sparsity) complicate the matter and result in the deterioration of the OLS estimator, even in a relatively low-dimensional setting.
To see why the non-negative constraint can help, let us assume that the true signal P (Ω) is within the convex cone of the estimators. Figure 1a illustrates an ideal situation. When the two individual estimators form a large angle, and the perturbation range of A around P is not too large, the observed adjacency matrix is likely to stay within or close to the convex cone, so the OLS and NNL estimators perform similarly. When the perturbation range is large, as in Figure 1b (which happens in sparse networks [Le et al., 2017]), A may be far from the cone (and P ), and the projection to the cone can result in a better estimator. On the other hand, when the two estimators align well, the convex cone is smaller (Figure 1c), and A is unlikely to belong to the cone. In this case, the NNL estimator can also outperform the OLS estimator. Overall, these observations show that sparse networks and strong alignment between estimators tend to make the OLS estimator vulnerable to increasing m, even when m is much smaller than n 2 . We now present the theoretical support for the NNL estimator and its claimed advantage over linear mixing. Let be the linear subspace and the cone generated by {P (r) (Ω), 1 ≤ r ≤ m}. Denote by Π L and Π C the projections onto L and C, respectively.
Theorem 2 (Non-negative mixing estimator, positive inner products). Let A be the adjacency matrix of a random network drawn from the inhomogeneous Erdős-Rényi model with probability matrix P = EA. Denote byP (r) , 1 ≤ r ≤ m, the estimates of P obtained by applying the methods from M to A(Ω c ), and letP (nnl) be the NNL estimator defined by (7). Assume that Then there exists a constant C > 0 such that with high probability, where C > 0 is a constant and Theorem 2 states that the proposed estimatorP (nnl) (Ω) is nearly as accurate as the non-negative linear oracle Π C P (Ω). Given a fixed δ, the extra error term ε can be at most of the order log(n + m), although it will grow much more slowly if the entries of P andP (r) are of similar orders and the network is relatively sparse.
As a consequence of Theorem 2, the next corollary shows that the NNL estimator can be strictly more accurate than the OLS estimator Π L A(Ω), especially when m is large (a typical scenario in model aggregation).
Corollary 1 (Comparison of NNL and OLS estimators). Let A be the adjacency matrix of a random network drawn from the inhomogeneous Erdős-Rényi model with probability matrix P = EA. Denote byP (r) , 1 ≤ r ≤ m, the estimators of P obtained by applying the methods from M to A(Ω c ), and letP (ols) be the OLS estimator defined by (5). With high probability, Furthermore, consider the setting in Theorem 2 and denote Assume m ≥ C P (Ω) ∞ ρ −1 log n for some sufficiently large constant C > 0. Then with high probability, where Note that the term Π C Π L P (Ω) − Π L P (Ω) F is the distance from the linear oracle Π L P (Ω) to the cone C, measuring the degree of violation of the non-negative cone assumption. Corollary 1 demonstrates the effects of the main factors on performance. It can be seen that the advantage of NNL mixing over OLS mixing increases with m and also with the correlation δ, as we have discussed above. Next, we illustrate the effect of network density.
Example 1. Assume a sparse network and that all the entries of P are of the same order d/n = o(1), where d is the average degree. Also, assume that the m individual estimators are reasonable, so all theP (r) ij s are also of the same order as the P ij s. When the conical assumption holds and m n log n/d, taking the relative error to cancel the scaling effect of P , we have For sparser networks with small d, the advantage of the NNL estimator is more dramatic. This theoretical prediction is further supported by our empirical evidence in Section 3.
Theorem 2 requiresP (r) (Ω), 1 ≤ r ≤ m, to have positive inner products. This assumption is reasonable because they are all estimators of P (Ω) and must have non-negative entries. More importantly, we can directly calculate δ from data.
Next, we consider an even weaker assumption, which is true in all reasonable settings we can think of. We replace the assumption (8) in Theorem 2 with the following assumption onP (r) (Ω), 1 ≤ r ≤ m, known as the self-regularizing property [Slawski and Hein, 2011]. Denote by Σ ∈ R m×m the matrix with entries Σ rs = P (r) (Ω),P (s) (Ω) and Σ ∞ = max 1≤r,s≤m |Σ rs |. We say that Σ satisfies the self-regularizing property with constant κ, Notice that, like (8), this condition can be numerically verified from the data by solving a convex problem. More importantly, in our current setting of network mixing, this condition will almost always hold due to the following property, taken directly from the discussion of Slawski and Hein [2011].
then Σ is self-regularizing with constant κ/T .
To see why (13) is weaker than (8), notice that without loss of generality we can assume that all theP r (Ω)s have the same norm, because rescaling does not change our linear fitting. So when (8) holds, Σ also satisfies (13), with κ = ρ. Moreover, even if we are in the extreme situation where P (r) (Ω),P (s) (Ω) = 0 for some r, s, Proposition 1 indicates that (13) can still hold as long as we separate such pairs in different groups. Even in the worst case, the m-way partition guarantees the self-regularizing property, with κ = 1/m.
Theorem 3. Consider the setting of Theorem 2 with condition (8) replaced by the self-regularizing property (13). Then there exists a constant C > 0 such that with high probability, and π * is the non-negative linear coefficient of the oracle estimator such that Compared with Theorem 2, the price we pay for the weaker assumption in Theorem 3 is greater additive error. As a simple demonstration, consider the setting of Example 1 and assume that all thê P r (Ω)s have the same norm. In this case, κ = δ. At best, the of Theorem 3 has ≈ d 3/2 n 1/2 log(n + m) + d 5/4 n 1/4 log 1/2 (n + m).

Candidate set and other practical considerations
The mixing strategy involves determining a set of candidate estimators and the hold-out proportion, p. For the OLS and NNL mixing methods, the scales of the individual estimators do not matter, and they do not need to be accurate individually. This property is essential for a mixing method to be applicable for general link prediction problems (see Section 4). But exponential mixing does not make sense if theP (r) s are in the wrong scale or if they are all inaccurate.
In determining the candidate set M, there is a necessary tradeoff between computational efficiency and expressiveness. We recommend using spectral clustering together with SBM fitting for 1 ≤ k ≤ K max , spherical spectral clustering [Rohe et al., 2011, Lei andRinaldo, 2014] and DCBM fitting for 1 ≤ k ≤ K max , and the universal singular value thresholding (USVT) of Chatterjee [2015], with the improvement mentioned in , by taking the first n 1/3 singular components. This set of candidate estimators can be calculated efficiently for large networks, with the main computational burden being the one-time calculation of SVD. The block models, despite their simplicity, have been shown to have great approximation power if they are used properly [Airoldi et al., 2013], and USVT also comes with good expressiveness. When networks are sparse, the block components likely receive higher mixing weights, helping stabilize the estimator. In contrast, when networks are sufficiently dense, USVT can be more effective. Empirically, we observe that this setup gives very effective estimation performance. In the above recommendation, K max is a reasonably large positive integer, and in Sections 3 and 4, we show that our method is very robust to K max and ρ within reasonable ranges, which is again due to the regularizing property of the non-negative constraint.
The discussion so far has focused on one random split, Ω. In practice, one can also randomly repeat the procedure multiple times and take the average as the output. This may slightly improve accuracy in our limited evaluation but will also increase the computational cost by a multiplicative factor. For this reason, we only consider single splits in this study.
The mixing method, and any aggregation methods, do have limitations. The primary performance metric of the mixing approach is estimation accuracy. Combining multiple estimators may destroy the structural interpretations of each individual estimator, such as block structures or smoothness.
However, when special structures are desirable, one can always apply the same structural extraction strategy, such as community detection, to the estimatedP . A more accurate estimate of P may lead to more accurate structure extraction.

Simulation examples
We now evaluate the mixing methods using simulation examples. We focus on the task of estimating the network connection probability matrix, P . In all the examples below, we set the network size to n = 1000, with varying density. We first evaluate our method in the graphon model settings. In particular, we use the three connection graphon connection matrices W from , for which we set P = α · W and use α to control the expected average degree of the networks (Figure 2). The mixing methods are based on the recommended setting discussed in Section 2.3, using SBM and DCBM fitting for 1 ≤ k ≤ K max , and USVT. We first show that the proposed method is robust to K max within a reasonable range. After that, we compare the mixing estimator with a few benchmarks. We use K max = 15 and p = 0.1 as the default setting in all the benchmark comparisons. Given an estimatorP , performance is measured by the relative Frobenius error, defined as

Comparison of mixing aggregation strategies
We first want to compare different aggregation strategies for the mixing method, including linear mixing (OLS-m), NNL mixing (NNL-m), exponential mixing (EXP-m), and edge cross-validation model selection (ECV-m). Since K max needs to be specified, we investigate its impact on the mixing method's performance. To demonstrate that the method is almost tuning-free, we want to show that it remains stable for a reasonable range of K max values. Figure 3 shows the estimation performance with varying K max for networks with expected average node degree 20. As can be seen, small values of K max result in underfitting and a large error. As K max increases, the error drops until overfitting kicks in. However, overall, exponential mixing, NNL mixing, and crossvalidation do not suffer from overparameterization of block approximations. Linear combination, in contrast, degrades when K max is large. This observation matches our theoretical understanding of the linear estimator as we are fitting the model with an increasing number of variables. Exponential mixing significantly outperforms cross-validation, while NNL mixing is much better than both.
Next, we evaluate the robustness of the different strategies with respect to the hold-out proportion, p. We evaluate their performance when p varies from 0.1 to 0.5 and K max is fixed at 15 (Figure 4). Overall, all the aggregation methods tend to improve as p decreases within this range. This indicates that the aggregation weight determination step is easier and requires a smaller sample size than the individual model estimation step. Of all the aggregation methods, non-negative linear combination is the most robust. Linear combination also delivers competitive estimation accuracy. Based on this observation, we recommend using K max = 15, p = 0.1 as the default configuration for network mixing estimation. Now we fix K max = 15 and p = 0.1 and evaluate the estimation performance for a range of network sparsity levels, varying the expected average degree from 5 to 45 ( Figure 5). In the sparse regime, NNL mixing outperforms the others, while linear mixing catches up as the network becomes denser. This is predicted in Corollary 1 as well. Exponential mixing eventually coincides with the crossvalidation method, which also matches our theory. Overall, NNL mixing is preferable compared to the other methods. Linear mixing is less robust to the choice of K max and is inferior to non-negative combination in the sparse setting. However, when the network is denser, it outperforms the others. In all of the experiments to follow, we use K max = 15 and p = 0.1, and we believe this configuration can be used in almost all applicable tasks.

Comparison with benchmark network estimation methods
Now we compare the mixing method with a few benchmark network estimation methods. We split this into two parts. First, we consider the three graphon models in Figure 2, and we compare the mixing method to the graphon methods, which have theoretical guarantees and reasonable computational cost. These include the USVT estimator of Chatterjee [2015], the neighborhood smoothing (NS) method of , and the sort-and-smooth (SAS) method of block approximation from Airoldi et al. [2013] and Chan and Airoldi [2014]. Second, we generate networks from two special parametric models: the SBM and the latent space model [Hoff et al., 2002]. Under these models, oracle estimations can be achieved by parametric model fitting, and they are included for comparison. The mixing, USVT, NS, and latent space model fitting  are all based on the R package randnet , while SAS is based on the R package graphon [You, 2020]. Figure 6 shows the model estimation performance of the proposed mixing strategies and the three graphon estimation benchmarks. All four variant mixing strategies uniformly outperform the benchmark methods with all three models. The advantage of the mixing method is clear when the networks are sparse. Of the three graphon estimation methods, SAS is more accurate for the first two models, while USVT and NS are better for the third one. The difference between various mixing strategies is negligible. Next, we generate networks from two parametric models. The first one is the SBM, with six communities of equal size, following the configuration of . Under the SBM, we include two oracle estimators. Oracle1 requires the true community labels and estimates the probability matrix P by averaging entries within corresponding blocks of the adjacency matrix A. Oracle2 knows the true SBM with six clusters but not the true community labels; it uses spectral clustering to find the node labels and averages entries within estimated blocks of A. In particular, notice that Oracle2 itself is automatically included as one of the individual models in the mixing.
The second model is the latent space model (LSM), with where Z i , Z j are latent vectors in R 4 , generated by N (0, I 4 ) and then centralized according to . Oracle1 in this case is the oracle version of the model that uses the true model structure and also the true latent dimension, estimated by the gradient descent method of  with the recommended initialization. Oracle2 still assumes the correct model structure, but uses the wrong dimension (3 instead of 4), representing the possibility of dimensionality mis-specification. For reference, we also include an oracle version of the NNL mixing method, which has Oracle1 as an individual component. Figure 7 shows the performance of all the methods under the two parametric models. In the SBM setting, the mixing methods are inferior only to the unbeatable Oracle1 and are even better than Or-acle2. Since the mixing procedure includes Oracle2 as an individual component, this result demonstrates the effects of ensembling multiple models. In the difficult regime, including multiple models may help stabilize the estimation and further improve the estimate of the true model when it is fitted separately. In the LSM setting, the mixing method is again inferior only to the perfect oracle parametric estimation (Oracle1) in the dense setting and, again, it is even better than it in the sparse setting. USVT is also very effective in this setting, as indicated by the theory of Chatterjee [2015], and the mixing method adapts to it. NNL-m and USVT are even better than Oracle2 estimate. The oracle NNL-m adapts to Oracle1 in the dense setting and outperforms it in the sparse setting. The comparison between the mixing method and the oracle methods under these two parametric models highlights two advantages of the mixing approach: • In the sparse regime, the estimation accuracy may be poor even if the true model is known.
The mixing strategy provides a mechanism for combining simpler models to obtain a more stable estimate.
• In the dense regime, the mixing approach may remain effective even if the mixing ensemble does not include the true model and can match the oracle otherwise.

Link prediction in real-world networks
Link prediction is a widely used procedure in network science and many scientific domains to predict missing links/non-links based on a partially observed network [Liben-Nowell and Kleinberg, 2007, Lichtenwalter et al., 2010, Lü and Zhou, 2011. A link prediction algorithm usually provides a set of scores corresponding to the missing links/non-links. A higher score means that it is more likely that a link exists between the corresponding nodes. When the data are generated from a random network model, the estimated edge probabilities, if available, provide natural prediction scores. One can predict that edges exist for the node pairs with scores higher than a given threshold. Ghasemian et al. [2020] test a large body of link prediction algorithms on 550 real-world networks from six application domains (biological, economic, social, technological, information and transportation). They observe that no method provides superior prediction accuracy across all domains. To achieve reasonable adaptivity, they propose the optimal link prediction (OLP) strategy, which aggregates a large ensemble of link prediction algorithms using random forest. They first train random forest models and then refit the link prediction algorithms using all the training data. For this data set, OLP gives nearly optimal performance across all domains. In this section, we explore the adaptivity of the mixing method for this task and compare it with OLP.
We first explain how to apply the mixing method to the link prediction task. In Ghasemian et al. [2020], three categories of link prediction algorithms are used in the ensemble: topological methods, model-based methods, embedding methods. The linear and NNL mixing methods we have introduced can include the link prediction scores from topological/embedding algorithms as individual predictions and then use either linear combination or non-negative linear combination to aggregate them. Therefore, the mixing method in this situation is similar to OLP, except that we use linear methods instead of random forest to aggregate the available estimates. Conceptually, the advantages of the mixing method are straightforward. First, it is much more computationally efficient than OLP, especially for large networks and validation sets; this is because the sample size for random forest fitting in OLP scales as O(n 2 ). Second, theoretical guarantees are available for the mixing method.
The data set of all the networks and the Python implementation of OLP are provided by Ghasemian [2020]. Ghasemian et al. [2020] use 42 topological network features for link prediction. We expand this set by adding the model estimates from the canonical mixing procedure -the SBM and DCBM estimates with K = 1, 2, ..., 15, and the USVT estimate. We use three different configurations of features to investigate the effectiveness of different categories of features: link prediction based on all 73 features (42 topological plus 31 model-based); link prediction based on the 42 topological features; and prediction based on the 31 model-based features only. In the model-based category, we also include the three graphon methods and the latent space model studied in Section 3. Notice that the model-based methods are computationally efficient, while calculating some of the 42 topological features is much more time-consuming. Therefore, the topological approaches (using all 42 features) can be computationally prohibitive for large networks.
For a stable evaluation of performance, we consider only networks with more than 200 nodes, giving in 385 networks in total. Based on their origins, they are labeled as biological (72), economic (112), informational (10), social (108), technological (55), and transportation networks (28). Given each network, we randomly sample 10% (capped at 20,000) of the node pairs to form the test set and use the complement subnetwork for training. This sampling scheme is different from the one used by Ghasemian et al. [2020], where the test set is sampled in a balanced manner to maintain similar amounts of edges and non-edges. Although that can provide better predictive performance, it may be unrealistic in practice when one has no control over the test set. We therefore take the more natural approach to sample the test set randomly. We measure the link prediction performance of a method by the predictive area under the ROC curve (AUC) based on the test data and average it over ten independent repetitions. The predictive AUC is a widely used performance metric for the link prediction task [Huang and Ling, 2005]. Figure 8 shows the results for the three categories of link prediction features. The overall link prediction accuracy varies widely across domains. Social networks are easier for link prediction, while economic, technological, and transportation networks are more difficult for the same task. When topological algorithms and mode-based algorithms are implemented, the mixing method with non-negative linear combination and OLP perform similarly in five of the six categories. OLP gives slightly worse results for the sixth category, informational networks. The mixing method with linear combination is inferior to the other two in five domains but is better for economic networks. The inferiority of linear mixing may be because 73 predictors are used, and linear combination suffers from this large number of predictors. The comparison remains similar for the methods based on topological algorithms only. Overall, for model-based algorithms, the two mixing methods are still better than the other methods. The LSM and graphon estimators might be good for one category (e.g., NS for social networks) but inferior in others. The mixing methods and OLP are reasonably good in all categories. In particular, the mixing methods are still comparable to or even better than the OLP method.  Figure 9 compares the performance of each mixing strategy for different categories of features. Both OLP and non-negative linear mixing show the same pattern: the model-based prediction is worse than the topology-based prediction, while the version using both the model and topological features has similar performance as the topology-based algorithms. Linear mixing using all features is never superior in all domains. This is likely because it becomes unstable when many mixing features are used.
As practical guidance, we also want to briefly mention the speed of the above (model-based) methods. USVT is only based on the SVD of the adjacency matrix, which is usually very sparse. It is usually the most efficient one. The mixing methods need additional modeling fitting after SVD (note that the mixing needs only one round of SVD as well), but that additional cost is relatively small. So with an efficient SVD implementation [Baglama et al., 2019, Qiu andMei, 2019], SVD and USVT can easily handle networks of size O(10 5 ) on a single laptop. In principle, SAS has a similar speed, though we observe that it is slightly slower than mixing. OLP can generate all the features in the same way as the mixing methods, but random forest fitting can be much slower, and model-fitting cost becomes the major bottleneck. In our experiments, while OLP is feasible for networks with a few thousand nodes, it is much slower than mixing. The gradient method for the latent space model and neighborhood smoothing can handle networks of moderate size but may become too slow if the size is larger than 3000.
In summary, though the OLP method is claimed to be nearly optimal for link prediction, we find that the NNL mixing strategy delivers comparable or slightly better accuracy and adaptivity across all domains. Given its additional theoretical guarantees and much higher computational efficiency, we believe that the mixing method can generally serve as an off-the-shelf link prediction algorithm.

Discussion
We have introduced a mixing strategy for network estimation. It can be used as an off-the-shelf method for network modeling with both flexibility and scalability. The method is designed according to geometrical insights into the network modeling problems, and we have shown the advantage of the design theoretically. We have also demonstrated its competitive performance empirically in link prediction problems.
The current study could be extended in several directions for future follow-up. In this paper, we focus on the accuracy of the mixing estimator for a network model. In many network analysis problems, the crucial quantities can differ from the properties of the network model P itself. For example, in regression inference of network-linked data [Le and Li, 2020], the recovery of the network projection operator is critical. Intuitively, using the mixing estimator should render a more robust inference. The theory of this type of inference warrants further study. Another direction is to extend the mixing strategy to dynamic network modeling problems [Kim et al., 2018]. Flexible dynamic network models tend to be computationally challenging to fit. Given its flexibility and computational efficiency, it would be interesting to see whether the mixing method can help alleviate this difficulty.

A Exponential weight mixing
Proof of Theorem 1. From (2), we have where ∝ denotes "proportional to". Rewrite S (r) as follows: ij .
Denote ρ = max i,j P ij . Conditioning on A(Ω c ) (thus onP (r) ), S (r) is the sum of independent random variables X (r) ij with mean and variance given by Since X (r) ij ≤ 8, by Bernstein's inequality, Since m = o(n 2 ), we have for all 1 ≤ r ≤ m with high probability for a sufficiently large constant C. In other words, S (r) concentrates well around the mean (i,j)∈Ω P ij − P (r) ij 2 , and the mixture of P (r) puts exponentially large weight on the estimator P (r) with the smallest error P (r) (Ω) − P (Ω) 2 F .

Consider the index set
By the triangle inequality, m r=1 π r P (r) (Ω) − P (Ω) From the definition of I we get r∈I π r x r ≤ x + 2C ρ log n + 2C log(nm) x .
For the second sum, consider r ∈ I. Recall that π r ∝ exp(−S (r) ) and by (14), S (r) and S (r * ) concentrate around x 2 r and x 2 , respectively. Therefore with high probability we have Since r ∈ I, and consequently, Since x r ≤ n, by choosing C > 1, we get r ∈I π r P (r) (Ω) − P (Ω) F ≤ r ∈I n(nm) −2C = o(1).
In summary, we have proved that with high probability, m r=1 π r P (r) (Ω) − P (Ω) In case min 1≤r≤m P (r) (Ω) − P (Ω) F ≤ 1, we can replace log(nm)/x in the definition of I by log 1/2 (nm) and repeat the above argument. Again, for constant C > 1, the same derivations can still go through and the resulting error would be δ = 3C log n · max i,j P i,j 1/2 + 3C log 1/2 (nm).
Rewriting 3C as a new constant C completes the proof.
The proof is complete.
Proof of Theorem 3. We will mainly follow the proof strategy of Slawski and Hein [2011]. However, there is one key step in their proof (right after B.5) that may not go through. So our proof can be seen as a corrected version of that.
Note that the constraint π 0 is equivalent to δ π * , and only the second and the last terms of the expression above depend on δ, it follows that δ = argmin δ π * Xδ 2 + 2E T Xδ .
Since the objective function on the right-hand side is zero when δ = 0, its value at the minimizer δ =δ is at most zero. Equivalently, we have Let x P be the vector obtained from x by setting all entries within S − (x) to zero and x N = x − x P . The inequality above implies where δ P 1 ≤ π * 1 becauseδ = π * −π ≤ π * forπ 0. It remains to bound δ N 1 .
The proof is complete.