Inferring Degrees from Incomplete Networks and Nonlinear Dynamics

1 Introduction

Many networks, although large in sizes, are substantially incomplete [Barab´asi et al., 2002; Szilagyi et al., 2005; Woodward et al., 2010; Chen et al., 2016] when compared with the real systems. To learn the topological characteristics of a complex network, we often need powerful accessing of the vertices and edges of the whole system. However, the associated cost is prohibitive, especially for large-scale networks, e.g., social networks with millions or even billions of nodes. Therefore, we seek a reliable and cost-effective way to estimate graph characteristics from incomplete (sampled) data [Ribeiro and Towsley, 2012; Eden et al., 2018].

The individual vertex degrees provide a cornerstone to have reliable estimations of many graph statistics, including network size [Bunge and Fitzpatrick, 1993; Stumpf et al., 2008; Katzir et al., 2011; Kurant et al., 2012], average degree [Goldreich and Ron, 2008; Dasgupta et al., 2014] and motif counting [Klusowski and Wu, 2018; Kashtan et al., 2004]. The vertex degree measures the lethality of a node in protein-protein interaction networks [Jeong et al., 2001], which is usually unknown if the network is incomplete. Consequently, estimating individual degree has its practical meaning in understanding the cellular function. Nevertheless, it is a challenging problem. Although many studies are concentrating on graph sampling techniques to infer some macroscopic properties of networks, such as the degree distribution or the associated moments of degree distribution [Frank, 1980; Stumpf and Wiuf, 2005; Achlioptas et al., 2009; Ribeiro and Towsley, 2012; Zhang et al., 2015], very few work on estimating the individual vertex degrees [Ganguly and Kolaczyk, 2017] so far. Ganguly and Kolaczyk [2017] derived a class of estimators from a risk-theoretic perspective and provided analytical analysis for them. These estimators take full advantage of the local topology information in the sampled sub-network to infer the individual degrees of the sampled vertices. However, the wide variety in the network structure can make the attempt in vain if we do not have any knowledge of specific global network characteristics, such as species abundance in ecological networks, or gene expression in biological networks. Such crucial information is available, but was not used in the risk-minimization framework.

In our framework, assuming that we are given a set of coupled nonlinear equations governs dynamics of the network, and the steady-states of all the nodes or a subset of nodes under the dynamics have been observed. We propose to estimate the individual vertex degrees from its steady-state without knowing any topology of the network, and we further improve the accuracy of the prediction if we know a tiny percentage of sampled graph regardless of sampling methods, as illustrated in Figure 1. The motivation is that the steady-state of some process (e.g., an epidemic spreading, gene expression) may be observable even though the full network is not observable. For example, in social networks, we may know who got sick, but we may do not know who interacts with whom; in biological networks, we may know the gene expression of all the genes of a patient, but we may not completely know the wiring diagram of them. It is possible to infer individual vertex degrees because the steady-state is the outcome of a complex interplay between the network's topology and the dynamics. Our framework does not rely on any prior knowledge of the sampling methods or the degree distribution. However, it offers excellent inference with as low as 10% of known edges, and moderate level of noises in the measurements. Another advantage is its ability to recover individual degrees for all vertices with only partial topology, (in)accurate measurements of the state information, though the dynamical system may be misspecified. To demonstrate the impacts of sampling scheme and the measurement noise in steady-states, and how the local topology and the dynamics contribute to the overall accuracy, we conduct a series of empirical experiments on various real networks. Experimental results indicate that our approaches can accurately infer individual vertex degrees and very robust to the sampling scheme. The inference quality can be guaranteed as long as the measurements of the steady-state are not distorted too much.

Problem Setup

Let the implicit full network be G = (V,E). It is only partially observable from a sub-network G(s) = (V (s) ,E(s) ) with V (s) ⊂ V and E(s) ⊂ E. Let δ = (δ1, δ2, . . . , δn) be the true individual vertex degrees, and n = |V |. To be convenient, bold letters hereinafter in mathematical formulas represent vectors.

We consider a general dynamics in which each vertex of G evolves according to a self-driving force and a sum of interaction forces over the neighbors

Both f(·) and g(·, ·) are general and usually nonlinear, and the adjacency matrix A ≥ 0 modulates the interactions between vertices. The coupled ordinary differential equations (ODEs) which are based on the complete topology of G drive the states of the vertices to an equilibrium x∗ . Given the complete network, this equilibrium can be approximated efficiently using a mean-field approximation derived in [Gao et al., 2016]

