# Finite Expression Method for Learning Dynamics on Complex Networks

Zezheng Song Email<EMAIL_ADDRESS>Department of Mathematics, University of
Maryland, College Park, MD 20742, USA  Chunmei Wang Email:
<EMAIL_ADDRESS>Department of Mathematics, University of Florida,
Gainesville, FL 32611, USA  Haizhao Yang Email<EMAIL_ADDRESS>Department of
Mathematics, University of Maryland, College Park, MD 20742, USA  Department
of Computer Science, University of Maryland, College Park, MD 20742, USA

###### Abstract

Complex network data pervades various real-world domains, including physical,
technological, and biological systems. Despite the prevalence of such data,
predicting trends and understanding behavioral patterns in complex systems
remains challenging due to poorly understood underlying mechanisms. While
data-driven methods have made strides in uncovering governing equations from
time series data, efforts to extract physical laws from network data are
limited and often struggle with incomplete or noisy data. To address these
challenges, we introduce a novel approach called the Finite Expression Method
(FEX) and its fast algorithm for this learning problem on complex networks.
FEX represents dynamics on complex networks using binary trees composed of
finite mathematical operators. The nodes within these trees are trained
through a combinatorial optimization process guided by reinforcement learning
techniques. This unique configuration allows FEX to capture complex dynamics
with minimal prior knowledge of the system and a small dictionary of
mathematical operators. Our extensive numerical experiments demonstrate that
FEX excels in accurately identifying dynamics across diverse network
topologies and dynamic behaviors.

Keywords Complex Networks, Network Dynamics, Finite Expression Method, High
Dimensions, Symbolic Learning

##  1 Introduction

Complex network data is prevalent in various real-world domains, spanning from
transportation [44, 62, 72, 74] to biology [13, 10, 43] and epidemiology [49,
15, 55]. For instance, proteins frequently interact with one another to
execute a wide range of functions. As with many complex systems, comprehending
the underlying dynamics within these networks is crucial for analyzing and
predicting their behaviors. In practice, time-series data for individual nodes
within a network are often accessible. Still, extracting the dynamic
interactions among these nodes to gain insights into the system’s behavior
proves to be a challenging endeavor. Consequently, the physical laws governing
the majority of real-world complex networks remain elusive and underexplored.

In recent years, machine learning has emerged as a promising avenue for data-
driven discovery of physical laws. Notably, deep learning has demonstrated its
potential to identify governing partial differential equations (PDEs) from
noisy data [50, 51, 48, 65, 39, 40, 17, 57, 3, 8, 28]. The success of these
methods is attributed to the expressive power of deep learning models in
capturing relationships among different variables. However, there has been
limited focus on the discovery of physical laws within complex networks.
Discovering dynamics on network data involves various challenges:

  * •

Complexity and scale: Real-world network data can be extremely large and
complex, making analysis computationally expensive. In the full interaction
model with $N$ nodes in the network, each node interacts with every other
particle, leading to $N(N-1)$ interaction terms that need to be computed at
each time, resulting a $O(N^{2})$ complexity.

  * •

Data quality and availability: The quality of network data can vary greatly.
Issues like missing data, noise, and topological inaccuracies can
significantly impact the analysis.

In Subsection 1.1, we summarize some of the recent work on learning dynamics
from network data, and we also introduce some of the recent progress on fast
algorithms. In Subsection 1.2, we introduce the overall idea of our proposed
method and how it tackles these challenges.

###  1.1 Related work

In this Subsection, we introduce related work on identifying governing
equations on network data and fast algorithms for scientific computing.

####  1.1.1 Discoverning physical laws on network

SINDy: The Sparse Identification of Nonlinear Dynamics (SINDy) [11] algorithm
aims to extract underlying governing equations from data in a interpretable
and efficient way, with theoretical convergence analysis [75]. Taking time-
series data as inputs, SINDy infers the underlying physical laws of the system
by using Sequentially Thresholded Least Squares algorithm to choose a
parsimonious model from a large library of candidate functions, usually
including polynomials and trignometric functions. SINDy has been used to
identify models in fluid dynamics [23], control systems [37], and neuroscience
[14], etc. However, it has been shown that the vanilla SINDy performs poorly
in identifying dynamics on network data.

Two-Phase SINDy: In [25], Gao and Yan proposed a two-phase procedure
consisting of global regression and local fine-tuning for graph dynamics
inference. This method works by building two comprehensive function libraries,
$L_{F}$ and $L_{G}$, for self and interaction dynamics. Then, in the first
phase, this method performs a sparse regression on the normalized input data
and the libraries $L_{F}$ and $L_{G}$ to identify potential terms in the
network dynamics. In the second phase, it performs topological sampling to
fine-tune the coefficients of functions from the first phase. This method can
identify the explicit functional forms of the system dynamics, given
sufficient large function libraries. However, users often lack a priori
knowledge of the candidate functions, so they assume a prior solid knowledge
of the system. In addition, this method is sensitive to noisy and low-
resolution data because topological samplings are used in the fine-tuning
stage.

Graph Neural Networks: Graph neural networks (GNNs) are popular models to
learn the dynamics on the networks and make predictions on the systems [54,
73, 22, 24]. These models typically use an encoder neural network to denote
node states and create a latent representation. Following this, Neural
Ordinary Differential Equations (Neural ODEs) are applied on the latent space
to learn the dynamics. Finally, a decoder is used to project the latent
embedding back to the original high-dimensional state. While this approach is
flexible enough to be applied to a wide variety of network data, it
essentially functions as a “black box”, making it difficult to understand the
underlying dynamics. Additionally, it struggles to generalize well to out-of-
distribution data.

ARNI: Algorithm for Revealing Network Interactions (ARNI) [12] is a model-free
method originally proposed to identify direct interactions of particles in a
network data. It works by decomposing the dynamics of each node as a linear
combination of basis functions in pairwise and higher-order interactions with
other nodes in the system, and recast the reconstruction problem into a
mathematical regression problem with grouped variables. Authors of [25]
modified ARNI appropriately to make it capable of identifying network dynamics
instead of only inferring the graph topology.

####  1.1.2 Fast algorithms for scientific computation

In the realm of physical sciences, it is important to develop fast and
efficient algorithms to handle complex systems while keeping a high accuracy.
For example, Jin et al. [36] proposed the Random Batch Methods (RBM), which
splits particles in a physical system into several batches at each time to
reduce the computational cost from $O(N^{2})$ to $O(N)$, where $N$ is the
number of particles. In machine learning, stochastic gradient descent (SGD)
[2, 9, 32] takes a random sample of the entire dataset in each iteration to
compute the gradient and update model parameters. Similar to SGD, stochastic
coordinate descent [46, 47] randomly samples some of the coordinates to update
at each time step. In addition, randomized algorithms have been widely applied
to fast matrix factorizations [18, 31, 53] and Bayesian inference [71, 52].

###  1.2 Our contributions

To tackle the problems seen in previous methods, we propose the finite
expression method (FEX) that represents dynamics on network data by binary
trees, whose nodes are mathematical operators (e.g., sin, cos, exp, etc). FEX
selects the operator of each node via combinatorial optimization (CO) by a
novel reinforcement learning (RL) approach. In particular, our method has
major characteristics. (1) Compared to the widely used SINDy approach and its
variants, FEX does not require a large library of candidate functions. In
contrast, thanks to the binary tree representation, FEX can produce
complicated composition of functions with a few simple mathematical operators
provided by the user. (2) Our model employs a RL approach to perform operator
selection in a data-driven manner (discussed in Section 3). This design has
been demonstrated to promote the expressiveness of our method of identifying
network dynamics. (3) To mitigate the intensive computational burden when
calculating the pair-wise interaction dynamics on network data, we employ RBM
[36] to accelerate the computation while keeping a reasonable accuracy.
Throughout this work, we demonstrate the advantage of our method by applying
it to several synthetic complex network data to extract governing equations.
An overview figure of FEX to infer physical laws on network data is provided
in Fig. 1.

Figure 1: (A) Based on the input data of the time series and adjacency matrix,
FEX performs combinatorial optimization to infer the mathematical structures
of the dynamics and then implements fine-tuning on the top-performing
candidates in the pool to identify the optimal dynamics. (B) FEX implements
combinatorial optimization by training a controller $\chi$ to output
probability mass function for each node in the FEX binary tree.

##  2 Problem setup

###  2.1 Learning dynamics on complex networks

A common model [25, 45] that describes the evolution of the states of the
network of size $N$ is given by

$\frac{\mathrm{d}\mathbf{x}_{i}(t)}{\mathrm{d}t}=\mathbf{F}\left(\mathbf{x}_{i}(t)\right)+\sum_{j=1}^{N}\mathbf{A}_{ij}\mathbf{G}\left(\mathbf{x}_{i}(t),\mathbf{x}_{j}(t)\right),$
(1)

where
$\mathbf{x}_{i}(t)\equiv\left(x_{i,1}(t),\ldots,x_{i,d}(t)\right)^{\top}$
represents node $i$’s $d$-dimensional activity and $i=1,\ldots,N$. The
dynamics of $\mathbf{x}_{i}(t)$ are driven by
$\mathbf{F}(\mathbf{x}_{i})\equiv\left(F_{1}\left(\mathbf{x}_{i}\right),\ldots,F_{d}\left(\mathbf{x}_{i}\right)\right)^{\top}$,
denoting the self-dynamics, and
$\mathbf{G}(\mathbf{x}_{i}(t),\mathbf{x}_{j}(t))\equiv\left(G_{1}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right),\ldots,G_{d}\left(\mathbf{x}_{i},\mathbf{x}_{j}\right)\right)^{\top}$,
representing the interaction dynamics between any pair of nodes
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$. $\mathbf{A}_{ij}$ denotes the $ij$-th
entry in the $N\times N$ adjacency matrix $\mathbf{A}$, where
$\mathbf{A}_{ij}=1$ if there is a connection between node $i$ and node $j$,
and $\mathbf{A}_{ij}=0$ otherwise. Therefore, the dynamics are determined by
self-dynamics, interaction dynamics and the topology of the network. However,
in many real-world complex network systems, we do not have the information of
the explicit forms of $\mathbf{F}$ and $\mathbf{G}$. Therefore, it is
important to distill the functional forms of $\mathbf{F}$ and $\mathbf{G}$
from node activity data, which is the focus of our work.

