1 Introduction↩︎

Identifying network dynamics is a critical approach to understand the mechanism of real-world social systems. The network dynamics is defined as \[\dot{x_i} = F(x_i) + \sum_{j=1}^N A_{ij}G(x_i, x_j), \label{netdyns}\tag{1}\] where \(x_i\in \mathbb{R}^d, i = 1, \dots, N\) are the \(d\)-dimensional activities of \(N\) nodes. The function \(F\) is the self-dynamics of node \(i\), and the function \(G\) is the interactive dynamics of node \(i\) and node \(j\). \(A\) is the weighted adjacency matrix, and \(A_{ij} \geq 0\). The interactive property is fundamental and widely existing in many natural physical systems [1]–[5]. By understanding the underlying mechanism, one can further investigate the resilience and synchronization of social systems, predict long-term evaluation of social behaviors, and then guide social regulation [6], [7].

To identify the network dynamics, [8] propose a two-phase inference algorithm. The algorithm assumes that the network \(A\) is known, and then sparsely identifies the basis matrix of \(F\) and \(G\) that composed of two groups of basis functions. By contrast, [9] propose a two-step prediction algorithm, that assumes \(F\) and \(G\) are known, but \(A\) is unknown, which is exactly the opposite of the two-phase inference algorithm. Similarly, [10] propose an algorithm to identify the network \(A\) for dynamics of Rulkov maps. Thus, these algorithms must either know \(F\) and \(G\), or know \(A\) in prior, but most natural dynamics do not come with such knowledge.

Then, another type of algorithm, e.g., SINDy [11], tries to sparsely identify the nonlinear dynamics without such knowledge. Then, [12] expand that algorithm to the identification of partial differential equations. [13] sparsely identify the nonlinear dynamics with Bayesian principle. Moreover, there are also some research works [14]–[16] that uses artificial neural networks to universally approximate nonlinear dynamics, but the artificial neural networks lack interpretability for the system mechanism. And more importantly, this makes the functional space too extensive to search network dynamics efficiently.

Thus, to identify the network dynamics without any knowledge, and then interpret the network system mechanism, Alternating Sparse Identification of Network Dynamics (ASIND) algorithm is proposed here. The proposed algorithm focuses on the sparse identification of network dynamics, and does not assume the known \(F\), \(G\) and \(A\) in prior. It parameterizes \(F\), \(G\) and \(A\) together, and identifies \(A\) firstly, then identifies \(F\), \(G\) afterwards, and iterate this process until convergence.

Thus, the main contributions of this work are as follows:

1. ASIND algorithm is proposed to identify network dynamics without the knowledge of \(F\), \(G\) and \(A\).

2. Experiments are conducted to validate the state-of-the-art performance on the network dynamics identification, compared to the popular SINDy algorithm [11]. Moreover, the results also show that the network structure \(A\) is weakly identifiable, that means different networks are likely to generate highly similar trajectories.

2 Alternating Sparse Identification of Network Dynamics↩︎

The sparse identification issue of network dynamics in Eq. (1 ) can be defined as \[\begin{align} &\min_{w, A} \sum_{i=1}^N \sum_{m=1}^{M_1 + M_2} \left\Vert w_{im}\right\Vert_1 + \sum_{i=1}^N \sum_{j=1}^N \left\Vert A_{ij}\right\Vert_1 \\ \text{s.t.}~&\dot{x}_i = \sum_{m=1}^{M_1} w_{im} F_m(x_i) + \sum_{m=1}^{M_2} w_{im} \sum_{j=1}^N A_{ij}G_m(x_i, x_j), \\ & A_{ij} \geq 0, i = 1, \dots, N, j = 1, \dots, N. \label{original} \end{align}\tag{2}\] Here, for each node \(i\), \(F_m, m = 1, \dots, M_1\) and \(G_m, m = 1, \dots, M_2,\) are the basis function, and the \(w_{im}, m = 1, \dots, M_1 + M_2\) are the corresponding coefficients. Moreover, \(\left\Vert\cdot\right\Vert_1\) is the \(l_1\)-normalization operator that enables the \(w_{im}\) and \(A_{ij}\) to be sparse.

