Variational Bayesian Weighted Complex Network Reconstruction

TL;DR

This paper introduces a variational Bayesian framework for reconstructing weighted complex networks, addressing noise issues where traditional lasso methods underperform, and demonstrating improved accuracy and efficiency.

Contribution

It proposes a novel variational Bayesian approach for network reconstruction that outperforms lasso in noisy environments, with faster computation and higher accuracy.

Findings

Outperforms lasso in reconstruction accuracy

Faster inference speed than traditional methods

Effective on both synthetic and real data

Abstract

Complex network reconstruction is a hot topic in many fields. Currently, the most popular data-driven reconstruction framework is based on lasso. However, it is found that, in the presence of noise, lasso loses efficiency for weighted networks. This paper builds a new framework to cope with this problem. The key idea is to employ a series of linear regression problems to model the relationship between network nodes, and then to use an efficient variational Bayesian algorithm to infer the unknown coefficients. The numerical experiments conducted on both synthetic and real data demonstrate that the new method outperforms lasso with regard to both reconstruction accuracy and running speed.

Tables

Table 1. Table 1: The list of notations in the paper.

Notation	Description
$𝑨,𝑾$	The adjacency and weighted matrix of a network
$\widehat{𝑨},\widehat{𝑾}$	The reconstructed network
$𝒚$	Response vector
$𝑿$	Design matrix
$𝒘$	Regression coefficient vector
$𝒂$	Binary regression coefficient vector
$N,M$	The number of nodes and observed data, respectively
$ϵ,{\sigma }^2$	Noise item and its variance
$c_0,d_0,g_0,h_0$	Parameters in Gamma distribution
$e_0,f_0$	Parameters in Beta distribution
$\rho$	A parameter in Bernoulli distribution
${\lambda }_j$	A parameter in the distribution of $w_j$

Table 2. Table 2: Results on empirical networks (Experiment 4). ECT is simulated on the first two networks and communication is simulated on the rest ones.

Networks	$N$	Method	TPR	TNR	Error	Time
HB494BUS	494	Lasso	0.972	1.000	0.028	238.549
		VBR	0.987	1.000	0.002	73.710
HB1138BUS	1138	Lasso	0.948	1.000	0.027	3381.924
		VBR	0.977	1.000	0.001	1553.495
Jazz	198	Lasso	0.981	0.998	0.031	49.361
		VBR	0.969	1.000	0.003	0.913
Karate	34	Lasso	0.983	0.999	0.031	0.665
		VBR	0.997	1.000	0.000	0.015

Equations

$$\sum _{j=1}^N\frac{a_{ij}}{r_{ij}}(V_i-V_j)=I_i(i=1,2,\cdots ,N),$$

$$\left[\begin{array}{c}{y}_{t_1}^{(i)}\\ y_{t_2}^{(i)}\\ ⋮\\ y_{t_M}^{(i)}\end{array}\right]=\left[\begin{array}{cccc}{x}_{1,t_1}^{(i)}& x_{2,t_1}^{(i)}& \cdots & x_{N,t_1}^{(i)}\\ x_{1,t_2}^{(i)}& x_{2,t_2}^{(i)}& \cdots & x_{N,t_2}^{(i)}\\ ⋮& ⋮& & ⋮\\ x_{1,t_M}^{(i)}& x_{2,t_M}^{(i)}& \cdots & x_{N,t_M}^{(i)}\end{array}\right]\left[\begin{array}{c}{a}_{i1}{w}_{i1}\\ a_{i2}{w}_{i2}\\ ⋮\\ a_{iN}{w}_{iN}\end{array}\right].$$

$${\mathbf{y}}^{(i)}={\mathbf{X}}^{(i)}\mathbb{D}({\mathbf{a}}^{(i)}){\mathbf{w}}^{(i)},$$

$$\mathbf{y}=\mathbf{X}\mathbb{D}(\mathbf{a})\mathbf{w}+\mathbf{ϵ},$$

$$\mathbf{y}\mid \mathbf{a},\mathbf{w},\tau \sim \mathcal{N}\left(\mathbf{y}\mid \mathbf{X}\mathbb{D}(\mathbf{a})\mathbf{w},{\tau }^{-1}{\mathbf{I}}_M\right),$$

$$w_j\sim \mathcal{N}(w_j\mid 0,{\lambda }_j^{-1}),j=1,2,\cdots ,N,$$

$$a_j\mid \rho \sim \mathrm{Bernoulli}\left(a_j\mid \rho \right),j=1,\cdots ,N,$$

$$p\left(\mathbf{y},\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho \right)=p\left(\mathbf{y}\mid \mathbf{a},\mathbf{w},\tau \right)p(\tau )\left[\prod _{j=1}^{N}p(a_j\mid \rho )p(w_j\mid {\lambda }_j)p({\lambda }_j)\right]p(\rho ).$$

$$p\left(\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho \mid \mathbf{y}\right)=p\left(\mathbf{y},\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho \right)/p\left(\mathbf{y}\right).$$

$$p\left(\mathbf{y}\right)$$

$$=\int p\left(\mathbf{y}\mid \mathbf{a},\mathbf{w},\tau \right)p(\mathbf{a}|\rho )p(\mathbf{w}|\mathbf{\lambda })p(\mathbf{\lambda })p(\tau )p(\rho )\mathrm{d}\mathbf{a}\mathrm{d}\mathbf{w}\mathrm{d}\mathbf{\lambda }\mathrm{d}\tau \mathrm{d}\rho$$

$$q^{\ast }=\underset{q(\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho )}{\min }KL\left(q(\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho )||p(\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho \mid \mathbf{y})\right).$$