However, inferring underlying dynamics from real-world time series data
presents a formidable challenge, primarily due to the often unclear or only
partially known governing laws. The endeavor to discover dynamics within a
complex network system introduces an additional layer of complexity,
attributable to the interactions among entities described by the graph
topology. Moreover, numerous oscillator systems exhibit varied coupling
behaviors contingent upon the coupling strength between entities and the
inherent graph topology, a prime example being synchronization phenomena. In
addition, the data in real-world are often sparse, noisy, and come from
incomplete network topology. Hence, to demonstrate the capability of FEX to
discover dynamics on complex networks, we subject FEX to a rigorous
examination against three classical network dynamics: Hindmarsh-Rose (HR,
$d=3$) [34, 64, 70], FitzHugh-Nagumo (FHN, $d=2$) neuronal systems [19], and
coupled heterogeneous Rössler oscillators ($d=3$) [7], where $d$ denotes the
dimension of feature for each node. Regarding to the graph topology, we
demonstrate that FEX is capable of identifying governing equations on various
synthetic networks, including Erdős-Rényi (ER) [27] and scale-free (SF) [1]
networks. Through this extensive examination, we aim to demonstrate the
effectiveness and robustness of FEX in decoding the intricate dynamics
inherent in complex network systems.

###  2.2 FEX for learning dynamics on networks

The motivation behind FEX is to develop a robust methodology for identifying
physical laws from network data. Our model depicts network dynamics using two
fixed-size binary trees: one representing the self-dynamics
$\mathbf{F}(\mathbf{x}_{i})$ while another one representing the interaction
dynamics $\mathbf{G}(\mathbf{x}_{i},\mathbf{x}_{j})$. Each tree node denotes a
simple mathematical operator—either unary (sin, cos, and exp …) or binary
($+,-,\times$…). The aim of the learning process is to identify the most
suitable mathematical operator for each node of the tree, ensuring accurate
representation of the dynamics. The input to FEX algorithm consists of time
series data $\mathbf{x}_{i}$ for each node in the network, and the adjacency
matrix $\mathbf{A}$ representing the network topology, as illustrated by Fig.
2.

Figure 2: Representation of the components of our implementation with FHN
dynamics. (A) In the offline stage, we generate time series data according to
the governing law of the system, potentially with noise or low resolution. (B)
In the online phase, based on the time series data, we use two FEX binary
trees to infer the dynamics of the system.

###  2.3 Learning dynamics from sparse and noisy data

In real-world scenarios, data are inevitably contaminated by noise or by the
incompleteness of network topology. Thus, evaluating the robustness of the
proposed method across diverse scenarios—such as low resolution, observational
noises, and spurious or missing links—is paramount. Moreover, a comparative
analysis with baseline models is conducted to ascertain the relative efficacy
of the proposed method.

Low resolution: Practical constraints or technological limitations often
result in the collection of sparse raw data. To emulate such circumstances, we
simulate nonlinear dynamics with a time step size of $t=0.01$, followed by a
uniform down-sampling of the time series data. The outcomes reveal that FEX
retains a high degree of accuracy, i.e., relative errors of coefficients
identified by FEX and the true coefficients are less than $1\%$, even with
merely $5\%$ or $10\%$ of the original data.

Spurious and missing links: In the real world, network topology may undergo
evolution, making the capture of all node interactions a challenging, if not
unfeasible, task. To assess the robustness of FEX against incomplete topology,
we randomly either add or remove a proportion of links from the authentic
network topology. The results indicate that FEX can accommodate approximately
$15\%$ missing or spurious links, showcasing its robustness in handling
topological inaccuracies.

Observational noise: Observational noise is another common issue in data
measurement. We replicate such noise by introducing Gaussian noises to the
time series data for every node within the network. The noise intensity is
modulated via the signal-to-noise ratio (SNR). The findings underscore that
FEX sustains a high accuracy level, even amidst a 40-dB observational noise
environment. In contrast, other baseline methods such as Two-Phase SINDy [25]
and ARNI [12] identified erroneous and redundant terms with low-accurate
coefficients.

Through this thorough evaluation across varied scenarios and a comparative
analysis with baseline models, we aim to demonstrate the robustness and
stability of FEX in extracting dynamics from real-world, noisy data with
potential topological inaccuracies. Further details are provided in Section 4.

##  3 Methods

###  3.1 Overview of FEX

FEX [42, 63] introduces a novel strategy for solving PDEs and identifying
physical laws from data [35]. FEX represents the system’s dynamics using a
learnable binary tree to identify dynamics on network data. This tree
comprises unary operators (e.g., sin, cos, etc.) and binary operators
($+,-,\times$, etc.). Additionally, each unary operator incorporates learnable
weights and biases, augmenting the expressiveness of the expression. Given a
specific set of operators in the tree, input variables
$\mathbf{x}\in\mathbb{R}^{d}$ are passed through the FEX tree via the leaves.
The complete FEX expression is then derived through a preorder traversal of
the tree. Unlike many symbolic regression (SR) methods that generate equations
autoregressively, FEX allows users to predetermine the tree’s depth and size.
This distinctive feature has been demonstrated to enhance FEX’s accuracy in
solving PDEs and discovering governing equations. Subsequently, with a
suitable loss function, FEX conducts CO using a RL strategy to perform
operator selection, and continuous optimization to refine the parameters
associated with unary operators.

To address this mixed optimization problem, FEX employs a CO technique to
determine the best selection of operators while leveraging a continuous
optimization method to finetune the corresponding parameters of tree nodes.
Regarding the CO, FEX utilizes a RL strategy to transform CO into a continuous
optimization over probability distributions. To achieve this transformation, a
controller network (modeled as a fully connected neural network) is integrated
into FEX, which produces the probability mass function for each operator’s
selection at every tree node. Consequently, selecting the most suitable
operators for tree nodes is recast as the problem of training the optimal
controller network capable of sampling top-performing operators. The
refinement process of identifying the best controller is fundamentally a
continuous optimization task. The controller network’s parameters are refined
through policy gradient techniques to maximize the anticipated reward,
consistent with RL conventions. As a result, the best operators can be
efficiently identified by sampling from the probability distribution of the
learned controller network.

In particular, analyzing dynamics in complex networks brings forth several
challenges, primarily the computational complexity involved. Computing
particle interactions on a graph requires $O(N^{2})$ complexity, where $N$ is
the number of nodes in a network. This is prohibitive for large-scale systems.
We’ve employed the RBM [36] when calculating the loss function to mitigate
this. Simply put, nodes in the network are grouped into random batches and
particle interactions are limited to their respective batches. By integrating
RBM into the CO segment of FEX training, we’ve achieved faster computational
speeds, i.e., $O(N)$ computational complexity, without sacrificing the
precision required to identify the accurate structure of the FEX tree.

###  3.2 Functional space of finite expressions

FEX models the dynamics of network data within a functional space
characterized by a finite set of operators. As such, it is crucial to clearly
define this functional space where the solution resides.

######  Definition 3.1 (Mathematical expression [42]).

A mathematical expression is a combination of symbols, which is well-formed by
syntax and rules and forms a valid function. The symbols include operands
(variables and numbers), operators (e.g., “+”, “sin”, integral, derivative),
brackets, and punctuation.

######  Definition 3.2 ($k$-finite expression [42]).

A mathematical expression is called a $k$-finite expression if the number of
operators in this expression is $k$.

######  Definition 3.3 (Finite expression method).

The finite expression method aims to numerically capture dynamics by using a
limited size expression, so that the produced function closely matches the
actual dynamics.

We denote $\mathbb{S}_{k}$ the functional space that consists of functions
formed by finite expressions with the number of operators less than or equal
to $k$.

###  3.3 The combinatorial optimization problem in FEX

In FEX, the loss functional $\mathcal{L}$ varies based on the problem. For
this study, we adopt the least squares loss functional, as presented in [25],
which is described as follows:

$\mathcal{L}(F,G)=\int_{0}^{T}\left(\left\|\text{FEX\\_Tree}_{F}(\mathbf{x},t)+{\mathbf{A}}\text{FEX\\_Tree}_{G}(\mathbf{x},t)-\dot{\mathbf{x}}(t)\right\|^{2}\right)\mathrm{d}t,$
(2)

where $\mathbf{A}$ is the adjacency matrix of the network, and
$\text{FEX\\_Tree}_{F}(\mathbf{x},t)$ denotes the self-dynamics modeled by a
FEX tree, $\text{FEX\\_Tree}_{G}(\mathbf{x},t)$ denotes the interaction
dynamics by another FEX tree. In FEX, the solution is found by solving the
combinatorial optimization problem

$\displaystyle\min_{F,G\in\mathbb{S}_{\sf FEX}}\mathcal{L}(F,G),$ (3)

where the solution space $\mathbb{S}_{\sf FEX}\subset\mathbb{S}_{k}$ will be
elaborated in Section 3.4.1.

###  3.4 Implementation of FEX

As mentioned previously, FEX’s computational process begins with the formation
of a mathematical binary tree, where each node represents either a unary or a
binary operator. The candidate solution $u$ is derived from evaluating the
function represented by this tree. Subsequently, the CO denoted by (3) is
employed to dynamically choose the best operators for all tree nodes. This
optimization aims to find operators that recover the structure of the genuine
solution.

####  3.4.1 Finite expressions with binary trees

FEX employs a binary tree structure $\mathcal{T}$ to represent finite
expressions. The user defines the sets of potential unary and binary
operators, represented by $\mathbb{U}$ and $\mathbb{B}$. Common choices of
unary and binary operators are respectively,

$\sin,\exp,\text{Id},(\cdot)^{2},\int\cdot\text{d}x_{i},\frac{\partial\cdot}{\partial
x_{i}},\cdots\ \quad{\rm and}\quad+,-,\times,\cdots.$

Every unary operator in the leaf level operates on the inputs in an element-
wise manner, followed by scaling parameters $\alpha_{i}$, where
$i=1,\ldots,d$, and a bias parameter $\beta$. For instance:

$\alpha_{1}\exp(x_{1})+\ldots+\alpha_{d}\exp(x_{d})+\beta.$

Figure 3: Computational rule of a FEX binary tree. Each unary operator is
equipped with a weight $\alpha$ and bias $\beta$ to enhance expressiveness.
For tree depth greater than one (i.e., $L>1$), the computation is implemented
by recursion.

All parameters associated with unary operators are collectively represented by
$\bm{\theta}$. The complete expression is derived through a preorder traversal
of the operator sequence, $\bm{e}$, in the binary tree $\mathcal{T}$.
Consequently, this finite expression is represented as
$u(\bm{x};\mathcal{T},\bm{e},\bm{\theta})$, as a function of input $\bm{x}$.
Given a fixed $\mathcal{T}$, the maximum number of operators remains
constrained and is upper-bounded by a constant, labeled $k_{\mathcal{T}}$. In
FEX, the functional space wherein the dynamics of the system is solved is
defined as:

$\mathbb{S}_{\sf
FEX}=\\{u(\bm{x};\mathcal{T},\bm{e},\bm{\theta})~{}|~{}\bm{e}\in\mathbb{U}\cup\mathbb{B}\\}$