This mean-field approximation interprets the impacts of a neighborhood using an average impact determined by xeff, which is equivalent to a homogeneity assumption. The total interactions of the neighbors of a vertex amounts to multiplying this average interaction by δi. To compute the meanfield steady-states z∗ in (2), one needs vertex degrees δi and the effective network impact, xeff, which can be computed given the complete network. Our contribution is to show how to compute these steady-states given only the incomplete network, thereby establishing a relationship between the observed steady-states and the vertex degrees, from which we back-out the individual degrees. To fill in the details, we first describe the mean-field approximation.

2 The Effective Network Impact

Now consider a vertex i and the interaction term ∑j Aijg(xi, xj) in (1), where Aij is the influence j has on i. Similarly, i influences j with a weight Aji. We define the in-degree δin i = ∑j Aij and the out-degree δout i = ∑j Aji. Assuming Aij ≥ 0, the interaction can be written as

Here, the in-degrees δin captures the idiosyncratic part, and the average g(·, ·) captures the network effect. The mean-field approximation is to replace local averaging with global averaging, which approximates the network-impact on a vertex as nearly homogeneous. Specifically, we have

where the vector g(xi, x) has the jth component g(xi, xj). We define a linear operator

to produce a weighted average of the entries in z. Then, in our mean-field approximation,

In the first order linear approximation, we can take the LA-average inside g. The average of external interactions is approximately the interaction with its average. That is LA[g(xi, x)] ≈ g(xi, LA(x)) and

where LA(x) is a global state. Let xav ≜ LA(x). Applying LA to both sides of (7) gives

According to the extensive discussion and tests in Ref. [Gao et al., 2016], the in-degrees δin and the interaction with the external xav are roughly uncorrelated, so the LA-average of the product is roughly the product of LA-averages. Thus, we have LA[δin g(x, xav)] ≈ LA(δin )LA[g(x, xav)]. Using the first order linear approximation, we take the LA-average inside f and g

Now we have

where the resilience β = LA(δin ). For undirected graphs, δin = δout , so β = ∑i δ2 i / ∑i δi = 〈δ2 〉/〈δ〉. The steadystate of (10) is the effective network impact, xeff. Plugging it into (7) gives an uncoupled ODE for xi,

In the mean-field approximation, g(xi, xj) in (1) is replaced by an interaction with a mean-field external world g(xi, xeff) and the number of neighbors impacting xi is captured by δin i . To approximately obtain the steady-states of the system, one first solves the ODE in (10) to get xeff, and then n uncoupled ODEs at each vertex to get xi, which only depends on δin i if given xeff. The method works well because the mean-field approximations only need to hold at the steady-state. Hence, we can recover the steady-state for any vertex (for example the sampled vertices) from accurate estimates δin i (δi in the undirected case) and β.

3 Our Approaches

Our approaches need the observed steady-states x∗ of nodes to estimate the individual node degrees δ, which can be improved if we also know the partial topology of the network. We consider the observed x∗ as the ground truths and aim to search an optimal combination of degrees such that the predicted steady-states z∗ of the approximated system (2) are close to x∗ . In this sense, the general techniques underlying our work can be formulated as an optimization problem:

where δ(s) i is the degree of the sampled vertex i in the incident/induced subgraph, di is the estimated missing degree for i (see more details in Section 3.2).

First, we propose the ZeroTopo algorithm to estimate the degree of each node without knowing any topology of the network, i.e., A(s) has only zero elements. Then, we design an efficient 2-step procedure - TopoPlus and Round - to improve the accuracy of degree prediction by introducing the observed topology.

• ZeroTopo estimates δ and xeff with only the observed steady states x∗ and with ZERO known topology.
• TopoPlus estimates δ and d with both the partially observed topology and x∗ . Thus, ˆ δ = δ(s) + d, meaning that the estimated degree of a node is equal to the sum of its sampled degree and the missing degree.

Algorithm 1 ZeroTopo

Input: x∗

• 1: Let t = 0 and x(0) eff = 〈x∗ 〉
• 2: Solve ˆ δ(1) i from (2) with x(0) eff
• 3: repeat
• 4: Compute ˆ β(t) = 〈ˆ δ(t)2 〉/〈 ˆ δ(t) 〉
• 5: Solve x(t) eff from (10) with ˆ β(t)
• 6: Solve ˆ δ(t+1) i from (2) with x∗ and x(t) eff
• 7: until ˆ β(t) and x(t) eff do not change

Output: ˆ

δ = ⌊ˆ δ(t) ⌉

Algorithm 2 TopoPlus

Input: x∗ , δ(s) ,A(s)