$$q(\mathbf{a},\mathbf{w},\mathbf{\lambda },\tau ,\rho )=q(\mathbf{w})q(\rho )q(\tau )\prod _{j=1}^{N}q({\lambda }_j)q(a_j).$$

$$q(\mathbf{w})$$

$$q(\tau )$$

$$q(\rho )$$

$$q(a_j)$$

$$q({\lambda }_j)$$

$$\mathbf{\Omega }$$

$$\mathbf{\Sigma }$$

$$g_j$$

$$h_j$$

$$c$$

$$d$$

$${\theta }_j$$

$$u_j$$

$$e$$

$$\mathrm{TPR}=\frac{\mathrm{TP}}{\mathrm{P}}=\frac{{\sum }_{i=1}^N{\sum }_{j=1}^N{a}_{ij}{\hat{a}}_{ij}}{{\sum }_{i=1}^N{\sum }_{j=1}^N{a}_{ij}},$$

$$\mathrm{TNR}=\frac{\mathrm{TN}}{\mathrm{N}}=\frac{{\sum }_{i=1}^N{\sum }_{j=1}^N(1-a_{ij})(1-{\hat{a}}_{ij})}{{\sum }_{i=1}^N{\sum }_{j=1}^N(1-a_{ij})},$$

$$\mathrm{Error}=\frac{\sqrt{{\sum }_{i=1}^N{\sum }_{j=1}^N({\hat{w}}_{ij}-w_{ij}{)}^2}}{\sqrt{{\sum }_{i=1}^N{\sum }_{j=1}^N{w}_{ij}^2}}.$$

$$s_{it}=\sum _{k\ne i}{a}_{ik}{w}_{ik}{s}_{kt}+{ϵ}_t,$$

$${\mathrm{CI}}_i=\frac{{\sum }_{k\ne i}(a_{ik}+a_{ki})I(l_i=l_k)}{{\sum }_{k\ne i}(a_{ik}+a_{ki})I(l_i\ne l_k)}.$$

$$\mathrm{NMI}=\frac{2I(\hat{\mathbf{l}};\mathbf{l})}{H(\hat{\mathbf{l}})+H(\mathbf{l})},$$

$$q^{\ast }(\mathbf{z})=\min KL(q(\mathbf{z})\Vert p(\mathbf{z}\mid \mathbf{y})).$$

$$KL(q(\mathbf{z})||p(\mathbf{z}\mid \mathbf{y}))=$$

$$\log p(\mathbf{y})=KL(q(\mathbf{z})||p(\mathbf{z}\mid \mathbf{y}))+\mathrm{ELBO}.$$

$$\mathrm{ELBO}=-\int \prod _{j}q(z_j)\log \frac{{\prod }_{j}q(z_j)}{p(\mathbf{z},\mathbf{y})}\mathrm{d}\mathbf{z}-\int q(z_j)\{\log q(z_j)-{\mathbb{E}}_{-q(z_j)}\log p(\mathbf{z},\mathbf{y})\}\mathrm{d}{z}_j+const.,$$

Full text

1 Introduction

The networked-systems are ubiquitous in many fields, including social-tech science [16, 38], bioinformatics [27, 25, 3, 26], epidemic dynamics [19, 37, 39] and power grid [15, 17]. As is often the case, it is unable to observe the topology of a network, while data generated by this network is available. Therefore, in interdisciplinary science, one of the most important but challenging problems is to reconstruct the complex network from observed data or time series [28].

Suppose that a complex network consists of $N$ nodes, and we are given the time series of the states for the $N$ nodes. A decade ago, the main network reconstruction technique was causality analysis (CA) [6, 8, 9, 34, 35] which infers the causal influence between two variables via a pair of linear regression models. With the rapid development of variable selection and feature learning [18, 32, 40], CA has been replaced with lasso [22] (a.k.a. compressive sensing). Wang et al. [29] made the first attempt to apply lasso to reconstructing networks with game-theoretic dynamics. Thereafter, Wang’s group applied lasso to networks with epidemic [20], geospatial [21] and electrical power [10] data. They also studied the application to related problems (e.g., the prediction of catastrophes in nonlinear dynamical systems [30]). Recently, Wu’s group has extended this framework to multilayer networks [14] and the networks with time-varying nodal parameters[31]. To some extent, lasso is a “black-box” tool for the complex network reconstruction task. Recently, an increasing number of researchers start to develop alternative techniques to lasso. For example, Ma et al. and Xiang et al. cast the problem into a statistical inference issue [13, 36].

Essentially, lasso [22] solves an $L_{1}$-norm penalized least squares problem, i.e., $\min_{\bm{w}\in R^{N}}\|\bm{y}-\bm{Xw}\|^{2}_{2}+\lambda_{\rm lasso}||\bm{w}||_{1}$. Owing to the property of the $L_{1}$-norm penalty [7], the solution of $\bm{w}$ is sparse, that is, many entries of $\bm{w}$ are zero. The parameter $\lambda_{\rm lasso}>0$ controls the sparsity, and larger $\lambda_{\rm lasso}$ makes $\bm{w}$ sparser. In the context of (unweighted) network reconstruction, $\bm{y}$ and $\bm{X}$ are observed data, and the element of regression coefficient vector $\bm{w}$ indicates whether a pair of nodes are connected or not. In other words, $\bm{w}$ corresponds to a column of the adjacent matrix $\bm{A}$ of a network. Therefore, the network topology can be recovered column-by-column via applying lasso $N$ times.