The computation within the binary tree $\mathcal{T}$ proceeds recursively,
beginning at the leaf nodes. Unary operators at these leaf nodes act on the
input $\bm{x}$ in an element-wise manner, with the scaling factor $\alpha$,
transforming the dimension from $\mathbb{R}^{d}$ to $\mathbb{R}$.
Subsequently, computations proceed in a bottom-up recursive pattern until
reaching the tree’s root in the top level. A representation of the computation
flow of such tree $\mathcal{T}$ is presented in Fig. 3.

####  3.4.2 Implementation of FEX

The FEX solution represented as $u(\mathbf{x};\mathcal{T},\bm{e},\bm{\theta})$
is derived by tackling the CO problem as:

$\displaystyle\min_{\bm{e},\bm{\theta}}\mathcal{L}(u(\cdot;\mathcal{T},\bm{e},\bm{\theta})).$
(4)

To achieve this, FEX adopts a two-phase approach:

Operator sequence optimization: The primary goal here is to refine the
operator sequence $\bm{e}$ such that it reflects the structure of the
authentic dynamics.

Parameter optimization: Following the sequence refinement, FEX then optimizes
the parameter set $\bm{\theta}$ aiming to minimize the functional described in
(4). The FEX framework is organized into four components:

  * •

Score Computation: To discover the solution’s structure, FEX deploys a mix-
order optimization algorithm. This step assesses the score of the operator
sequence $\bm{e}$. RBM is employed to accelerate the score’s calculation over
the network data.

  * •

Operator Sequence Generation: Here, a neural network (NN) is leveraged to
conceptualize a controller. This controller subsequently outputs a probability
mass function, which is used for sampling optimal operator sequences.

  * •

Controller Update: This step refines the parameters of the controller based on
the reward outcomes of the generated operator sequences, guiding it to yield
more effective sequences over training iterations.

  * •

Candidate Optimization: FEX dynamically keeps a pool of the operator sequences
with the highest rewards during the search process. Upon completion of
training, each sequence candidate undergoes a fine-tuning phase. The ultimate
objective is to identify the best operator sequence as the approximation of
system dynamics.

Subsequent sections will elaborate on the details of each of these components.

####  3.4.3 Score computation

The evaluation of an operator sequence $\bm{e}$ through its score is pivotal
during training. This score guides the controller’s parameter updates, aiming
to produce optimal probability mass functions for the sampling of efficient
operators to approximate the network dynamics. The score of $\bm{e}$, denoted
as $S(\bm{e})$, is defined as follows:

$\displaystyle S(\bm{e}):=\big{(}1+L(\bm{e})\big{)}^{-1},\quad{\rm where}\quad
L(\bm{e}):=\min\\{\mathcal{L}(u(\cdot;\mathcal{T},\bm{e},\bm{\theta}))|\bm{\theta}\\}.$
(5)

For an efficient computation of the score $S(\bm{e})$, we employ a hybrid mix-
order optimization strategy to update the parameter $\bm{\theta}$: We first
take $T_{1}$ iterations using a first-order optimization algorithm, such as
stochastic gradient descent [59] or Adam [41], and this yields
$\bm{\theta^{\bm{e}}}_{T_{1}}$. $\bm{\theta^{\bm{e}}}_{T_{1}}$ will be used as
a warm start, followed by a $T_{2}$ iterations of a second-order optimization
method, such as Newton’s method [4] or BFGS [20], resulting in
$\bm{\theta^{\bm{e}}}_{T_{1}+T_{2}}$.

Random Batch Method To address the computational burden when performing the
total $T_{1}+T_{2}$ optimization steps, we propose to approximate the
interaction terms using only a subset of the particles at each step. In
particular, we divide the $N$ particles into $n$ batches and only consider the
interactions within each batch, then the number of interactions that need to
be computed is significantly reduced. Assuming that each batch contains
approximately $p$ nodes ($n\approx\frac{N}{p}$), then the number of
interactions to be computed for each batch is on the order of $p^{2}$. Since
there are $n$ such batches, the total number of interactions at each step is
on the order of $np^{2}=\frac{N}{p}p^{2}=Np$, which is $O(N)$ since $p$ is a
constant independent of $N$. This significantly reduces the computational cost
when evaluating the loss function during the $T_{1}+T_{2}$ steps of
optimization.

Upon completion of these steps, the score associated with the operator
sequence $\bm{e}$ is denoted as:

$\displaystyle
S(\bm{e})=\big{(}1+\mathcal{L}(u(\cdot;\mathcal{T},\bm{e},\bm{\theta}_{T_{1}+T_{2}}^{\bm{e}}))\big{)}^{-1}.$
(6)

####  3.4.4 Operator Sequence Generation

The primary objective of the controller is to generate operator sequences that
achieve high scores throughout the training process. We denote the controller
as $\bm{\chi}_{\Phi}$, which is parameterized by a fully connected neural
network with parameters $\Phi$. Given an operator sequence $\bm{e}$ comprising
$s$ nodes, the controller $\bm{\chi}_{\Phi}$ produces probability mass
functions $\bm{p_{\Phi}}^{i}$ for $i=1,\ldots,s$. The operator $e_{j}$ is then
sampled based on the probability mass function $\bm{p_{\Phi}}^{j}$. To
encourage exploration within the set of operators, the $\epsilon$-greedy
strategy [66] is implemented. To be more specific, the operator $e_{i}$ is
selected from a uniform distribution over the entire operator set with a
probability of $\epsilon<1$.

####  3.4.5 Controller Update

As previously discussed, we introduce a controller represented by the neural
network $\bm{\chi}_{\Phi}$, to generate a probability mass function for each
node in the binary tree. This transformation changes the problem from a CO
problem into a continuous optimization one. The primary goal now becomes the
optimization of the controller parameters $\Phi$ so that the controller
$\bm{\chi}_{\Phi}$ is capable of producing optimal probability mass function
for each node. To achieve this, we consider the expected value of scores of
operator sequences sampled from the controller $\bm{\chi}_{\Phi}$, i.e.,

$\displaystyle\mathcal{J}(\Phi):=\mathbb{E}_{\bm{e}\sim\bm{\chi}_{\Phi}}S(\bm{e}).$
(7)

The derivative of (7) with respect to $\Phi$ is

$\displaystyle\nabla_{\Phi}\mathcal{J}(\Phi)=\mathbb{E}_{\bm{e}\sim\bm{\chi}_{\Phi}}\left\\{S(\bm{e})\sum_{i=1}^{s}\nabla_{\Phi}\log\left(\bm{p}_{\Phi}^{i}(e_{i})\right)\right\\},$
(8)

where $\bm{p}_{\Phi}^{i}(e_{i})$ is the probability of the sampled $e_{i}$. To
efficiently approximate the expectation (8), we sample a batch of operator
sequences $\\{\bm{e}^{(1)},\bm{e}^{(2)},\cdots,\bm{e}^{(M)}\\}$ each time and
compute

$\displaystyle\nabla_{\Phi}\mathcal{J}(\Phi)\approx\frac{1}{M}\sum_{k=1}^{M}\left\\{S(\bm{e}^{(k)})\sum_{i=1}^{s}\nabla_{\Phi}\log\left(\bm{p}_{\Phi}^{i}(e_{i}^{(k)})\right)\right\\},$
(9)

where $M$ is the batch size. Subsequently, the model parameter $\Phi$ is
updated using gradient ascent, denoted as
$\Phi\leftarrow\Phi+\eta\nabla_{\Phi}\mathcal{J}(\Phi)$. Nonetheless, the
practical aim is to identify the operator sequence $\bm{e}$ with the highest
score, rather than optimizing the average scores of all generated operator
sequences. Therefore, we follow the approach in [56] and consider

$\displaystyle\mathcal{J}(\Phi)=\mathbb{E}_{\bm{e}\sim\bm{\chi}_{\Phi}}\left\\{S(\bm{e})|S(\bm{e})\geq
S_{\nu,\Phi}\right\\},$ (10)

where $S_{\nu,\Phi}$ represents the $(1-\nu)\times 100\%$-quantile of the
score distribution generated by $\bm{\chi}_{\Phi}$. In a discrete form, the
gradient computation becomes

$\displaystyle\nabla_{\Phi}\mathcal{J}(\Phi)\approx\frac{1}{M}\sum_{k=1}^{M}\left\\{(S(\bm{e}^{(k)})-\hat{S}_{\nu,\Phi})\mathbf{1}_{\\{S(\bm{e}^{(k)})\geq\hat{S}_{\nu,\Phi}\\}}\sum_{i=1}^{s}\nabla_{\Phi}\log\left(\bm{p}_{\Phi}^{i}(e_{i}^{(k)})\right)\right\\},$
(11)

where $\mathbf{1}$ denotes an indicator function which evaluates to $1$ when
the condition is met, and $0$ otherwise. Meanwhile, $\hat{S}_{\nu,\Phi}$
represents the $(1-\nu)$-quantile of the scores within the set
$\\{S(\mathbf{e}^{(i)})\\}_{i=1}^{N}$.

####  3.4.6 Candidate Optimization

During the training phase of optimizing the controller parameters $\Phi$,
efficient evaluation of expressions sampled from the controller is crucial. As
introduced in Section 3.4.3, this is achieved by computing the score
$S(\mathbf{e})$ through $T_{1}+T_{2}$ “coarse-tune” iterations. However, this
coarse score might not reflect the performance of operator sequence $\bm{e}$
owing to the optimization of a nonconvex function. Therefore, we maintain a
pool $\mathbb{P}$ with a fixed size of $K$ during the training of controller,
which dynamically holds the top-performing candidate expressions. When the
controller’s training is completed, we further refine the coefficients of each
candidate expression $\bm{e}$ within the pool $\mathbb{P}$ to attain high
accuracy. To enhance the robustness of our method against link perturbations,
we employ a random sampling strategy. In particular, for each candidate
function, we perform random sampling of $L$ times, where we randomly sample
$S$ nodes from the graph each time. For each batch, we perform $T_{3}$
iterations of optimizations based on the objective function
$\mathcal{L}(u(\cdot;\mathcal{T},\bm{e},\bm{\theta}))$ with a first-order
algorithm with a small learning rate, and the final coefficients
$\hat{\bm{\theta}}$ of the expression are obtained by averaging the
coefficients over $L$ times.

##  4 Dynamics learning results