• 1: Let t = 0 and ˆ δ(t) = δ(s)
• 2: repeat
• 3: Compute ˆ β(t) = 〈ˆ δ(t)2 〉/〈 ˆ δ(t) 〉
• 4: Solve x(t) eff from (10) with ˆ β(t)
• 5: Solve d(t) i from (14) with A(s) , x∗ and x(t) eff
• 6: Update ˆ δ(t+1) i ← δ(s) i +d(t) i
• 7: until ˆ β(t) and x(t)

Output: ˆ

eff do not change

δ = ⌊ˆ δ(t) ⌉, d = ⌊d(t) ⌉ and ˆxeff = x(t) eff

• Round further optimizes (12) for integer missing degrees ˆ d, and the final estimate becomes ˆ δ∗ = δ(s) + ˆ d.

3.1 ZeroTopo

ZeroTopo refers to the estimation of vertex degree based on observed steady-states x∗ and the dynamics f(·) and g(·, ·), without any known topology information, as shown in Algorithm (1). The average steady-state is set as the initial value of the effective state xeff = 〈x∗ 〉 when ZeroTopo warms up (line 1). Based on the mean-field approximation, we can repeatedly estimate the resilience index β from estimated degrees and update the effective state such that the entire system is self-consistent. The algorithm is simple but very powerful to recover all individual degrees simultaneously.

3.2 TopoPlus

In the ZeroTopo algorithm, we assume that we have no prior information about the topology of the network. In this section, we aim to enhance the accuracy of the degree prediction when we observe a subgraph G(s) . Note that the TopoPlus algorithm only requires the topology of the observed subgraph regardless of how one samples it. We incorporate the adjacent matrix A(s) into (2) to identify the missing degree di in the discarded local topology. Hence, we get an enhanced meanfield approximation

with the steady-states x∗ of the observed vertices, and the effective network impact xeff. Once xeff is found, it is easy to

Algorithm 3 Round

Input: x∗ , δ(s) , d, ˆxeff Output: ˆ δ∗ = δ(s) + d 1: repeat 2: Let z∗ i (di) be a solution to ˙zi = f(zi) + ∑j A(s) ij g(zi, x∗j ) +dig(zi, ˆxeff) 3: Let εi(di) = Err{z ∗ i (di), x∗i } 4: Compute Q+i = εi(di +1) - εi(di), Q-i = εi(di -1) - εi(di),∀i ▷ isolated nodes are exclusive in calculating Q5: Obtain nodes I + with Q+i ≤ 0, I - with Q-i ≤ 0 6: Let I = I + ∩ I -, I + ← I + - I and I - ← I - - I 7: for each (i, j) ∈ (I +, I -) do 8: di ← di + 1, dj ← dj - 1 9: end for 10: for each i ∈ I do 11: if |Q+i | > |Q-i | then 12: di ← di + 1 13: else if |Q+i | < |Q-i | then 14: di ← di - 1 15: end if 16: end for 17: Let ˆ δ∗ = δ(s) + d, and calculate ˆ β = 〈ˆ δ∗2 〉/〈ˆ δ∗ 〉 18: Solve ˆxeff from (10) with ˆ β 19: until d does not change

derive di, that is

The term di is a compensation to δ(s) i , i.e., ˆ δi = δ(s) i + di. The idea lies in the mean-field approximation (2) that ∑j Aijg(x∗i , x∗j ) ≈ δig(x∗i , xeff). We do not have all connections δi associated with vertex i, but we have some, δ(s) i ≤ δi. Based on the same approximation (2), we know ∑j A(s) ij g(x∗i , x∗j ) ≈ δ(s) i g(x∗i , xeff), motivating us to search the missing di = δi - δ(s) i from (13) and (14). Note that for vertices not in the sampled vertices, δ(s) i = 0, ∀i ̸∈ V (s) . As a result, the enhanced approximation (13) is equivalent to (2), then computing di in (14) is exactly solving ˆ δi from (2) in ZeroTopo. TopoPlus is formulated in Algorithm (2).

3.3 Round