Unfortunately, the current researchers pay little attention to the weighted network reconstruction task. In the cases of weighted networks, besides determining whether each pair of nodes are connected or not, it is also required to estimate the connection strength. To meet the demands, lasso based framework has to carefully select $\lambda_{\rm lasso}$. In general, a good $\lambda_{\rm lasso}$ can be obtained by cross validation (CV) technique [1]. But it is a time-consuming strategy. Additionally, lasso in theory cannot precisely estimate the connection strength owing to lasso’s biased estimation property [22]. A toy example shown in Fig. 1 illustrates the drawback of lasso. To the best of our knowledge, there are only two references related to weighted network reconstruction, but they focus on small-scale networks [12, 41] (the number of nodes is less than 10).

By putting the network reconstruction task into the hierarchical Bayesian modeling framework, we in this paper propose an elegant solution to weighted complex network reconstruction, especially for large-scale networks. This work contains the following threefold contributions.

(i) Firstly, to overcome the shortcoming of lasso, we consider a special regression problem with multiple regression coefficients, that is, $\min_{\bm{w},\bm{a}}\|\bm{y}-\bm{X}(\bm{a}\odot\bm{w})\|^{2}$, where $\odot$ denotes the element-wise product, $\bm{w}\in R^{N}$ plays the same role as the traditional regression coefficient vector, and $\bm{a}\in\{0,1\}^{N}$ indicates variables are active or not. In the context of complex networks, the real-valued $\bm{w}$ is used to estimate the connection strength, while the binary-valued $\bm{a}$ represents whether nodes are connected or not. In later discussions, it will be found that this formulation is very suitable for handling large-scale network reconstruction problems.

(ii) Secondly, a full hierarchical Bayesian inference technique is used to estimate the unknown parameters ${\bm{w}}$ and ${\bm{a}}$. By the virtue of hierarchical Bayesian inference, our model overcomes the shortcoming of lasso, i.e., tuning hyper-parameters. To speed up the inference, a variational Bayesian technique rather than Markov chain Monte Carlo (MCMC) sampling is employed to approximate the posterior distributions of unknown variables [5, 11, 24]. To the best of our knowledge, this is the first attempt to utilize the variational Bayesian technique to reconstruct networks.

(iii) The experiments conducted on simulated networks (scale-free and small-world networks with electrical current transportation (ECT) and communication dynamics) and a stock network show that our framework significantly outperforms lasso.

The rest of this paper is organized as follows. In Section 2, we provide the definition of the network reconstruction task. In Section 3, we formulate the new framework and show how to infer the network topology with it. Then, some experiments are conducted in Section 4. Finally, the paper ends with conclusions and future works.

2 Problem Statement

To facilitate later discussions, Table 1 lists the main notations and symbols used in the paper. By following the common practice of many references [10, 12, 14, 20, 21, 29, 30, 31, 41], we consider a weighted complex network without loops. Suppose that its weighted matrix is $\bm{W}=(w_{ij})$, where $w_{ij}$ represents the connection strength from node $i$ to $j$. We set $w_{ij}=0$ if nodes $i$ and $j$ are not connected. Its adjacency matrix is $\bm{A}=(a_{ij})$, where $a_{ij}=1$ if $w_{ij}\neq 0$ and 0 otherwise. Moreover, we assume that the nodes are governed by a specific dynamics. As defined in the reference [28], the network reconstruction task is to infer the network topology according to observed data or time series from nodal dynamics.

Here, to facilitate illustration, we take the ECT in a power network consisting of resistors as a special example to introduce the network reconstruction task. The resistance of a resistor between nodes $i$ and $j$ is denoted by $r_{ij}$, where $r_{ij}=\infty$ if they are not connected. Based on the Kirchhoff’s law, there is

[TABLE]

where $V_{i}$ and $I_{i}$ denote the voltage and the electrical current of node $i$, respectively. Generally speaking, in the real world, we are able to observe the nodes’ voltage and electrical current, while the network structure is invisible. If we record the voltage and electrical current at different time points, the reconstruction task for this power network is to infer the network topology (including $a_{ij}$ and $r_{ij}$) based on the observed data.

3 Bayesian complex network reconstruction

3.1 Reconstructing network by regression

In this subsection, we show that the network reconstruction task can be accomplished by solving a series of regression problems. In what follows, we define $x_{j}^{(i)}=V_{i}-V_{j}$ and $w_{ij}=1/r_{ij}$. At the same time, we use $y^{(i)}$ to refer to $I_{i}$. Since the network simulates a dynamic system, it is reasonable to assume that data are collected at $M$ different time points. To facilitate description, let $\{t_{m}\}_{m=1}^{M}$ denote the time index set. With the above assumptions, according to Eq. (1), the variables should satisfy

[TABLE]

It is noteworthy that the item $y_{t_{m}}^{(i)}=I_{i,t_{m}}$ represents the electrical current of node $i$ and $x_{j,t_{m}}^{(i)}=V_{i,t_{m}}-V_{j,t_{m}}$ denotes the difference of voltage between nodes $i$ and $j$ at time $t_{m}$. In the matrix form, there is

[TABLE]

where $\mathbb{D}(\bm{a}^{(i)})$ is a diagonal matrix whose main diagonal is $\bm{a}^{(i)}$, and $\bm{a}^{(i)}=\left(a_{i1},\cdots,a_{iN}\right)^{\rm T}$ is the $i^{\rm th}$ row of matrix $\bm{A}$. As stated above, the nodal states are observable. In other words, $\bm{y}^{(i)}$ and $\bm{X}^{(i)}$ in Eq. (3) are observed data; $\bm{a}^{(i)}$ and $\bm{w}^{(i)}$ are unknown items. The network’s weighted matrix can recovered column-by-column via solving Eq. (3) with $i=1,\cdots,N$. In this manner, the network reconstruction can be cast into a series of regression problems.