In this section, we present the results of identification of HR, FHN, and
coupled Rössler dynamics with FEX on the Scale-Free (SF) network with clean
data. The set up of SF network follows the description in Appendix 6.1. We
choose SF network since it demonstrates more complex phenomena than the ER
model, and many real-world networks have scale-free property. During the
offline stage, we generate the time series data with the true dynamics (ODEs).
During the inference stage, we use FEX only with the time series data to
identify the governing equations of the system, without accessing to the
genuine ODEs. Throughout all numerical experiments, we set the batch size in
RBM as 32. In addition, we provide the depths of FEX trees used in each of the
considered dynamics in the numerical examples in Table 1.

Dynamics | Depth of Tree
---|---
HR | 4
FHN | 3
Coupled Rössler | 3
Table 1: Tree depth of FEX tree for each dynamics.

###  4.1 Numerical derivatives and error metric

To approximate the time-varying derivative of each node, we apply the five-
point approximation [60]

$\dot{x}_{t}\approx\frac{x_{t-2\delta t}-8x_{t-\delta t}+8x_{t+\delta
t}-x_{t+2\delta t}}{12\delta t},$ (12)

where $\delta t$ is the time step.

To evaluate the accuracy of inferred dynamics, we use the metric of symmetric
mean absolute percentage error (sMAPE) [21]:

$\mathrm{sMAPE}=\frac{1}{m}\sum_{i=1}^{m}\frac{\left|D_{i}-R_{i}\right|}{\left(\left|D_{i}\right|+\left|R_{i}\right|\right)},$
(13)

where $m$ is the total number of terms of both true and inferred functions,
$D_{i}$ is the inferred coefficient and $R_{i}$ is the true coefficient. sMAPE
is an appropriate metric in evaluating the performance of baselines since it
not only captures the coefficient discrepancies but also penalizes incorrect
inferred terms. sMAPE ranges from 0 to 1, and the lower the value the higher
the accuracy, and vice versa.

###  4.2 Hindmarsh-Rose dynamics

The HR model [34, 64, 70] was introduced by Hindmarsh and Rose in the early
1980s to describe neuronal activity. The HR model is a simplified model of
Hodgkin-Huxley model [29] while still reproducing many of the dynamic
behaviors observed in real neurons. The governing equations of HR model is
given as