Both ZeroTopo and TopoPlus contain an easily unnoticed but crucial step: Round. A rounding algorithm involves the approximation of fractional values to integer ones. In specific, Eq.(14) gives a fractional di, and in our setup we consider the networks are unweighted. Therefore, we need the rounding to get an integer approximation of the degree of each node ˆ δi = δ(s) i + di. It is a complicated and hard problem. In ZeroTopo and TopoPlus (marked with ⌊·⌉), we applied the simple plain rounding algorithm because of its ease to use. However, rounding is beyond a purely numerical problem, and we should also consider its afterward impact on the predicted steady-states z∗ . Therefore, we design a heuristic combinatorial optimization based rounding scheme (see Algorithm 3) to not only ensure the accuracy in estimating individual degrees, but also minimize the error between the observed steady-states and the simulated steady-state based on the estimated degrees. The idea is to repeatedly pick a node pair (i, j), then adjust their estimated missing degrees di and dj, keeping the total number of links unchanged: increase di by 1 and decrease dj by 1 (line 8, Algorithm 3). The selection rule is that the two steady-states z∗ i and z∗ j solved from Eq. (13) do not yet match the observed states x∗ well with the current estimates di and dj. Nevertheless, after the adjustments, the steady-states z∗ becomes closer to x∗ . Though computing seems to be very intensive, each distinct missing degree can involve in solving Eq. (13) (line 2, Algorithm 3) for only once, meanwhile, the function only includes one variable. Thus, it is embarrassingly parallelizable.

Table 1: Examples of real systems governed by nonlinear dynamics. We test our algorithms on two real networks for each dynamic. The notations n and m represent the number of vertices and edges, respectively. While 〈δ〉 and β denote the average degree and resilience index respectively.

4 Experimental Results

We evaluated the performance of our approaches on six real networks, governed by three dynamical equations (ecological, regulatory, and epidemic) [Gao et al., 2016] (Table 1). For each network, we sample a fraction of the edges and repeat 10,000 times to get statistically significant results. On average, it takes only a few of seconds to converge. We employ the fraction of vertices with a relative error of not greater than 5% as a proxy metric of the estimation accuracy, and the relative error is defined as Err[δi, ˆ δi] = | log(ˆ δi/δi)|. As reported in Table 2, we show the performance of each of our algorithms. Though ZeroTopo does not take advantage of the sampled local topology information, the dynamics itself can give a decent estimate of vertex degrees even without the sampled sub-network. The TopoPlus, either equipped with plain rounding algorithm or the designed Round algorithm, can improve the accuracy of our estimation significantly, indicating that it is vital to utilize some available topology and rounding carefully.

4.1 Sampling Technique

Prior structure-based degree estimation methods [Ganguly and Kolaczyk, 2017] rely on strong assumptions about how to sample the sub-network. Specifically, these methods are designed for the induced graph from uniform sampling. In other words, it does not work for other sampling methods, for example, random walk sampling. Our approaches do not, and they are universal, regardless of sampling methods. To demonstrate this, we compare two very different ways to sample edges: uniform sampling, and random walk sampling [Hu and Lau, 2013]. Figure 2 shows that our approaches are very robust to the sampling scheme, as is expected.

4.2 Measurement Noise andModel Misspecification

Two fundamental assumptions we made are that the observed steady-states are measured precisely, and the true dynamics are known. We now show that our approaches can tolerate high level of noise added to the observed steady-states, as well as model misspecification in the form of noise added to the parameters in the dynamical system. We add multiplicative noise εx ∼ N(0, σ) to the observed states, so we observed the states x∗ (1 + εx) ≥ 0 in simulation. Similarly, the parameters in the dynamical system are perturbed by multiplying by a factor 1 + r for r ∈ [-0.5, 0.5]. Figure 3 shows the robustness to observation noise, and Figure 4 shows the impact of errors on three parameters in the ecological system (see Table 2). The parameter B describes the incoming migration rate from neighboring ecosystems, C controls the critical abundance of the Allee effect, and H depicts the positive contribution from one vertex to another. The estimates are very sensitive to the error in the interaction parameters C and H, but is very robust to the error in B.

4.3 Local Topology and Dynamics

We, for the first time, estimate the structure characteristics of graphs using the network dynamics. Then, we compare the performance of our estimators and that of the structurebased estimators [Ganguly and Kolaczyk, 2017]: the univariate risk minimizer (URM), the multivariate risk minimizer (MRM) and the Bayes estimator as the baselines. Note that our estimator is regardless of the sampling methods, and the structure-based methods are only suitable for induced subgraph sampling. To ensure the fairness of the comparison, for each run, (1) we uniformly sample the same 10% of vertices in our approaches and the structure-based estimators; (2) we only utilize the steady-state of the observed nodes in our approaches and the induced graph of the observed nodes in the structure-based estimators. We perform 10,000 runs and show the average accuracies of the baselines and our approaches in Table 3.