3.2 Model formulation

Since the network reconstruction actually corresponds to solving some special regression problems, in this subsection we cast it into a new framework based on Bayesian statistics.

In the following discussions, we reformulate the estimation problem in Eq. (3) as a linear regression model via

[TABLE]

where $\bm{\epsilon}$ denotes the noise item. The response vector $\bm{y}\in R^{M}$ and the design matrix $\bm{X}\in R^{M\times N}$ are observable, while the binary coefficient vector $\bm{a}$ and the continuous coefficient vector $\bm{w}$ need to be estimated. Here $\mathbb{D}(\bm{a})$ denotes a diagonal matrix whose main diagonal is the vector $\bm{a}$. We emphasize that this framework can handle weighted networks with many kinds of dynamics [28], including communication dynamics [10], evolutionary game dynamics [29], epidemics [20], synchronization dynamics [14], to name but a few. Besides ECT, we introduce how to apply our framework to communication dynamics in Appendix B.

In general, it can be assumed that $\epsilon_{t_{m}}$ are independently and identically distributed (i.i.d.) as Gaussian, namely, $\epsilon_{t_{m}}\mid\sigma\sim\mathcal{N}(\epsilon_{t_{m}}\mid 0,\sigma^{2}),m=1,\cdots,M$, where $\sigma^{2}$ denotes the variance. Note that the popular reconstruction method lasso also falls into this category. The significant advantage of lasso over the simple least-squares method lies in that it can provide a sparse solution. However, its performance highly depends on the tuning of its parameter $\lambda_{\rm lasso}$. In contrast, Bayesian methods can provide satisfactory estimates for $\bm{w}$ and $\bm{a}$ while avoiding the tedious parameter adjustment. The core idea of Bayesian methods is to impose a prior distribution on each unknown variable (i.e., $\bm{w}$ and $\bm{a}$), then MCMC sampling or a variational technique is employed to approximate the posterior distribution according to the famous Bayes theorem. There are hyper-parameters in the prior distribution sometimes, and some proper hyper-prior distributions can be hypothesized. In what follows, we will adopt a full hierarchical Bayesian inference process to infer our interested items in Eq. (4).

To facilitate illustration, we let $\tau^{-1}\equiv\sigma^{2}$ and $\tau$ is often called precision. Thereafter, the conditional distribution of $\bm{y}$ is given by

[TABLE]

where $\bm{I}_{M}$ denotes an identity matrix of order $M$. For each regression coefficient (connection strength) $w_{j}$, we place the following Gaussian prior

[TABLE]

where $\lambda_{j}^{-1}$ is the variance. Since each entry $a_{j}$ of the coefficient vector $\bm{a}$ takes binary value $\{0,1\}$, it is natural to consider a Bernoulli prior for $a_{j}$, that is,

[TABLE]

where $\rho$ denotes the probability of $a_{j}$ taking value 1. At last, to make a full Bayesian inference, we impose conjugate priors on the parameters $\lambda_{j}$, $\tau$ and $\rho$, namely, $\tau\sim{\rm Gamma}(\tau\mid c_{0},d_{0}),\rho\sim{\rm Beta}(\rho\mid e_{0},f_{0}),\lambda_{j}\sim{\rm Gamma}(\lambda_{j}\mid g_{0},h_{0}),$ where $c_{0},d_{0},e_{0},f_{0},g_{0},h_{0}$ are hyper-parameters. To facilitate the understanding of the Bayesian framework, Fig. 2 shows the hierarchical Bayesian graph and the probability density functions of related random variables.

Subsequently, it is able to write the joint distribution of all variables, viz.,

[TABLE]

Our aim is to obtain the posterior distribution of $\bm{a},\bm{w},\bm{\lambda},\tau$ and $\rho$. Based on the Bayes theorem, we have

[TABLE]

However, the margin distribution

[TABLE]

is computationally infeasible, because our model is complicated. We have to seek an alternative so as to dispense with the computation of $p\left(\bm{y}\right)$. In the literature of Bayesian methods [18, 40, 4], there are mainly two types of algorithms to deal with the inference of complex posterior distributions. One is MCMC sampling which approximates the posterior distribution $p\left(\bm{a},\bm{w},\bm{\lambda},\tau,\rho\mid\bm{y}\right)$ via iteratively drawing samples from the full conditional distributions of each variable. Although MCMC sampling behaves very well in many cases, it is time-consuming when the number of unknown variables is large (i.e., large $N$). The other one is approximation-based techniques [5, 11, 24] such as variational Bayes or expectation propagation (EP) which works by directly utilizing another easily estimated distribution to approximate the desired posterior distribution. The prominent advantage of these methods is their good performance at the much lower computational cost.

Here, we employ a variational Bayesian method to infer this model. The basic idea is to use a variational distribution $q(\bm{a},\bm{w},\bm{\lambda},\tau,\rho)$ to approximate the posterior one. The Kullback-Leibler (KL) divergence is utilized to measure the difference between two distributions [5]. Hence, the original problem is converted into the following optimization issue, namely,

[TABLE]

Specifically, we hypothesize that the variational distributions for each item are independent, that is,

[TABLE]

Then, the optimal variational posterior distribution can be obtained by the Theorem 1.

Theorem 1

The optimal variational posterior distributions of our model are

[TABLE]

where

[TABLE]

{pf}

Please see Appendix A for details.

3.3 Algorithm and implementation details