Then, Eq. (2 ) can be rewritten as \[\begin{align} &\min_{w, A, \lambda} \mathcal{L}(w, A, \lambda) \\ \text{s.t.}~&A_{ij} \geq 0, i = 1, \dots, N, j = 1, \dots, N. \label{new} \end{align}\tag{3}\] where, \[\small \begin{align} &\mathcal{L}(w, A, \lambda) \\ = &\sum_{i=1}^N\sum_{m=1}^{M_1 + M_2} ||w_{im}||_1 + \sum_{i=1}^N \sum_{j=1}^N ||A_{ij}||_1 \\ & + \sum_{i=1}^N \lambda_i \left(\dot{x}_i - \sum_{m=1}^{M_1} w_{im} F_m(x_i) - \sum_{m=1}^{M_2} w_{im} \sum_{j=1}^N A_{ij}G_m(x_i, x_j) \right) \\ & + \frac{\rho}{2} \sum_{i=1}^N \left\Vert\dot{x}_i - \sum_{m=1}^{M_1} w_{im} F_m(x_i) - \sum_{m=1}^{M_2} w_{im} \sum_{j=1}^N A_{ij}G_m(x_i, x_j) \right\Vert_2^2. \end{align}\] Here, \(\left\Vert\cdot\right\Vert_2\) is the \(l_2\)-normalization operator for augmented Lagrangian term. \(\lambda_i\) is the Lagrange multiplier, and \(\rho\) is the augmented Lagrangian penalty parameter. The Eq. (2 ) and Eq. (3 ) have the same optimal solution, and the Eq. (3 ) can be alternatively solved by the following steps: \[\begin{align} &A^{(k+1)} = \mathop{\arg\min}\limits_{A} \mathcal{L}\left(w^{(k)}, A, \lambda^{(k)}\right), \\ \text{s.t.}~&A_{ij} \geq 0, i = 1, \dots, N, j = 1, \dots, N, \end{align} \label{Ak}\tag{4}\] and then, \[w^{(k+1)} = \mathop{\arg\min}\limits_{w} \mathcal{L}\left(w, A^{(k+1)}, \lambda^{(k)}\right), \label{wk}\tag{5}\] and then, \[\lambda^{(k+1)} = \lambda^{(k)} + \alpha \times \frac{\partial \mathcal{L}\left( w^{(k+1)}, A^{(k+1)}, \lambda \right)}{\partial \lambda}, \label{lk}\tag{6}\] where \(\alpha\) is the step size to control the gradient descent of \(\lambda\).

Then, to update \(A^{(k+1)}\) in Eq. (4 ), independently for \(A_{i}^{(k+1)} = \left[ A_{i1}^{(k+1)}, \dots, A_{iN}^{(k+1)} \right], i=1,\dots, N\), \[\small \begin{align} &\mathop{\arg\min}\limits_{A_{i}} \mathcal{L}\left(w^{(k)}, A_{i}, \lambda^{(k)}\right) \\ = & \mathop{\arg\min}\limits_{A_{i}} \sum_{j=1}^N A_{ij} \\ & + \lambda_i^{(k)}\left( \dot{x}_i - \sum_{m=1}^{M_1} w_{im}^{(k)} F_m(x_i) - \sum_{m=1}^{M_2} w_{im}^{(k)} \sum_{j=1}^N A_{ij}G_m(x_i, x_j) \right) \\ & + \frac{\rho}{2} \left\Vert \dot{x}_i - \sum_{m=1}^{M_1} w_{im}^{(k)} F_m(x_i) - \sum_{m=1}^{M_2} w_{im}^{(k)} \sum_{j=1}^N A_{ij}G_m(x_i, x_j) \right\Vert_2^2. \\ &\text{s.t.}~A_{ij} \geq 0, j = 1, \dots, N. \end{align} \label{update95A}\tag{7}\]