$\left\\{\begin{array}[]{l}\frac{dx_{i,1}}{dt}=x_{i,2}-ax_{i,1}^{3}+bx_{i,1}^{2}-x_{i,3}+I_{ext}+\epsilon\left(V_{\mathrm{syn}}-x_{i,1}\right)\sum_{j=1}^{N}\mathbf{A}_{ij}\sigma\left(x_{j,1}\right),\\\
\frac{dx_{i,2}}{dt}=c-ux_{i,1}^{2}-x_{i,2},\\\
\frac{dx_{i,3}}{dt}=r\left[s\left(x_{i,1}-x_{0}\right)-x_{i,3}\right],\end{array}\right.$
(14)

where the coupling term is

$\sigma\left(x_{j,1}\right)=\frac{1}{1+e^{-x_{j,1}}}.$

Here $x_{i,1}$ is the membrane potential of neuron $i$, $x_{i,2}$ is the
transport rate of ions across the membrane through the ion channels, and
$x_{i,3}$ is the adaptation current. Parameters are set as:
$a=1,b=3,c=1,u=5,s=4,r=0.004,x_{0}=-1.6$, $\epsilon=0.15,V_{syn}=2$, and
$I_{ext}=3.24$ is external current. Notice that the dynamic behavior of the HR
neuron can be adjusted by varying these parameters. By doing so, one can
observe various neuronal activities, such as tonic spiking, phasic spiking and
chaotic bursting. We generate the time series data with Equation (14), with
the adjacency matrix of SF network, terminal time $T=500$ and time step
$\delta t=0.01$.

Applying FEX to the neuronal activities data generated by HR dynamics on a
directed SF network, we obtain the inferred equations as

$\left\\{\begin{aligned}
\frac{d\hat{x}_{i,1}}{dt}&=0.9961x_{i,2}-1.0014x_{i,1}^{3}+2.9758x_{i,1}^{2}-0.9917x_{i,3}+3.2392\\\
&+\sum_{j=1}^{N}\mathbf{A}_{ij}\left(0.2985\sigma\left(x_{j,1}\right)-0.1477x_{i,1}\sigma\left(x_{j,1}\right)\right),\\\
\frac{d\hat{x}_{i,2}}{dt}&=1.0004-5.0001x_{i,1}^{2}-1.0001x_{i,2},\\\
\frac{d\hat{x}_{i,3}}{dt}&=0.0255+0.0182x_{i,1}-0.0050x_{i,3},\end{aligned}\right.$
(15)

From Eq. (15), we observe that FEX not only identifies all the correct terms
in the dynamics in each dimension of the node, but also learns the coefficient
of each term with relative errors lower than $1\%$. The neuronal activities
generated by FEX and true governing equations are shown in Fig. 4.

Figure 4: True HR dynamics and the dynamics identified by FEX.

###  4.3 FitzHugh-Nagumo dynamics

We further test our method to the neuronal activities data generated by FHN
dynamics [19]. The FHN model is another mathematical representation used to
describe the spiking behavior of neurons. It is defined by a system of two
ordinary differential equations, capturing the main features of excitability
in nerve membrane dynamics. The equations governing the FHN neuronal network
dynamics are

$\left\\{\begin{array}[]{l}\frac{dx_{i,1}}{dt}=x_{i,1}-x_{i,1}^{3}-x_{i,2}-\epsilon\sum_{j=1}^{N}\mathbf{A}_{ij}\frac{\left(x_{j,1}-x_{i,1}\right)}{k_{i}^{in}},\\\
\frac{dx_{i,2}}{dt}=a+bx_{i,1}+cx_{i,2},\end{array}\right.$ (16)

where the first component $x_{i,1}$ represents the membrane potential
containing self and interaction dynamics, $k_{i}^{in}$ is the in-degree of
neuron $i$ (representing the number of incoming connections to node $i$), and
$\epsilon=1$. The second component $x_{i,2}$ represents a recovery variable
where $a=0.28,b=0.5$ and $c=-0.04$. We generate the time series data with
Equation (16), with the adjacency matrix of SF network, terminal time $T=300$
and time step $\delta t=0.01$.

The equations inferred by our approach from the neuronal activities data
generated on a directed SF network are shown below:

$\left\\{\begin{array}[]{l}\frac{d\hat{x}_{i,1}}{dt}=0.9942x_{i,1}-0.9998x_{i,1}^{3}-0.9999x_{i,2}-1.0022\sum_{j=1}^{N}A_{ij}\frac{\left(x_{j,1}-x_{i,1}\right)}{k_{i}^{\text{in
}}}\\\
\frac{d\hat{x}_{i,2}}{dt}=0.2801+0.5000x_{i,1}-0.0400x_{i,2}.\end{array}\right.$

The trajectories generated by physical laws inferred by FEX and the true
governing equations are shown in Fig. 5. FEX has shown to accurately infer
both dimensions of FHN system and coincide closely with the true dynamics.

Figure 5: True FHN dynamics and the dynamics identified by FEX.

###  4.4 Coupled Rössler dynamics

We also consider the coupled Rössler oscillators [58, 67], which is a
classical model often used to study chaotic dynamics and synchronization
phenomena in complex networks. We generate chaotic activities data according
to the true equations governing heterogeneous Rössler dynamics [7, 61],

$\left\\{\begin{array}[]{l}\frac{dx_{i,1}}{dt}=-\omega_{i}x_{i,2}-x_{i,3}+\epsilon\sum_{j=1}^{n}A_{ij}\left(x_{j,1}-x_{i,1}\right),\\\
\frac{dx_{i,2}}{dt}=\omega_{i}x_{i,1}+ax_{i,2},\\\
\frac{dx_{i,3}}{dt}=b+x_{i,3}\left(x_{i,1}+c\right),\end{array}\right.$ (17)

where coupling strength $\epsilon=0.15$, $a=0.2,b=0.2,$ and $c=-5.7$. The
natural frequencies of the oscillators, denoted by $\omega_{i}$, follow a
normal distribution with a mean value of 1 and a standard deviation of 0.1. We
generate the time series data with Equation (17), with the adjacency matrix of
SF network, terminal time $T=100$ and time step $\delta t=0.01$.

The equations inferred by FEX from the data generated on a directed SF network
with $\epsilon=0.15$ are

$\left\\{\begin{array}[]{l}\frac{d\hat{x}_{i,1}}{dt}=-1.0093x_{i,2}-1.0027x_{i,3}+0.1491\sum_{j=1}^{n}A_{ij}\left(x_{j,1}-x_{i,1}\right),\\\
\frac{d\hat{x}_{i,2}}{dt}=0.9909x_{i,1}+0.2030x_{i,2},\\\
\frac{d\hat{x}_{i,3}}{dt}=0.1967+0.9987x_{i,3}x_{i,1}-5.6653x_{i,3},\end{array}\right.$
(18)

Note that the coefficient -1.0093 of $x_{i,2}$ in the first equation of Eq.
(18) and 0.9909 of $x_{i,1}$ in the second equation of Eq. (18) are estimated
average value of nodes’ natural frequencies, otherwise inferring a distinctive
frequency for each node would lead to an $n-$fold increase in the
dimensionality of model space. The trajectories generated by physical laws
inferred by FEX and the true governing equations are shown in Fig. 6.

Figure 6: True coupled Rössler dynamics and the dynamics identified by FEX.

###  4.5 Details of performance of random batch method

As introduced before, computation on a large complex network is demanding or
even infeasible. Therefore, we incorporate the RBM strategy into the FEX
algorithm to accelerate the identification of the optimal mathematical
structure of network dynamics while maintaining a high accuracy.

To demonstrate the efficiency of RBM in accelerating CO in the training
process of FEX, we create a ER model of network size $64,128,256$, and $512$,
respectively. Then, we use the time series data of FHN model on the four ER
models, and implement FEX to infer the dynamics. We record the clock time of
computation of $T_{1}+T_{2}$ coarse-tune in the FEX training algorithm by full
interaction method and RBM (with a batch size of 32), respectively. The log-
log plot of time vs. number of particles is shown in Fig. 7. We observe that
the fitting line clock time of the full interaction model on different network
sizes has a slope of approximately 2. The RBM’s slope is roughly 1, indicating
that the RBM demonstrates a linear scaling performance.

Figure 7: Clock time comparisons of full interaction model and RBM. The
$x$-axis shows 64, 128, 256, and 512 particles of the ER network,
respectively. The $y$-axis shows the clock time of one iteration computation
for each network size.

Furthermore, we comment that we only use RBM strategy in the CO part of FEX to
identify the optimal structure of the dynamics. The CO stage is the
computational bottleneck of FEX, and RBM has demonstrated its capability to
determine the correct structure of dynamics while accelerating the computation
significantly. After the structures of top-performing candidates are learned,
we use the full interaction model in the $T_{3}$ iterations of coefficients
fine-tuning stage to achieve high accuracy of dynamics inferred by FEX. The
computational burden at this stage is manageable since we only have to fine-
tune the coefficients of the top $K$ candidates in the pool $\mathbb{P}$.

##  5 Conclusion

In this work, we proposed FEX as a novel methodology to learn the dynamics on
complex networks. FEX has shown its capability to accurately discover physical
laws on various synthetic networks, even with noise and low resolution of
data. Furthermore, we incorporate RBM into our proposed FEX algorithm, which
scales up FEX to handle large data when the network size increases
significantly. Therefore, we have demonstrated that FEX has the potential to
identify real-world dynamics in an interpretable and reliable way.

There are also questions that remained to be addressed. First, throughout our
work, we consider only the pairwise interaction dynamics in the physical laws.
However, it is common for complex systems to exhibit three-point or higher-
order couplings among nodes in the network. Yet, it is straightforward to
adapt such a setting into the FEX methodology, with the extra computational
cost of including more binary trees to represent higher-order dynamics.
Second, dynamics often demonstrate implicit stochasticity in real-world data,
which can be described by stochastic differential equations (SDEs) [16, 26].
Lastly, this work focuses on classical dynamics with an explicit form of
governing equations. However, dynamics are largely unknown or underexplored in
most situations. Hence, it is interesting to observe how accurately FEX can
approximate empirical dynamics and the insights it provides from inferring
complex network dynamics.

## References

  * [1] Réka Albert and Albert-László Barabási.  Statistical mechanics of complex networks.  Reviews of modern physics, 74(1):47, 2002.
  * [2] Shun-ichi Amari.  Backpropagation and stochastic gradient descent method.  Neurocomputing, 5(4-5):185–196, 1993.
  * [3] Hassan Arbabi, Judith E Bunder, Giovanni Samaey, Anthony J Roberts, and Ioannis G Kevrekidis.  Linking machine learning with multiscale numerics: Data-driven discovery of homogenized equations.  Jom, 72:4444–4457, 2020.
  * [4] Mordecai Avriel.  Nonlinear programming: analysis and methods.  Courier Corporation, 2003.
  * [5] Albert-László Barabási.  Scale-free networks: a decade and beyond.  science, 325(5939):412–413, 2009.
  * [6] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  science, 286(5439):509–512, 1999.
  * [7] Baruch Barzel, Yang-Yu Liu, and Albert-László Barabási.  Constructing minimal models for complex system dynamics.  Nature communications, 6(1):7186, 2015.
  * [8] Tom Bertalan, Felix Dietrich, Igor Mezić, and Ioannis G Kevrekidis.  On learning hamiltonian systems from data.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(12), 2019.
  * [9] Léon Bottou.  Large-scale machine learning with stochastic gradient descent.  In Proceedings of COMPSTAT’2010: 19th International Conference on Computational StatisticsParis France, August 22-27, 2010 Keynote, Invited and Contributed Papers, pages 177–186. Springer, 2010.
  * [10] Cameron P Bracken, Hamish S Scott, and Gregory J Goodall.  A network-biology perspective of microrna function and dysfunction in cancer.  Nature Reviews Genetics, 17(12):719–732, 2016.
  * [11] Steven L Brunton, Joshua L Proctor, and J Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the national academy of sciences, 113(15):3932–3937, 2016.
  * [12] Jose Casadiego, Mor Nitzan, Sarah Hallerberg, and Marc Timme.  Model-free inference of direct network interactions from nonlinear collective dynamics.  Nature communications, 8(1):2192, 2017.
  * [13] Dong-Yeon Cho, Yoo-Ah Kim, and Teresa M Przytycka.  Chapter 5: Network biology approach to complex diseases.  PLoS computational biology, 8(12):e1002820, 2012.
  * [14] Bryan C Daniels and Ilya Nemenman.  Automated adaptive inference of phenomenological dynamical models.  Nature communications, 6(1):8133, 2015.
  * [15] Leon Danon, Ashley P Ford, Thomas House, Chris P Jewell, Matt J Keeling, Gareth O Roberts, Joshua V Ross, Matthew C Vernon, et al.  Networks and the epidemiology of infectious disease.  Interdisciplinary perspectives on infectious diseases, 2011, 2011.
  * [16] Gustavo Deco, Edmund T Rolls, and Ranulfo Romo.  Stochastic dynamics as a principle of brain function.  Progress in neurobiology, 88(1):1–16, 2009.
  * [17] Felix Dietrich, Alexei Makeev, George Kevrekidis, Nikolaos Evangelou, Tom Bertalan, Sebastian Reich, and Ioannis G Kevrekidis.  Learning effective stochastic differential equations from microscopic simulations: Linking stochastic numerics to deep learning.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(2), 2023.
  * [18] Petros Drineas, Ravi Kannan, and Michael W Mahoney.  Fast monte carlo algorithms for matrices ii: Computing a low-rank approximation to a matrix.  SIAM Journal on computing, 36(1):158–183, 2006.
  * [19] Richard FitzHugh.  Impulses and physiological states in theoretical models of nerve membrane.  Biophysical journal, 1(6):445–466, 1961.
  * [20] Roger Fletcher.  Practical methods of optimization.  John Wiley & Sons, 2013.
  * [21] Benito E Flores.  A pragmatic view of accuracy measurement in forecasting.  Omega, 14(2):93–98, 1986.
  * [22] Cornelius Fritz, Emilio Dorigatti, and David Rügamer.  Combining graph neural networks and spatio-temporal disease models to improve the prediction of weekly covid-19 cases in germany.  Scientific Reports, 12(1):3930, 2022.
  * [23] Kai Fukami, Takaaki Murata, Kai Zhang, and Koji Fukagata.  Sparse identification of nonlinear dynamics with low-dimensionalized flow representations.  Journal of Fluid Mechanics, 926:A10, 2021.
  * [24] Junyi Gao, Rakshith Sharma, Cheng Qian, Lucas M Glass, Jeffrey Spaeder, Justin Romberg, Jimeng Sun, and Cao Xiao.  Stan: spatio-temporal attention network for pandemic prediction using real-world evidence.  Journal of the American Medical Informatics Association, 28(4):733–743, 2021.
  * [25] Ting-Ting Gao and Gang Yan.  Autonomous inference of complex network dynamics from incomplete and noisy data.  Nature Computational Science, 2(3):160–168, 2022.
  * [26] Mikhail Genkin, Owen Hughes, and Tatiana A Engel.  Learning non-stationary langevin dynamics from stochastic observations of latent trajectories.  Nature communications, 12(1):5986, 2021.
  * [27] Jesús Gómez-Gardeñes and Yamir Moreno.  From scale-free to erdos-rényi networks.  Physical Review E, 73(5):056124, 2006.
  * [28] Raul González-García, Ramiro Rico-Martìnez, and Ioannis G Kevrekidis.  Identification of distributed parameter systems: A neural net based approach.  Computers & chemical engineering, 22:S965–S968, 1998.
  * [29] John Guckenheimer and Ricardo A Oliva.  Chaos in the hodgkin–huxley model.  SIAM Journal on Applied Dynamical Systems, 1(1):105–114, 2002.
  * [30] Aric Hagberg, Pieter Swart, and Daniel S Chult.  Exploring network structure, dynamics, and function using networkx.  Technical report, Los Alamos National Lab.(LANL), Los Alamos, NM (United States), 2008.
  * [31] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM review, 53(2):217–288, 2011.
  * [32] Moritz Hardt, Ben Recht, and Yoram Singer.  Train faster, generalize better: Stability of stochastic gradient descent.  In International conference on machine learning, pages 1225–1234. PMLR, 2016.
  * [33] Tong He, Zhi Zhang, Hang Zhang, Zhongyue Zhang, Junyuan Xie, and Mu Li.  Bag of tricks for image classification with convolutional neural networks.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 558–567, 2019.
  * [34] James L Hindmarsh and RM Rose.  A model of the nerve impulse using two first-order differential equations.  Nature, 296(5853):162–164, 1982.
  * [35] Zhongyi Jiang, Chunmei Wang, and Haizhao Yang.  Finite expression methods for discovering physical laws from data.  arXiv preprint arXiv:2305.08342, 2023.
  * [36] Shi Jin, Lei Li, and Jian-Guo Liu.  Random batch methods (rbm) for interacting particle systems.  Journal of Computational Physics, 400:108877, 2020.
  * [37] Eurika Kaiser, J Nathan Kutz, and Steven L Brunton.  Sparse identification of nonlinear dynamics for model predictive control in the low-data limit.  Proceedings of the Royal Society A, 474(2219):20180335, 2018.
  * [38] Mihyun Kang and Zdenek Petrasek.  Random graphs: Theory and applications from nature to society to the brain.  Internationale Mathematische Nachrichten, 227:1–24, 2014.
  * [39] George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.  Physics-informed machine learning.  Nature Reviews Physics, 3(6):422–440, 2021.
  * [40] Felix P Kemeth, Sergio Alonso, Blas Echebarria, Ted Moldenhawer, Carsten Beta, and Ioannis G Kevrekidis.  Black and gray box learning of amplitude equations: Application to phase field systems.  Physical Review E, 107(2):025305, 2023.
  * [41] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [42] Senwei Liang and Haizhao Yang.  Finite expression method for solving high-dimensional partial differential equations.  arXiv preprint arXiv:2206.10121, 2022.
  * [43] Gipsi Lima-Mendez and Jacques Van Helden.  The powerful law of the power law and other myths in network biology.  Molecular BioSystems, 5(12):1482–1493, 2009.
  * [44] Jingyi Lin and Yifang Ban.  Complex network topology of transportation systems.  Transport reviews, 33(6):658–685, 2013.
  * [45] Bing Liu, Wei Luo, Gang Li, Jing Huang, and Bo Yang.  Do we need an encoder-decoder to model dynamical systems on networks?  arXiv preprint arXiv:2305.12185, 2023.
  * [46] Ji Liu and Stephen J Wright.  Asynchronous stochastic coordinate descent: Parallelism and convergence properties.  SIAM Journal on Optimization, 25(1):351–376, 2015.
  * [47] Ji Liu, Steve Wright, Christopher Ré, Victor Bittorf, and Srikrishna Sridhar.  An asynchronous parallel stochastic coordinate descent algorithm.  In International Conference on Machine Learning, pages 469–477. PMLR, 2014.
  * [48] Yuxuan Liu, Zecheng Zhang, and Hayden Schaeffer.  Prose: Predicting operators and symbolic expressions using multimodal transformers.  arXiv preprint arXiv:2309.16816, 2023.
  * [49] Alun L Lloyd and Steve Valeika.  Network models in epidemiology: an overview.  Complex population dynamics: nonlinear modeling in ecology, epidemiology and genetics, pages 189–214, 2007.
  * [50] Zichao Long, Yiping Lu, and Bin Dong.  Pde-net 2.0: Learning pdes from data with a numeric-symbolic hybrid deep network.  Journal of Computational Physics, 399:108925, 2019.
  * [51] Bethany Lusch, J Nathan Kutz, and Steven L Brunton.  Deep learning for universal linear embeddings of nonlinear dynamics.  Nature communications, 9(1):4950, 2018.
  * [52] Yi-An Ma, Tianqi Chen, and Emily Fox.  A complete recipe for stochastic gradient mcmc.  Advances in neural information processing systems, 28, 2015.
  * [53] Per-Gunnar Martinsson, Vladimir Rokhlin, and Mark Tygert.  A randomized algorithm for the decomposition of matrices.  Applied and Computational Harmonic Analysis, 30(1):47–68, 2011.
  * [54] Charles Murphy, Edward Laurence, and Antoine Allard.  Deep learning of contagion dynamics on complex networks.  Nature Communications, 12(1):4720, 2021.
  * [55] Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani.  Epidemic processes in complex networks.  Reviews of modern physics, 87(3):925, 2015.
  * [56] Brenden K Petersen, Mikel Landajuela Larma, Terrell N. Mundhenk, Claudio Prata Santiago, Soo Kyung Kim, and Joanne Taery Kim.  Deep symbolic regression: Recovering mathematical expressions from data via risk-seeking policy gradients.  In International Conference on Learning Representations, 2021.
  * [57] Yorgos M Psarellis, Seungjoon Lee, Tapomoy Bhattacharjee, Sujit S Datta, Juan M Bello-Rivas, and Ioannis G Kevrekidis.  Data-driven discovery of chemotactic migration of bacteria via machine learning.  arXiv preprint arXiv:2208.11853, 2022.
  * [58] Michael G Rosenblum, Arkady S Pikovsky, and Jürgen Kurths.  Phase synchronization of chaotic oscillators.  Physical review letters, 76(11):1804, 1996.
  * [59] David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  nature, 323(6088):533–536, 1986.
  * [60] Timothy Sauer.  Numerical solution of stochastic differential equations in finance.  In Handbook of computational finance, pages 529–550. Springer, 2011.
  * [61] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  science, 324(5923):81–85, 2009.
  * [62] Harold Soh, Sonja Lim, Tianyou Zhang, Xiuju Fu, Gary Kee Khoon Lee, Terence Gih Guang Hung, Pan Di, Silvester Prakasam, and Limsoon Wong.  Weighted complex network analysis of travel routes on the singapore public transportation system.  Physica A: Statistical Mechanics and its Applications, 389(24):5852–5863, 2010.
  * [63] Zezheng Song, Maria K Cameron, and Haizhao Yang.  A finite expression method for solving high-dimensional committor problems.  arXiv preprint arXiv:2306.12268, 2023.
  * [64] Marco Storace, Daniele Linaro, and Enno de Lange.  The hindmarsh–rose neuron model: bifurcation analysis and piecewise-linear approximations.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 18(3), 2008.
  * [65] Yifan Sun, Linan Zhang, and Hayden Schaeffer.  Neupde: Neural network based ordinary and partial differential equations for modeling time-dependent data.  In Mathematical and Scientific Machine Learning, pages 352–372. PMLR, 2020.
  * [66] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [67] Longkun Tang, Xiaoqun Wu, Jinhu Lü, Jun-an Lu, and Raissa M D’Souza.  Master stability functions for complete, intralayer, and interlayer synchronization in multiplex networks of coupled rössler oscillators.  Physical Review E, 99(1):012304, 2019.
  * [68] Remco Van Der Hofstad.  Random graphs and complex networks.  Available on http://www. win. tue. nl/rhofstad/NotesRGCN. pdf, 11:60, 2009.
  * [69] Ivan Voitalov, Pim Van Der Hoorn, Remco Van Der Hofstad, and Dmitri Krioukov.  Scale-free networks well done.  Physical Review Research, 1(3):033034, 2019.
  * [70] X-J Wang.  Genesis of bursting oscillations in the hindmarsh-rose model and homoclinicity to a chaotic saddle.  Physica D: Nonlinear Phenomena, 62(1-4):263–274, 1993.
  * [71] Max Welling and Yee W Teh.  Bayesian learning via stochastic gradient langevin dynamics.  In Proceedings of the 28th international conference on machine learning (ICML-11), pages 681–688, 2011.
  * [72] Zengwang Xu and Robert Harriss.  Exploring the structure of the us intercity passenger air transportation network: a weighted complex network approach.  GeoJournal, 73:87–102, 2008.
  * [73] Chengxi Zang and Fei Wang.  Neural dynamics on complex networks.  In Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining, pages 892–902, 2020.
  * [74] Massimiliano Zanin, Xiaoqian Sun, and Sebastian Wandelt.  Studying the topology of transportation systems through complex networks: handle with care.  Journal of Advanced Transportation, 2018, 2018.
  * [75] Linan Zhang and Hayden Schaeffer.  On the convergence of the sindy algorithm.  Multiscale Modeling & Simulation, 17(3):948–972, 2019.

Acknowledgments C. W. was partially supported by National Science Foundation
under award DMS-2206332. H. Y. was partially supported by the US National
Science Foundation under awards DMS-2244988, DMS-2206333, and the Office of
Naval Research Award N00014-23-1-2007.

Competing interests The authors declare no competing interests.

##  6 Appendix

###  6.1 Networks

In our work, we consider several synthetic networks to demonstrate the
capabilities of FEX to identify dynamics on various network topologies.

Erdős-Rényi (ER) model: ER model [68, 38] is one of the foundational models
for random graphs. It is used widely to understand the properties of complex
networks in various fields, such as sociology, biology, and technology. In a
ER model, a graph is constructed by connecting nodes randomly with a
probability of $p$. We generate ER networks by using the
erdos_renyi_graph(n,p) in networkX [30], where $n=100,p=0.05$. Here, $n$
denotes the number of nodes in the graph and $p$ represents the probability to
connect two arbitrary nodes, respectively. This set up results a graph with an
average degree of $\langle k\rangle=4.95$. To create a directed graph based on
ER model, we allocate the direction of the link, i.e., from $i$ to $j$, from
$j$ to $i$, or a reciprocal connection between them (both $i$ to $j$ and $j$
to $i$ directions exist), with an equal likelihood (each of probability
$\frac{1}{3}$).

Scale-Free (SF) network: SF networks are ubiquitous in real-world applications
[69]. This type of network demonstrates an interesting phenomenon: nodes with
more early connections tend to attract even more connections, leading to the
“rich get richer” scenario. Therefore, the network’s degree follows a power
law distribution, implying a few nodes with very high degrees (hubs) and many
with low degrees. Networks structured according to the scale-free model are
constructed utilizing the Barabási-Albert (BA) [1, 6, 5] framework. Initially,
we build such network by building a complete graph of a small number $m_{0}$
of nodes. Then at each time step, we add a new node with $m\leq m_{0}$ edges.
The probability that the new node will connect to node $i$ is
$P(i)=\frac{k_{i}}{\sum_{j}k_{j}}$, where $k_{i}$ is the degree of node $i$,
and the sum is over all existing nodes $j$. Therefore, when the new node is
added, it has a higher probability of connecting to nodes with high degrees,
known as the “preferential attachment”. We use the barabasi_albert_graph
function (with parameters $n=100$, $m=5$) in NetworkX to generate SF network,
where $n$ is the number of nodes in the graph, and $m$ denotes each node
connects to $m$ existing nodes. Initially, every link is established as
bidirectional. Subsequently, a select portion of these unidirectional links
are randomly removed, albeit with the stipulation of maintaining a weakly
connected network architecture. The result is a directed SF network with an
average total degree of $\langle k\rangle=5.1$. An illustration of the ER
model and the SF network is provided in Fig. 8.

Figure 8: Visualization of ER model and SF (BA) network. (a) In ER network,
the graph is essentially random without any structures or patterns. (b) In SF
(BA) network, nodes closer to the boundary have generally low degrees, while
there are some interior nodes with a high number of connections (hubs).

Despite being synthetic, the ER model and the SF network are classic and
essential network models within the scientific community. Given that numerous
real-world networks exhibit scale-free properties, the ability of FEX to
identify physical laws on such networks accurately holds significant potential
for applications in real-world data.

###  6.2 Algorithm and experimental details

In this section, we provide an overview of the implementation details of FEX
of identification of network dynamics. The training algorithm of FEX is
summarized in Algorithm 1. The details of hyperparameters in the algorithm are
introduced as follows.

(1) Score computation. The score of an operator sequence $\bm{e}$ is updated
with the Adam optimizer with a learning rate of 0.001 across $T_{1}$
iterations. Following this initial phase, refinement of the score is carried
out through the BFGS optimizer with learning rate of 1, for up to 20
iterations.

(2) Operator sequence generation. In FEX framework, a binary tree of variable
depths is utilized, where the binary operator set is chosen as
$\mathbbm{B}=\\{+,-,\times\\}$, and the unary operator set as
$\mathbbm{U}=\\{0,1,\text{Id},(\cdot)^{2},(\cdot)^{3},(\cdot)^{4},\exp,\sin,\cos,\tanh,\text{sigmoid}\\}$.
The controller $\bm{\chi}_{\Phi}$ is parameterized as a neural network with
parameter $\Phi$. It takes a constant input of dimension 20, while the output
size is determined by the expression $n_{1}|\mathbbm{B}|+n_{2}|\mathbbm{U}|$.
Within this expression, $n_{1}$ and $n_{2}$ denote the counts of binary and
unary operators, respectively, and $|\cdot|$ denotes the cardinality of a set.

(3) Controller update. The update of policy gradient is implemented with a
batch size of $N=10$ over 300 iterations with the Adam optimizer at a learning
rate of 0.002. To encourage the exploration of new operators, the
$\epsilon$-greedy strategy is employed, wherein the probability of random
sampling for an $\bm{e}_{i}$ is fixed at 0.1.

(4) Candidate optimization. The capacity of the candidate pool is denoted by
$K=15$. For every $\bm{e}$ belonging to the set $\mathbb{P}$, the parameter
$\theta$ is optimized via the Adam optimizer, initiated at a learning rate of
0.001 across 20,000 iterations. The learning rate follows a cosine decay
schedule as in [33]. The number of sampling is chosen as $L=5$, and each
random sample contains $S=20$ nodes.

Algorithm 1 Fast FEX algorithm with RMB to identify dynamics on networks

Input: Adjacency matrix $\mathbf{A}$ of network; Time series data $\mathbf{x}$
of time steps $T_{s}$ of $N$ nodes in the network; Loss functional
$\mathcal{L}$; self-dynamics tree $\mathcal{T}_{1}$; Interaction dynamics tree
$\mathcal{T}_{2}$; Searching loop iteration $T$; Coarse-tune iteration $T_{1}$
with Adam; Coarse-tune iteration $T_{2}$ with BFGS; Fine-tune iteration
$T_{3}$ with Adam; Pool size $K$; Batch size $M$; Coefficient filtering
threshold $\tau$; Number of sampling $L$; Sampling batch size $S$;

Output: The solution
$u(\bm{x};\mathcal{T}_{1},\mathcal{T}_{2},\hat{\bm{e}},\hat{\mathbf{f}},\hat{\bm{\theta}}_{1},\hat{\bm{\theta}}_{2})$.

1:Initialize the agent $\bm{\chi}_{\phi,1}$ and $\bm{\chi}_{\phi,2}$ for tree
$\mathcal{T}_{1}$, and $\mathcal{T}_{2}$, respectively

2:$\mathbb{P}\leftarrow\\{\\}$,

3:for $t$ from $1$ to $T$ do

4: Sample $M$ operator sequences
$\\{\bm{e}^{(1)},\bm{e}^{(2)},\cdots,\bm{e}^{(M)}\\}$ from
$\bm{\chi}_{\phi,1}$

5: Sample $M$ operator sequences
$\\{\mathbf{f}^{(1)},\mathbf{f}^{(2)},\cdots,\mathbf{f}^{(M)}\\}$ from
$\bm{\chi}_{\phi,2}$

6: for $m$ from $1$ to $M$ do

7: Divide $\\{1,2,\ldots,N\\}$ into $n$ batches randomly

8: for each batch $\mathcal{C}_{q}$ do

9: Compute the loss of particle $i$ ($i\in\mathcal{C}_{q}$) by
$\frac{1}{T_{s}}\sum_{t=1}^{T_{s}}\left(\left\|\mathcal{T}_{1}(x_{i},t)+\sum_{j\in\mathcal{C}_{q}}{\mathbf{A}_{ij}}\mathcal{T}_{2}(x_{i},t)-\dot{\mathbf{x}_{i}}(t)\right\|^{2}\right)$

10: end for

11: Optimize
$\mathcal{L}(u(\bm{x};\mathcal{T}_{1},\mathcal{T}_{2},\bm{e}^{(m)},\mathbf{f}^{(m)},\bm{\theta}_{1},\bm{\theta}_{2}))$
by coarse-tune with $T_{1}+T_{2}$ iterations.

12: Compute the reward $R(\bm{e}^{(m)},\mathbf{f}^{(m)})$ of $\bm{e}^{(m)}$
and $\mathbf{f}^{(m)}$

13: if $\bm{e}^{(m)}$ and $\mathbf{f}^{(m)}$ belongs to the top-$K$ of
candidates then

14: $\mathbb{P}$.append($\bm{e}^{(m)}$ and $\mathbf{f}^{(m)}$)

15: $\mathbb{P}$ pops some $\bm{e}$ and $\mathbf{f}$ with the smallest reward
when overloading

16: end if

17: end for

18: Update $\bm{\chi}_{\phi,1}$ from Section 3.4.5

19: Update $\bm{\chi}_{\phi,2}$ from Section 3.4.5

20:end for

21:for $\bm{e}$ and $\mathbf{f}$ in $\mathbb{P}$ do

22: for $l$ from $1$ to $L$ do

23: Randomly sample $S$ nodes from the total $N$ nodes;

24: Fine-tune
$\mathcal{L}(u(\bm{x}^{S};\mathcal{T}_{1},\mathcal{T}_{2},\bm{e},\mathbf{f},\bm{\theta}_{1}^{S},\bm{\theta}_{2}^{S}))$
with $T_{3}$ iterations, apply coefficient filtering with threshold $\tau$.

25: end for

26: average $\bm{\theta}_{1}^{S},\bm{\theta}_{2}^{S}$ over $L$ runs to obtain
$\hat{\bm{\theta}}_{1},\hat{\bm{\theta}}_{2}$;

27:end for

28:return the expression with the smallest fine-tune error.

###  6.3 Robustness analyses

In practical situations, data are often of low resolution and noisy. To
evaluate the robustness of our proposed method, we perform a series of tests
under various conditions: low data resolution, spurious or missing connections
within the graph structure, and noisy data. To mimic a low-resolution data
scenario, we down-sample the time series of all nodes within the complex
network regularly. A random set of edges within the graph is either inserted
or removed to emulate situations with spurious or missing connections.
Gaussian noise is integrated into the time series data when simulating noise
corruption. This comprehensive testing framework aims to show that our method
demonstrates robust performance across a range of realistic and noisy data
environments. In the following subsections, we elaborate on each of the tests
respectively. We use Two-Phase and ARNI as our baselines for inferring
dynamics on complex network. Both of them require two comprehensive libraries,
$L_{F}$ (self-dynamics) and $L_{G}$ (interaction dynamics). We summarize them
in Table 2 and Table 3, respectively.

Functions | $L_{F}$
---|---
Polynomial | $x_{i},x_{i}^{2},x_{i}^{3},x_{i}^{4}$
Exponential | $\exp{(x_{i}})$
Activation |  $\sigma(x_{i})$, $\tanh(x_{i})$
Rescaling | $\frac{x_{i}}{k_{i}^{\text{in}}}$
Trigonometric |  $\sin(x_{i})$, $\cos{x_{i}}$
Table 2: Library functions (Two-Phase and ARNI) for self-dynamics. Functions | $L_{G}(x_{j})$ | $L_{G}(x_{i}x_{j})$ | $L_{G}(x_{j}-x_{i})$ | $x_{i}L_{G}(x_{j})$
---|---|---|---|---
Polynomial |  $x_{j}$, $x_{j}^{2}$ | $x_{i}x_{j},(x_{i}x_{j})^{2}$ |  $x_{j}-x_{i}$, $(x_{j}-x_{i})^{2}$ | $x_{i}x_{j}^{2}$
Exponential | $\exp(x_{j})$ | $\exp(x_{i}x_{j})$ | $\exp(x_{j}-x_{i})$ | $x_{i}\exp(x_{j})$
Activation | $\tanh(x_{j}),\sigma(x_{j})$ | $\tanh(x_{i}x_{j}),\sigma(x_{i}x_{j})$ | $\tanh(x_{j}-x_{i}),\sigma(x_{j}-x_{i})$ | $x_{i}\tanh(x_{j}),x_{i}\sigma(x_{j})$
Rescaling | $\frac{x_{j}}{k_{i}^{\text{in}}}$ | $\frac{x_{i}x_{j}}{k_{i}^{\text{in}}}$ | $\frac{x_{j}-x_{i}}{k_{i}^{\text{in}}}$ | $\frac{x_{i}x_{j}}{k_{i}^{\text{in}}}$
Trigonometric | $\sin(x_{j}),\cos(x_{j})$ | $\sin(x_{i}x_{j}),\cos(x_{i}x_{j})$ | $\sin(x_{j}-x_{i}),\cos(x_{j}-x_{i})$ | $x_{i}\sin(x_{j}),x_{i}\cos(x_{j})$
Table 3: Library functions (Two-Phase and ARNI) for interaction dynamics.

###  6.4 Low resolution by experimental techniques

In this subsection, we test the robustness of FEX with low resolution data. We
test this on the Coupled Rössler dynamics on the directed SF network. To be
more specific, we systematically perform down-sampling on the time series data
of each node within the network, aiming to evaluate the impact of reduced data
resolution on our methodology. Specifically, we apply two distinct down-
sampling rates of $10\%$ and $5\%$. This process allows us to emulate
scenarios of considerably lower data resolution, enabling an assessment of our
method’s performance and robustness under such constrained data conditions.

The dynamics of the first dimension of Coupled Rössler dynamics identified by
Two-Phase, ARNI and FEX are summarized in Table 4 and Table 5 with down-
sampling rates $10\%$ and $5\%$, respectively. It can be observed that FEX
consistently achieves better coefficient identification accuracies under
different sampling rates, especially in the coefficients of interaction
dynamics $\mathbf{G}(x_{i},x_{j})$. The sMAPEs of all compared methods with
different down-sampling rates are summarized in Table 6.

| True Dynamics | Two-Phase | ARNI | FEX
---|---|---|---|---
$\mathbf{F}(x_{i})$ | $-\omega_{i}x_{i,2}$ | $-1.0113x_{i,2}$ | $-1.0123x_{i,2}$ | $-1.0045x_{i,2}$
| $-x_{i,3}$ | $-0.9783x_{i,3}$ | $-0.9625x_{i,3}$ | $-0.9991x_{i,3}$
$\mathbf{G}(x_{i},x_{j})$ | $0.15(x_{j,1}-x_{i,1})$ | $0.0973\left(x_{j,1}-x_{i,1}\right)$ | $0.1263\left(x_{j,1}-x_{i,1}\right)$ | $0.1431\left(x_{j,1}-x_{i,1}\right)$
Table 4: Numerical results for coupled Rössler dynamics (down-sampling rate
$10\%$) in self-dynamics $\mathbf{F}(x_{i})$ and interaction dynamics
$\mathbf{G}(x_{i},x_{j})$ term by term.

| True Dynamics | Two-Phase | ARNI | FEX
---|---|---|---|---
$\mathbf{F}(x_{i})$ | $-\omega_{i}x_{i,2}$ | $-1.0114x_{i,2}$ | $-1.0275x_{i,2}$ | $-1.0063x_{i,2}$
| $-x_{i,3}$ | $-0.9251x_{i,3}$ | $-0.9338x_{i,3}$ | $-0.9741x_{i,3}$
$\mathbf{G}(x_{i},x_{j})$ | $0.15(x_{j,1}-x_{i,1})$ | $0.0826\left(x_{j,1}-x_{i,1}\right)$ | $0.1045\left(x_{j,1}-x_{i,1}\right)$ | $0.1348\left(x_{j,1}-x_{i,1}\right)$

Table 5: Numerical results for coupled Rössler dynamics (down-sampling rate
$5\%$) in self-dynamics $\mathbf{F}(x_{i})$ and interaction dynamics
$\mathbf{G}(x_{i},x_{j})$ term by term.

.

Down-sampling rate | Two-Phase | ARNI | FEX
---|---|---|---
$10\%$ | 0.0765 | 0.0370 | 0.0087
$5\%$ | 0.1114 | 0.0755 | 0.0232
Table 6: sMAPE of Coupled Rössler dynamics with low resolution data.

###  6.5 Spurious and missing links

In real-world scenarios, accurately observing the topology of networks is
often impractical, so a thorough evaluation of our method’s robustness in the
presence of spurious or missing links is necessary. Therefore, we select
directed SF networks as our testbed, wherein we introduce perturbations by
randomly adding or removing certain edges. Consequently, during offline phase,
we generate the time series data with the true adjacency matrix of the SF
network. During the online phase, we use perturbed network adjacency matrix to
infer the underlying physical laws of the system. In particular, we randomly
add or remove $15\%$ links in the underlying SF networks, respectively. The
dynamics identified by Two-Phase, ARNI and FEX are summarized in Table 7 and
Table 8, respectively. We note that FEX exhibits robustness against graph
perturbations, in scenarios of both spurious and missing links, attributed to
its design which aims to limit the size of the generated expressions and the
fine-tuning of parameters. The Two-Phase method is comparable to FEX, and
shows decent robustness against graph topology perturbations, owing to its
topological sampling in the second phase. On the contrary, ARNI falls short in
performance on noisy network structures as it was initially proposed to infer
network topology, thereby potentially leading to overfitting. ARNI infers many
redundant terms in the dynamics with highly inaccurate coefficients of each
term. The sMAPEs of all compared methods with $15\%$ spurious and missing
links are summarized in Table 9.

| True Dynamics | Two-Phase | ARNI | FEX
---|---|---|---|---
$\mathbf{F}(x_{i})$ | $x_{i,1}$ | $0.9897x_{i,1}$ | $6.8282x_{i,1}$ | $0.9807x_{i,1}$
| $-x_{i,1}^{3}$ | $-0.9733x_{i,1}^{3}$ | $-0.1202x_{i,1}^{3}$ | $-0.9979x_{i,1}^{3}$
| $-x_{i,2}$ | $-0.9861x_{i,2}$ | $-1.3184x_{i,2}$ | $-0.9975x_{i,2}$
|  |  | $1.3564x_{i,1}^{2}$ |
|  |  | -88.3806 |
|  |  | $-1.8841\exp(x_{i,1})$ |
$\mathbf{G}(x_{i},x_{j})$ | $-(x_{j,1}-x_{i,1})$ | $-0.9862\left(x_{j,1}-x_{i,1}\right)$ | $-22.6880(x_{j,1}-x_{i,1})$ | $-1.0151\left(x_{j,1}-x_{i,1}\right)$
|  |  | $-18.4197\sin(x_{j,1}-x_{i,1})$ |
|  |  | $181.6294\sigma(x_{j,1}-x_{i,1})$ |

Table 7: Numerical results for FHN dynamics (randomly add $15\%$ edges) in
self-dynamics $\mathbf{F}(x_{i})$ and interaction dynamics
$\mathbf{G}(x_{i},x_{j})$ term by term (each term of true or inferred
$\mathbf{G}(x_{i},x_{j})$ is understood to be divided by $k_{\text{in}}$,
i.e., $\frac{\mathbf{G}(x_{i},x_{j})}{k_{\text{in}}}$, where $k_{\text{in}}$
is the in-degree of node $i$).

.

| True Dynamics | Two-Phase | ARNI | FEX
---|---|---|---|---
$\mathbf{F}(x_{i})$ | $x_{i,1}$ | $0.9828x_{i,1}$ | $0.6644x_{i,1}$ | $0.9835x_{i,1}$
| $-x_{i,1}^{3}$ | $-0.9910x_{i,1}^{3}$ | $-0.5332x_{i,1}^{3}$ | $-1.0038x_{i,1}^{3}$
| $-x_{i,2}$ | $-0.9837x_{i,2}$ | $-1.2174x_{i,2}$ | $-0.9916x_{i,2}$
|  |  | $-0.2063\exp(x_{i,1})$ |
|  |  | 0.4520 |
|  |  | $-0.2257x_{i,1}x_{i,2}$ |
|  |  | $0.2303x_{i,1}^{2}x_{i,2}$ |
|  |  | $-0.2061\exp(x_{i,1})$ |
$\mathbf{G}(x_{i},x_{j})$ | $-(x_{j,1}-x_{i,1})$ | $-0.9812\left(x_{j,1}-x_{i,1}\right)$ | $-0.9477(x_{j,1}-x_{i,1})$ | $-1.0247\left(x_{j,1}-x_{i,1}\right)$
|  |  | $-0.1134\sin(x_{j,1}-x_{i,1})$ |
|  |  | $0.6586\sigma(x_{j,1}-x_{i,1})$ |
|  |  | $0.4114x_{j,1}$ |

Table 8: Numerical results for FHN dynamics (randomly delete $15\%$ edges) in
self-dynamics $\mathbf{F}(x_{i})$ and interaction dynamics
$\mathbf{G}(x_{i},x_{j})$ term by term (each term of true or inferred
$\mathbf{G}(x_{i},x_{j})$ is understood to be divided by $k_{\text{in}}$,
i.e., $\frac{\mathbf{G}(x_{i},x_{j})}{k_{\text{in}}}$, where $k_{\text{in}}$
is the in-degree of node $i$).

.

Percent of edges perturbed | Two-Phase | ARNI | FEX
---|---|---|---
$15\%$ added | 0.0082 | 0.8437 | 0.0049
$15\%$ deleted | 0.0077 | 0.7192 | 0.0067
Table 9: sMAPE of FHN dynamics with $15\%$ perturbed edges added or removed.

###  6.6 Observational noises

In this subsection, we demonstrate the robustness of our method against
observational noises. We focus on the HR dynamics on a directed SF network. We
generate noisy data through

$X_{i}^{obs}(t)=\bm{x}_{i}(t)+a\beta_{X}(t),$

where $\beta_{X}(t)$ is observational noise, which follows a Gaussian
distribution with zero mean and standard deviation one, so the parameter $a$
is used to tune the level of noise.

| True Dynamics | Two-Phase | ARNI | FEX
---|---|---|---|---
$\mathbf{F}(x_{i})$ | $x_{i,2}$ | $0.9631x_{i,2}$ | $0.6002x_{i,2}$ | $0.9916x_{i,2}$
| $-x_{i,1}^{3}$ | $-1.0039x_{i,1}^{3}$ | $-1.3202x_{i,1}^{3}$ | $-1.0023x_{i,1}^{3}$
| $3x_{i,1}^{2}$ | $2.7820x_{i,1}^{2}$ | $2.3557x_{i,1}^{2}$ | $2.9659x_{i,1}^{2}$
| $-x_{i,3}$ | $-0.9663x_{i,3}$ | $-0.7982x_{i,3}$ | $-0.9871x_{i,3}$
| $3.24$ | $2.9673$ | $1.9699$ | $3.2453$
|  |  | $0.6597\exp(x_{i,1})$ |
$\mathbf{G}(x_{i},x_{j})$ | $0.3\sigma(x_{j,1})$ | $0.1529\sigma(x_{j,1})$ | $0.2804\sigma(x_{j,1})$ | $0.2963\sigma(x_{j,1})$
| $-0.15x_{i,1}\sigma(x_{j,1})$ | $-0.0995x_{i,1}\sigma(x_{j,1})$ | $-0.0847x_{i,1}\sigma(x_{j,1})$ | $-0.1262x_{i,1}\sigma(x_{j,1})$

Table 10: Numerical results for Hindmarsh-Rose dynamics (SNR = 45) in self-
dynamics $\mathbf{F}(x_{i})$ and interaction dynamics
$\mathbf{G}(x_{i},x_{j})$ term by term.

.

| True Dynamics | Two-Phase | ARNI | FEX
---|---|---|---|---
$\mathbf{F}(x_{i})$ | $x_{i,2}$ | $0.9341x_{i,2}$ | $0.3916x_{i,2}$ | $0.9633x_{i,2}$
| $-x_{i,1}^{3}$ | $-0.9948x_{i,1}^{3}$ | $-1.0638x_{i,1}^{3}$ | $-1.0058x_{i,1}^{3}$
| $3x_{i,1}^{2}$ | $2.6599x_{i,1}^{2}$ | $2.5882x_{i,1}^{2}$ | $2.8923x_{i,1}^{2}$
| $-x_{i,3}$ | $-0.9332x_{i,3}$ | $-0.3871x_{i,3}$ | $-0.9588x_{i,3}$
| $3.24$ | $3.1341$ | $1.8715$ | $3.2342$
|  |  | $0.1318x_{i,2}x_{i,3}$ |
|  |  | $0.2246\exp(x_{i,2})$ |
$\mathbf{G}(x_{i},x_{j})$ | $0.3\sigma(x_{j,1})$ | $0.1333\sigma(x_{j,1})$ | None | $0.2683\sigma(x_{j,1})$
| $-0.15x_{i,1}\sigma(x_{j,1})$ | None | None | $-0.1118x_{i,1}\sigma(x_{j,1})$
|  |  | $-0.2341\sigma(x_{j,1}-x_{i,1})$ |

Table 11: Numerical results for Hindmarsh-Rose dynamics (SNR = 40) in self-
dynamics $\mathbf{F}(x_{i})$ and interaction dynamics
$\mathbf{G}(x_{i},x_{j})$ term by term.

.

SNR | Two-Phase | ARNI | FEX
---|---|---|---
$45$ | 0.0924 | 0.2720 | 0.0158
$40$ | 0.2190 | 0.6252 | 0.0376
Table 12: sMAPE of HR dynamics with different levels of noise.

From Table 10, Table 11 and Table 12, FEX demonstrates its capability to
accurately identify all terms within the HR dynamics with their respective
coefficients with SNR = 45 and 40. The Two-Phase method can indeed infer the
correct terms in the dynamics in the case of SNR = 45. However, the accuracy
of inferred coefficients is low, particularly in the interaction dynamics. In
addition, in the case of SNR = 40, Two-Phase method fails to infer the
$-0.15x_{i,1}\sigma(x_{j,1})$ term in the interaction dynamics. ARNI can infer
most of the terms correctly when SNR = 45, with an extra term
$0.6597\exp(x_{i,1})$ in the self-dynamics. Furthermore, the coefficients in
the self-dynamics are poorly estimated, but it is interesting to notice that
the coefficients of interaction-dynamics are more accurate than those of Two-
Phase, owing again to ARNI’s capability to infer the network structure. As for
SNR = 40, ARNI infers more redundant terms in self-dynamics, and it fails to
identify both true terms in the original dynamics, and incorrectly identifies
a $-0.2341\sigma(x_{j,1}-x_{i,1})$ term. This shows that ARNI can reasonably
well infer most of the dynamics terms in low-noise case, while it fails to
work as noise level increases. Therefore, FEX has demonstrated its superior
performance compared to Two-Phase and ARNI to handle noisy data.