In this part, we describe how to apply our model to network reconstruction tasks. Algorithm 1 lists the main steps of the inference process of our model. In the non-informative fashion [4], the hyper-parameters are initialized as shown in line 1. Algorithm 2 summarizes the workflow of Bayesian complex network reconstruction, that is, repeatedly computing the response vector and the design matrix and then carrying out Algorithm 1 until all nodes’ structures are recovered. At last, the output of Algorithm 2, $\hat{\bm{W}}$, is the final reconstructed network. In what follows, we abbreviate our method (Variational Bayesian Reconstruction) as VBR.

4 Experiments

In this section, we will carry out experiments to study the behavior of our model. The source code of VBR is available at https://xsxjtu.github.io/Projects/VBR/index.html. Lasso with 5-fold CV to select its tuning parameter is used as the benchmark algorithm and it is implemented by the built-in function lasso in Statistics and Machine Learning Toolbox. All the experiments are carried out on a computer with Intel Core CPU 3.60 GHz, 8.00 GB RAM and Windows 10 (64-bit) system.

To evaluate the reconstruction accuracy of a method, two metrics, ${\rm TPR}$ (true positive rate) and ${\rm TNR}$ (true negative rate), are employed. They are defined as

[TABLE]

and

[TABLE]

respectively. As a matter of fact, TPR is the proportion that existed edges are correctly identified while TNR is the proportion that non-existed edges are correctly excluded. Then, error of connection strength is evaluated by

[TABLE]

In the meanwhile, we also utilize the computational time (in seconds) to compare the efficiency of each algorithm.

In the following subsections, we will carry out five groups of experiments. In the first experiment, we study the performance of VBR and lasso on Barabasi-Albert (BA) [2] and Watts-Strogatz (WS) [33] models. BA and WS are typical scale-free and small-world network models, respectively. Generally speaking, small-world networks have large cluster coefficient and small path length; the degree of scale-free networks follows the power-law distribution. Since the exponential-law coefficient of BA networks is a fixed value, in the second experiment, the random scale-free network model with adjustable exponential-law coefficient is employed. As for the third experiment, it aims at investigating the relationship between execution time and network scale. In the forth experiment, we apply VBR and lasso to four real-world networks rather than simulated ones. In the last experiment, we attempt to reconstruct a stock network based the opening price data.

4.1 Experiment 1: Performance on BA and WS networks

In this first case, we mainly compare the performance of VBR and lasso with regard to reconstruction accuracy and running speed on BA and WS networks with ECT and communication dynamics. The aim is to study their behaviors in the cases with different scales of noise. The detailed experimental settings can be found in Appendix C.

The results are shown in Figs. 3 and 4. From left to right, the mean values of TPR, TNR, Error and the time consumed by each algorithm are plotted as a function of $\sigma$, respectively. In each panel, the bar indicates the mean plus or minus one standard deviation. As a matter of fact, it is not surprising that the performance of both VBR and lasso weakens as $\sigma$ increases. Additionally, the following conclusions can be drawn: (1) For lasso, the fluctuation of TNR and Error curves is much greater. So, compared with VBR, lasso is more sensitive to the scale of noise $\sigma$. The reason is that, VBR uses the full Bayesian inference and the variance of noise is directly modeled by the latent variable $\tau$, while lasso does not consider this factor. Besides, this phenomenon is partially caused by lasso’s biased estimation property. (2) It is found that in most cases the standard deviation of lasso is higher than that of VBR. It means that lasso may be instable. (3) As for computational time, VBR takes a great advantage over lasso and its speed is very robust to $\sigma$. However, the computational time of lasso dramatically increases as $\sigma$ becomes larger. In general, VBR consumes around 0.05-0.2s and lasso takes 2-21s.

As a result, in the context of BA and WS networks, VBR outperforms lasso in both reconstruction accuracy and running speed.

4.2 Experiment 2: Performance on random scale-free networks

The previous investigations have shown that more heterogeneous networks are harder to be reconstructed [13, 36]. Note that the degree of a scale-free network follows the scale-free distribution, i.e., $p({\rm degree}=k)\propto k^{\gamma},$ where exponential-law coefficient $\gamma(<0)$ typically ranges from $-3$ to $-2$. In general, the smaller $|\gamma|$ is, the more heterogeneous a network is. Specifically, the exponential-law coefficient of a BA network is $-3$. In this part, we aim to study the methods’ behaviors on scale-free networks with adjustable $\gamma$. The detailed experimental settings can be found in Appendix C.

Fig. 5 reports the results. Firstly, it is shown that VBR outperforms lasso with different exponential-law coefficients in terms of all metrics. Secondly, both VBR and lasso perform worse as $|\gamma|$ becomes smaller. This phenomenon coincides with the former investigations [13, 36]. However, Fig. 5 shows that the pace of change for lasso is greater than that for VBR. For example, with communication dynamics, the TNR curve of lasso decreases significantly faster when $|\gamma|$ is going smaller. It indicates lasso tends to mistakenly identify null edges as active ones if the scale-free network is more heterogeneous. In conclusion, VBR is more robust to $|\gamma|$ than lasso.

4.3 Experiment 3: Running speed versus network scale

Here, we conduct experiments on relatively large-scale complex networks. Because lasso with CV is extremely slow in this case, we only implement VBR. And, the aim is to study how the running speed of VBR varies with the growing of network scale.