Then, to update \(w^{(k+1)}\) in Eq. (5 ), independently for \(w_{i}^{(k+1)} = \left[ w_{i1}^{(k+1)}, \dots, w_{i(M_1+M_2)}^{(k+1)} \right], i=1,\dots, N\), \[\fontsize{7.9pt}{10pt}\selectfont \begin{align} &\mathop{\arg\min}\limits_{w_i} \mathcal{L}\left( w_i, A^{(k+1)}, \lambda^{(k)} \right) \\ =& \mathop{\arg\min}\limits_{w_i} \sum_{m=1}^{M_1 + M_2} \left\Vert w_{im} \right\Vert_1 \\ & + \lambda_i^{(k)}\left( \dot{x}_i - \sum_{m=1}^{M_1} w_{im} F_m(x_i) - \sum_{m=1}^{M_2} w_{im} \sum_{j=1}^N A_{ij}^{(k+1)} G_m(x_i, x_j) \right) \\ & + \frac{\rho}{2} \left\Vert \dot{x}_i - \sum_{m=1}^{M_1} w_{im} F_m(x_i) - \sum_{m=1}^{M_2} w_{im} \sum_{j=1}^N A_{ij}^{(k+1)} G_m(x_i, x_j) \right\Vert_2^2. \end{align} \label{update95w}\tag{8}\] To solve the Eq. (8 ), we set \(u_i = \left[ u_{i1}, \dots, u_{i(M_1 + M_2)} \right]\), and \(u_{im} = \max(0, w_{im}), m = 1, \dots, M_1 + M_2\). And then, we set \(v_i = \left[ v_{i1}, \dots, v_{i(M_1 + M_2)} \right]\), and \(-v_{im} = \min(0, w_{im}), m = 1, \dots, M_1 + M_2\). Thus, it is easy to obtain, \(u_{im} \geq 0, v_{im} \geq 0\), and \(w_{im} = \max(0, w_{im}) + \min(0, w_{im}) = u_{im} -v_{im}\). Then, Eq. (8 ) can be rewritten as, \[\begin{align} & \mathop{\arg\min}\limits_{u_i, v_i} \sum_{m=1}^{M_1 + M_2} u_{im} + \sum_{m=1}^{M_1 + M_2} v_{im} \\ & + \lambda_i^{(k)} \bigg[\dot{x}_i - \sum_{m=1}^{M_1} (u_{im} - v_{im}) F_m(x_i) \\ &- \sum_{m=1}^{M_2} (u_{im} - v_{im}) \sum_{j=1}^N A_{ij}^{(k+1)} G_m(x_i, x_j) \bigg] \\ & + \frac{\rho}{2} \bigg\Vert\dot{x}_i - \sum_{m=1}^{M_1} (u_{im} - v_{im}) F_m(x_i) \\ & - \sum_{m=1}^{M_2} (u_{im} - v_{im}) \sum_{j=1}^N A_{ij}^{(k+1)} G_m(x_i, x_j)\bigg\Vert_2^2, \\ & \text{s.t.}~u_{im} \geq 0, v_{im} \geq 0, m = 1, \dots, M_1 + M_2. \end{align} \label{update95z}\tag{9}\] Then the \(l_1\)-normalization term in Eq. (8 ) is released, and Eq. (8 ) can be equivalently solved by Eq. (9 ). If \(u_i\) and \(v_i\) are both obtained, \(w_i\) can be ontained by \(w_i = u_i - v_i\). It is trivial to see, Eq. (7 ) and Eq. (9 ) are both quadratic programming problem, and they have the same form for \(A_i\), \(u_i\) and \(v_i\). Thus, we can use the same algorithm to solve them.

Then, to update \(\lambda^{(k+1)}\) in Eq. (6 ), we can obtain that \[\small \begin{align} &\frac{\partial \mathcal{L}(w^{(k+1)}, A^{(k+1)}, \lambda)}{\partial \lambda} \\ = &\dot{x}_i - \sum_{m=1}^{M_1} w_{im}^{(k+1)} F_m(x_i) - \sum_{m=1}^{M_2} w_{im}^{(k+1)} \sum_{j=1}^N A_{ij}^{(k+1)}G_m(x_i, x_j), \end{align}\] and we can use it to replace the corresponding term in Eq. (6 ). Thus, due to the convexity of Eq. (3 ), the unique solution can be obtained, and it is the same solution to Eq. (2 ).

3 Prediction Experiments↩︎

To test the prediction performance of ASIND algorithm, this work conducts multiple experiments to predict network dynamics over 100 steps, with different networks and different functions on Eq. (1 ) for universal social systems.

3.1 Models of Dynamics on Networks↩︎

Four universal network dynamics are used here to initialize Eq. (1 ), as shown in Table 1. Among them:

Kuramoto model [17], [18] is used to describe the synchronization between \(N\) coupled oscillators through phase coupling. In the model, each oscillator has its own natural frequency, and the coupling term promotes interactions between oscillators, which leads to the collective behavior. Thus, \(\omega_i\) is the frequency of oscillator \(i\), and \(c>0\) is the coupling coefficient.

Susceptible-Infected-Susceptible (SIS) model [19], [20] is used to describe the synchronization between \(N\) individuals, that the state of individuals switch cyclically between susceptible and infected states. Each individual can recover from the infected state, and later become infected again. The SIS model is widely used to study the epidemic dynamics in populations, thus here \(\delta_i > 0\) is recovery rate, and \(\gamma_i \geq 0\) is infection rate.

Lotka-Volterra (LV) model [21] is used to describe the dynamical evolution between \(N\) species competing for shared resources within an ecosystem. The growth of each species is influenced by its own population density and the population density of other competing species, which leads to coexistence and competitive exclusion between natural species. Thus, \(\alpha_i > 0\) and \(\theta_i > 0\) control the intrinsic growth rate of species \(i\) together, and \(\gamma_i\) is the competition coefficient.