We have demonstrated the promising performance of our approaches on the incident subgraphs. As we discussed in Sec. 4, they can perform well no matter how the sub-network are sampled. Therefore, they are expected to accurately estimate vertex degrees on induced subgraphs as well. Moreover, in Algorithm 1, we assume that the states of all the nodes are known. In this comparison, we only know the states of the observed nodes. To bridge this gap, we have to estimate the global β from the induced subgraph. Fortunately, Jiang et. al [2020] offered a methodology to estimate β from incomplete networks, it can be used here. Without the states of all vertices, as shown in Table 3, the performance of the ZeroTopo is still outstanding compared with structure-based methods, even with 10% noises in the measurement. It indicates that the local degrees or the neighborhood in the induced subgraph alone are not sufficient to obtain reliable estimates of individual vertex degrees. Vertex states should also be explored. There is an issue in the steady states: the measurement errors. It is ubiquitous in the real-world that the observed states are suffering from measurement noises. No measurement error (i.e., σ = 0) is an ideal case, and our approaches are significantly better than the baselines. After adding a moderate amount of noise (σ = 10% on N6, N8, ME, TY, and EM; 9% on FB), our approaches can still beat them. The TopoPlus and Round provide additional improvement because of the extra topology and careful rounding.

4.4 Application: Link Prediction

Estimating degrees itself is an essential and exciting issue in network science. It is also useful for many other applications, such as link prediction, network reconstruction, identification of essential genes, influencer discovery, and many more. Here, we take link prediction as an example to show how our approach can enhance the current link prediction methods. The most intuitive approach to build links is computing the affinity between a node pair, e.g., (a) Adamic-Adar (AA) [Adamic and Adar, 2003] su,v = ∑w∈Γu∩Γv 1/ log δw, where Γu is the neighbors of vertex u, and (b) Preferential Attachment (PA) [Barab´asi and Albert, 1999] su,v = δuδv. Both metrics contain the degree terms. Therefore, we improve these methods by replacing the observed degrees δ(s) by the estimated degrees ˆ δ using our estimators. We follow a simple rule to scale up the size of the common neighbors |Γu ∩ Γv| by αu,v = min{ˆ δu/δ(s) u , ˆ δv/δ(s) v }.

The area under the receiver-operating characteristic curve (AUC) [Hanley and McNeil, 1983] is used to measure the prediction accuracy. AA requires a joint neighborhood, but PA does not. A small sampling fraction gives a very sparse subnetwork such that joint neighborhoods are absent from G(s) . As a result, the affinities of many node pairs are distorted. To evaluate the performance, we randomly sample p =1% edges for an incident subgraph and predict the remaining 99%missing links. As reported in Table 4, the link prediction of AA is almost random. The scale-up factor can not exercise an impact on AA at all since no joint neighbors are observed. For PA, once the individual degrees are accurately revealed, the link prediction can be of high quality as well. It can give a promising improvement on all six real networks given as low as 1% observed edges. It demonstrates the value of our algorithm in link prediction. Note that our approach can not only improve the current link prediction methods, but it is also essential for many other inferring issues in network science.

5 Discussion

We addressed a prevalent problem. Consider this scenario. A biologist has measured the gene expressions of all genes and collected some relationships summarized in the adjacency matrix A(s) . The biologist also knows how genes interact, the dynamical system (1). The biologist would like to identify the essential genes using the degree centrality and to make some new discoveries. Then he/she tries to use some structure-based methods to estimate the missing degrees and finds the improvement is limited. The biggest issue is that both the gene expression data and the dynamical systems are known, but he/she does not use them. In the current-state-ofthe-art, no methods are available to utilize these information. The difficulty is rooted in lacking methods to theoretically solve the n coupled nonlinear ODE equations with incomplete graph. We bridge this gap by using a mean-field approximation, where the original system can be condensed to a single dimensional system with variable xeff and topological parameter β. We provide a framework to estimate the missing degrees accurately, which enables us to measure the essential genes. Moreover, our approaches are robust to sampling methods and work even without knowing any topology.

There are several directions for future work. First, we can go for a reconstruction of the full network topology, not just vertex degrees, which needs to solve the full link-prediction problem. This is ongoing work. Second, the mean-field approximation is the basis of our algorithm. However, the estimation itself relies on the assumption that the network is homogeneous. Therefore, when the steady-states of the neighbors are unknown, we use xeff as an approximation. All nodes are treated without any discrimination, although their neighborhoods may be different from each other. This assumption may not hold in real-world systems, and the estimation accuracy will be negatively affected. Another issue is the dynamical equations. In many scenarios, it is tough to calibrate the real dynamical equations. We can employ advanced machine learning approaches to extract the governing equations from an immense amount of experimental data.