Fig. 6 depicts the execution time of VBR versus $N$, an index indicating the network scale. It can be found in Fig. 6 that the running speed of VBR grows faster than network scale. The main crux lies in the update of $\bm{w}$, which requires the computation of the inverse of an $N$-dimensional matrix. In general, its computation complexity is $O(N^{3})$. Note that it is able to avoid the intensive computation if we assume that $w_{j}(j=1,2,\cdots,N)$ are independent in Eq. (12). However, we find that, in our experiments, this strategy is very likely to be tapped into the local minimum, and the values of TPR and TNR are small. Unfortunately, there is still no efficient approach to reduce the computation complexity and it is extremely interesting to study this issue in the future. At last, it is worth pointing out that the computation complexity of lasso is also $O(N^{3})$ if it is solved by alternating direction multipliers method.

4.4 Experiment 4: Performance on real-world networks

Here, we will conduct experiments on four real-world networks. Since there is no available network with the observed design matrix and response vector, we still simulate them in this subsection. The aim of this experiment is to study the performance of VBR and lasso on real networks rather than the simulated ones (that is, BA, WS and SF networks). HB494BUS and HB1138BUS are available at https://www.cise.ufl.edu/research/sparse/matrices/HB/index.html; Jazz and Karate are available at http://konect.uni-koblenz.de/networks/. In this experiment, we did not include noise item and set $\sigma=0$.

The results reported in Table 2 show that VBR always achieves higher TNR and lower Error than lasso. The values of TPR also manifest the superiority of VBR over lasso. Meanwhile, VBR is very efficient with regard to the execution time. For example, on HB1138BUS, lasso takes about one hour to accomplish the reconstruction task. As for our method VBR, it consumes only less than half an hour to recover the whole network. In summary, the experiments conducted with real-world topology still reveals the superiority of VBR over lasso in terms of all the metrics.

4.5 An empirical study on stock network

In the last experiment, we apply the reconstruction algorithm to a stock market. Tse et al. state that “the fluctuations of stock prices are not independent, but are highly inter-coupled with strong correlations with the business sectors and industries to which the stocks belong” [23]. The financial economists have made attempts to analyze the stock market by means of methodologies of complex network. We assume that the price of a stock is the combination of the prices of other stocks and the correlation strengths between stocks do not change in a short period, that is,

[TABLE]

where $s_{it}$ denotes the price of the $i^{\rm th}$ stock at time $t$, binary variable $a_{ik}$ denotes whether the $k^{\rm th}$ stock affects the $i^{\rm th}$ stock, and $w_{ik}$ is the strength of influence. And $\epsilon_{t}$ is the factor that is independent on $s_{kt}(k=1,2,\cdots,i-1,i+1,\cdots,N)$.

To reconstruct the stock network, we collect the data of stocks in Shanghai Stock Exchange and Shenzhen Stock Exchange from 2018-11-08 to 2019-09-27. It contains the opening prices of 50 stocks during 212 trading days. Each stock is assigned with a label $l_{i}$, that is, the industry to which it belongs. In this dataset, there are five industries, including agriculture, finance & insurance, information technology, real estate and transportation. Since there is no ground truth for the network topology, it is hard to evaluate the performances of VBR and lasso. Based on the prior knowledge that the stocks in the same industry tend to correlate with each other [23], we define the cohesion index (CI) for each node as follows,

[TABLE]

The numerator (or denominator) of CI counts the number of linkages between node $i$ and the one that belongs to the same (or different) industry. Therefore, greater CI indicates that the node tends to connect with ones from the same industry, and, in other words, the reconstructed network is better. Besides CI, we apply a state-of-the-art community detection algorithm, nonnegative matrix factorization, to the reconstructed network, and then use normalized mutual information (NMI) to assess whether the reconstructed network satisfies the prior knowledge. The NMI is defined by

[TABLE]

where $I(\hat{\bm{l}};\bm{l})=\sum_{\hat{l}\in\mathcal{L}}\sum_{l\in\mathcal{L}}p(\hat{l},l)\log[p(\hat{l},l)/p(\hat{l})p(l)]$ denotes the mutual information between the true industry label $\bm{l}$ and the cluster label $\hat{\bm{l}}$, $\mathcal{L}$ is the label space, and $H(\bm{l})=-\sum_{l\in\mathcal{L}}p(l)\log p(l)$ denotes the entropy. Note that the initialization of nonnegative matrix factorization may affect the result of community detection. To exclude this randomness, we apply nonnegative matrix factorization 100 times and compute the average NMI. The algorithm with larger NMI is better.

Fig. 7 displays the reconstructed stock networks. Each node is colored and labeled by its industry and CI value, respectively. The edge is colored by orange if it connects two nodes from the same industry and gray otherwise. The mean CI of VBR and lasso are 0.2977 and 0.2825, respectively. The NMI of VBR and lasso are 0.2722 and 0.2138, respectively. Both metrics demonstrate that the network recovered by VBR is better in the sense of whether the reconstructed network satisfies the prior knowledge.

5 Conclusion and future works

In this paper, we propose a general framework based on Bayesian statistics to reconstruct weighted complex networks. By reformulating the task as a series of regression problems, prior distributions are assigned to the unknown parameters. To efficiently infer from their posterior distributions, a variational Bayesian method is employed. Compared with lasso, the novel method does not need the fine tuning of its associated parameters. The experiments conducted with both synthetic and real networks show that, in the presence of noise, our method VBR outperforms lasso with regard to both reconstruction accuracy and running speed.

Finally, we would like to point out the future work. (a) As shown in Eq. (14), the update of $\bm{w}$ requires the computation of the inverse of an $N\times N$ matrix. The cost will be unacceptable if $N$ (that is, the number of nodes) is large. At the same time, our experiments show that the lasso based methods also suffer from the similar drawback. With the development of big data era, designing a fast network reconstruction algorithm is one of the key future works. (b) To the best of our knowledge, all papers evaluate the effectiveness of algorithm by simulated data (the design matrix and the response vector). To make the comparisons fairer, it is necessary to create and release a dataset with the observed design matrix, response vector, and the ground truth network topology as a benchmark. (c) Last but not least, both VBR and Lasso cannot handle networks with loops. It is worth paying attention to this problem.