Michaelis-Menten (MM) model [1], [22], [23] is used to describe the interactions between enzymes and substrates in biological networks, and interpret the signal transmission in metabolic pathways. The model demonstrates the mechanism of steady state behaviors and response characteristics of the biochemical systems. Thus, \(h\) is Hill coefficient, which represents the binding strength between the enzyme and its substrate.

17pt

Moreover, the adjacency \(A\) in Eq. (1 ) is initialized with Erdős-Rényi (ER) graph [24], Watts-Strogatz (WS) network [25] and Barabási-Albert (BA) network [26], respectively. In ER graph, the connection probability \(p=0.1\). In WS network, the average degree \(\left\langle k \right\rangle=4\) and the rewiring probability \(p=0.1\). In BA network, the three control parameters \(\alpha=0.41, \beta=0.54, \gamma=0.05\). The network size is 16.

9pt

3.2 Baseline↩︎

The prediction performance of ASIND is compared to that of the baseline SINDy algorithm [11]. Note that, the two algorithms are both conducted without any knowledge, and can be used for universal network dynamics identification.

SINDy is proposed to sparsely identify the governing equations in nonlinear dynamical systems. The algorithm uses a group of basis functions to make a basis matrix from observation data, and then identifies the weights of the basis functions by \(l_0\)-normalization optimization, to reconstruct the nonlinear dynamics finally. Here, the basis matrix is mainly composed of polynomial basis up to second order.

3.3 Evaluation Indices↩︎

To compare the prediction performance of ASIND algorithm with that of the other baseline, root mean square error (RMSE) and mean absolute percentage error (MAPE) are introduced to evaluate prediction performance, as follows: \[\begin{align} \text{RMSE} &= \sqrt{\frac{1}{T} \sum_{t=1}^T \left(\hat{x}_t - x_t \right)^2}, \\ \text{MAPE} &= \frac{1}{T} \sum_{t=1}^T \left| \frac{\hat{x}_t - x_t}{x_t}\right| \times 100\%, \end{align}\] where \(x_t\in \mathbb{R}^N\) is the observation activities of the \(N\) nodes at time \(t\), and \(\hat{x}_t\) is the corresponding prediction.

Moreover, to quantify and measure the identifiability of network \(A\) in Eq. (1 ), Jaccard index is introduced as follows: \[J(A, \hat{A}) = \frac{A \cap \hat{A}}{A \cup \hat{A}} \times 100\%,\] where \(\hat{A}\) is the estimation of \(A\). Note that, all elements of \(A\) and \(\hat{A}\) are converted to zero (i.e., 0) and non-zero elements (i.e., 1) when calculate Jaccard index.

3.4 Prediction Performance↩︎

As shown in Table 2, ASIND algorithm achieves the state-of-the-art performance on the task of 100-steps prediction, and the error indices are significantly low, especially for the network dynamics including SIS model, LV model and MM model. This means the prediction results of ASIND are close to the ground truths during long-term prediction, and the network dynamics are identified accurately.

Meanwhile, the prediction results of SINDy significantly deviate from the ground truths for 100-steps prediction, that means SINDy cannot identify the network dynamics accurately. However, it is found that SINDy perform well on the Kuramoto model. A possible reason is that high-order polynomial bases can well approximate the sine function, but hardly approximate the complex functions in other network dynamical models.

3.5 Weak Identifiability of Network Structure↩︎

As shown in Table 3, the network structure \(A\) in Eq. (1 ) is weakly identifiable. The ASIND algorithm can be trained on the observation data and then predict 100-steps evaluation accurately, but the Jaccard index \(J(A, \hat{A})\) is still low for different network dynamics. This means, the reconstructed networks and the ground-truth networks are equivalent to represent the network dynamics, and that makes the ground-truth networks difficult to be identified from the observation data.

6pt

4 Conclusion↩︎

ASIND algorithm is proposed here to sparsely identify the network dynamics without the knowledge of \(F\), \(G\) and \(A\). Experiments are conducted to show the state-of-the-art identification and 100-steps prediction performance, compared to the popular SINDy algorithm. Moreover, the results also show that the network structure \(A\) has weak identifiability, that means different networks can generate highly similar trajectories of network dynamics.

Of course, the experiments here are oversimplification for a fair comparison with SINDy, e.g., network size \(N=16\), which is not in line with the real-world systems. In fact, the maximum network size ASIND can support is \(N=100\). If the number goes beyond this point, the performance gradually declines.

References↩︎