Appendix A Proof of theorem 1

In this part, we provide more details about the variational inference (VI), which includes a brief introduction of VI and how the variational distributions for the items shown in Theorem 1 are derived.

A.1 Variational inference

At first, we show how to use VI to infer the optimal solution for a general model. Let $\bm{z}$ denote all the variables to be inferred, where $z_{j}$ is the $j^{\rm th}$ variable; and $\bm{y}$ represents the observed data. According to the main principle of VI, the optimal variational distribution is given by

[TABLE]

According to the definition of KL divergence, we have

[TABLE]

We define the evidence lower bound (ELBO) as ${\rm ELBO}=\mathbb{E}_{q(\bm{z})}[\log p(\bm{z},\bm{y})-\log q(\bm{z})]$, where $\mathbb{E}_{q(\bm{z})}$ represents the expectation operator with regard to (w.r.t.) variational distribution $q(\bm{z})$. Note that Eq. (19) can be reformulated as

[TABLE]

Because $\log p(\bm{y})$ is a constant which does not depend on the variational distribution $q(\bm{z})$, minimizing KL divergence is thus equivalent to maximizing ELBO. Because we assume variational distributions are independent, ELBO can be rewritten as

[TABLE]

where $\mathbb{E}_{-q(z_{j})}$ represents the expectation operator w.r.t. all variational distributions but $q(z_{j})$. The ${\rm const.}$ refers to all items that do not depend on $z_{j}$. Hence, the optimal variational distribution for $z_{j}$ should satisfy

[TABLE]

Because of this property, the variational distribution $q(\bm{z})$ can be efficiently attained by coordinately updating $q(z_{j})$. Remark that $p(\bm{z},\bm{y})$ is the joint distribution. For our model, it refers to Eq. (8). In the next, we show how to acquire the optimal solution to our model.

A.2 Inference of $\bm{w}$

According to Eq. (22), we have

[TABLE]

From Eq. (23), it can be seen that the $q(\bm{w})$ is still a Gaussian distribution. In what follows, we denote it by $\mathcal{N}(\bm{w}\mid\bm{\mu},\bm{\Sigma})$ with

[TABLE]

It should be mentioned that in $\bm{\mu}$ and $\bm{\Sigma}$, the expectations w.r.t. their corresponding variational distribution will be computed in subsection “The expectation computations”.

A.3 Inference of $\lambda_{j}$

According to Eq. (22), there is

[TABLE]

Therefore, $q(\lambda_{j})$ is still a Gamma distribution. In what follows, we denote it by ${\rm Gamma}(g_{j},h_{j})$, where

[TABLE]

A.4 Inference of $\tau$

As for $\tau$, we have

[TABLE]

according to Eq. (22). Hence, $q(\tau)$ is still a Gamma distribution. In what follows, we denote it by ${\rm Gamma}(\tau\mid c,d)$, where

[TABLE]

A.5 Inference of $a_{j}$

To facilitate the derivation, it is worthy mentioning one fact that, if $x\sim{\rm Bernoulli}(\rho)$, then $x^{2}\sim{\rm Bernoulli}(\rho)$. That is, $x$ and $x^{2}$ have the identical distribution. According to Eq. (22), the optimal variational distribution for $a_{j}$ should satisfy

[TABLE]

Here, $X_{j}$ represents the $j^{\rm th}$ column of matrix $\bm{X}$ and $\bm{r}_{j}\equiv\bm{y}-\sum_{n\neq j}X_{n}a_{n}$. Obviously, $q(a_{j})$ is still a Bernoulli distribution. In what follows, we denote it by ${\rm Bernoulli}(\theta_{j})$ with

[TABLE]

A.6 Inference of $\rho$

Similar to Eq. (21), the optimal variational distribution for $\rho$ should satisfy

[TABLE]

Therefore, $q(\rho)$ is still a Beta distribution. In the current paper, we denote it by ${\rm Beta}(e,f)$, where

[TABLE]

A.7 The expectation computations

So far, we have attained the variational distributions for all (hyper)parameters. However, there are some expectations remaining unsolved. In this part, we show how these expectations can be computed. For simplicity, the subscript ‘$q$’ is omitted from $\mathbb{E}_{q}$ in the following discussions.

$\mathbb{E}[\tau]=c/d$. 2. 2.

$\mathbb{E}[\bm{a}]=\bm{\theta}=(\theta_{1},\cdots,\theta_{N})^{\rm T}$. 3. 3.

$\mathbb{E}[\mathbb{D}(\bm{\lambda})]=\mathbb{D}(\bm{g}/\bm{h})=\mathbb{D}(g_{1}/h_{1},\cdots,g_{N}/h_{N})$. 4. 4.

Since $\mathbb{E}[a_{j}a_{n}]=\theta_{j}\theta_{n}(j\neq n)$ and $\mathbb{E}[a_{j}^{2}]=\theta_{j}=\theta_{j}^{2}+\theta_{j}(1-\theta_{j})$, we can thus obtain that $\mathbb{E}[\bm{a}\bm{a}^{\rm T}]=\bm{\Omega}=\bm{\theta}\bm{\theta}^{\rm T}+\mathbb{D}(\bm{\theta})\odot(\bm{I}_{N}-\mathbb{D}(\bm{\theta}))$, where $\mathbb{D}(\bm{\theta})={\rm diag}(\bm{\theta})$ and $\bm{I}_{N}$ is an $N$-dimensional identity matrix. 5. 5.

Note that there is $\bm{w}\sim\mathcal{N}(\bm{w}\mid\bm{\mu},\bm{\Sigma})$. Hence, we have $(\bm{w}-\bm{\mu})(\bm{w}-\bm{\mu})^{\rm T}\sim\mathrm{Wishart}(1,\bm{\Sigma})$, where 1 denotes the degree of freedom and $\bm{\Sigma}$ is the location parameter of the Wishart distribution. Then, there is

[TABLE]

Furthermore, $\mathbb{E}[w_{j}^{2}]=\Sigma_{jj}+\mu_{j}^{2}$. 6. 6.

According to 4, we have $\mathbb{E}[\mathbb{D}(\bm{a})\bm{X}^{\rm T}\bm{X}\mathbb{D}(\bm{a})]=\mathbb{E}\left[\left(\bm{X}^{\rm T}\bm{X}\right)\odot\left(\bm{a}\bm{a}^{\rm T}\right)\right]=\left(\bm{X}^{\rm T}\bm{X}\right)\odot\bm{\Omega}$. 7. 7.

According to 5 and 6, we have

[TABLE] 8. 8.

According to 5 and 6, we have

[TABLE] 9. 9.

$\mathbb{E}[\bm{r}_{j}]=\bm{y}-\sum_{n\neq j}X_{n}\mu_{n}\theta_{n}$; 10. 10.

$\mathbb{E}\left[\log\frac{\rho}{1-\rho}\right]=\psi(e)-\psi(e+f)-\psi(f)+\psi(e+f)=\psi(e)-\psi(f)$, where $\psi(x)$ denotes the digamma function defined as the logarithmic derivative of the gamma function.

By plugging these expectations into the variational distributions derived in previous subsections, and the final solutions shown in Theorem 1 can be easily obtained.

Appendix B Communication dynamics

In this subsection, we introduce how to apply VBR to communication dynamics described in literature [10]. The communication dynamics is used to capture communications in populations via phones or Emails. At time $t_{m}$, individual $i$ contact one of its neighbors $j$ with probability $w_{ji}$ by sending data packets. In this period, the total incoming flux of $i$ is

[TABLE]

where $o_{j,t_{m}}$ is the total outgoing flux from $j$ to its neighbors at time $t_{m}$ and $\sum_{j=1}^{N}a_{ji}w_{ji}=1$. Note the total outgoing flux fluctuates with time. Therefore, we have

[TABLE]

Obviously, VBR can be employed to deal with this case.

Appendix C Experimental Settings

C.1 Experiment 1

We mimic the electrical current transportation (ECT) on a network $\bm{W}$ that is described by following equation,

[TABLE]

where the node’s voltage is generated by alternating current $V_{i}=\bar{V}\sin[(\omega+\Delta\omega_{i})t]$ with $\bar{V}=1$, $\omega=10^{3}$ and $\Delta\omega_{i}$ and the resistance $w_{ij}$ are uniformly sampled from [0, 20] and [2, 3], respectively.

As for the communication dynamics on a network $\bm{W}$ that is described by Eq. (33), the communication probability $w_{ji}$ is uniformly sampled from [0,1], and then is normalized such that $\sum_{j=1}^{N}a_{ji}w_{ji}=1$. The total outgoing flux $o_{j,t_{m}}$ is uniformly sampled from [0,20].

The variance of the noise item is $\sigma^{2}$. In the experiments, we consider $\sigma=0.1,0.2,\cdots,1.0$. For each case, the experiment was conducted 100 times by generating different data. The network size $N=50$. For ECT dynamics, we set the number of observations $M=50$. Since network reconstruction of communication dynamics is harder than that of ECT, we set the number of observations $M=200$ for communication dynamics.

The experiments are conducted on simulated BA and WS networks. The BA networks are generated by the code written by Mathew George 111available at https://ww2.mathworks.cn/matlabcentral/fileexchange/11947-b-a-scale-free-network-generation-and-visualization. The WS networks are generated by the official code of Matlab 222available at https://ww2.mathworks.cn/help/matlab/math/build-watts-strogatz-small-world-graph-model.html?lang=en.

C.2 Experiment 2

The most settings of Experiment 2 are the same to those of Experiment 1. However, in this part, we set $\sigma=0.1$. The random scale-free network with exponential-law coefficient $\gamma$ is generated by mexGraphCreateRandomGraph, a function of Complex Networks Package 333available at http://www.levmuchnik.net/Content/Networks/ComplexNetworksPackage.html. This function is able to generate a graph of given size and with the given node’s degree distribution. In the experiments, $\gamma=-2,-2.2,\cdots,-3$. The number of node is $N=100$.

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

The authors would like to thank the editor, the associate editor and anonymous reviewers for their useful suggestions which greatly helped to improve the paper. The research of S. Xu is supported by the Fundamental Research Funds for the Central Universities [grant number xzy022019059]. The research of C.X. Zhang is supported by the National Natural Science Foundation of China [grant number 11671317] and the National Key Research and Development Program of China [grant number 2018AAA0102201]. The research of P. Wang is supported by the National Natural Science Foundation of China [grant number 61773153], the Supporting Plan for Scientific and Technological Innovative Talents in Universities of Henan Province [grant number 20HASTIT025], and the Training Plan of Young Key Teachers in Colleges and Universities of Henan Province [grant number 2018GGJS021]. The research of J.S. Zhang is supported by the National Key Research and Development Program of China [grant number 2018YFC0809001], and the National Natural Science Foundation of China [grant number 61976174].