Do We Need an Encoder-Decoder to Model Dynamical Systems on Networks?

1 Introduction

In deep learning, architectural design patterns proven successful for one task are often transplanted to tackle related tasks. One such design pattern is an encoder-decoder structure found in almost every competitive model, be it a transformer [Vaswani et al., 2017] or a graph neural network [Kipf and Welling, 2016]. With recent dramatic progress in neural dynamical model [Chen et al., 2018; Chamberlain et al., 2021], one would naturally assume that an encoder-decoder structure in such models is also helpful and necessary.

We challenge this assumption in the context of reconstructing network dynamics from observation time series. Many natural and social phenomena can be modelled as dynamical processes on complex networks, such as gene regulatory networks, social networks, and ecological networks [Newman, 2018]. Recovering the governing mechanism of such dynamical systems from observations is the first step towards understanding the global dynamical properties and eventually implementing interventions and controls. In a typical problem setup, a multivariate time series x(t1), x(t2), . . . , x(tN) is observed from a network dynamical system. From it, we aim to inductively infer a coupled dynamical system to model x(t). When the observations are evenly spaced, various solutions exist (e.g., [Takens, 1981; Brunton et al., 2016]). But as many complex systems generate unevenly spaced observations, recovering the dynamic remains a challenge.

∗Corresponding authors.

Many researchers exploited an encoder-decoder structure to recover coupled flows from time series data [Zang and Wang, 2020; Yan et al., 2021; Wang et al., 2022; Guo et al., 2022; Hwang et al., 2021]. Among these, the most representative model is NDCN [Zang and Wang, 2020]. As shown in Fig 1, instead of estimating a flow Φ from the original observations x, NDCN seeks to simultaneously estimate an encoder fe and a conjugate flow Φ′ such that Φ′ (ti, fe(x0)) ≃ (fe ◦ Φ)(ti, x0). As the new flow Φ′ now lives in a higherdimensional latent space, it potentially has a simpler form than Φ, in the same way that kernel methods allow linear algorithms to be used for nonlinear problems [Hofmann et al., 2008]. However, this seemingly innocent use of an encoder can result in an incorrect dynamical model, in the sense defined by Section 3.1. In particular, we observe that a flow estimated by NDCN is often dynamically unstable (with the largest Lyapunov exponent being positive), even when the time series originates from a stable flow. We analyse the potential causes. By concluding that an encoder is unnecessary, we show that it is possible and more desirable to directly model the dynamics in its original state space, by proposing an effective and versatile embedding-free architecture.

In this paper, we make the following contributions.

• We expose a severe limitation of encoder-decoder-based models for neural dynamics on networks: superfluous latent dimensions induce incorrect flows. In particular, we show that some widely used models do not even produce a valid dynamical system. We demonstrate this phenomenon on various dynamics and networks.
• We propose an encoder-free architecture Dy-Net Neural Dynamics (DNND) to learn network dynamics from time series data. DNND goes beyond mere observation

fitting and can reliably learn a model of the correct dynamical behaviours.

• We propose a loss warm-up procedure to address the numerical instability in training neural ODE on irregular time series data.
• We improve the existing evaluation protocol for neural dynamic models by including metrics aligned with common objectives of analysing nonlinear complex systems, including stability and limiting behaviours.

2 Background

2.1 Dynamics on Complex Networks

Consider an undirected network with n vertices. We represent the network with a symmetric adjacency matrix A. An element Aij in A is defined as

In this paper, we assume that the network contains no selfloop. Hence Aii = 0 for any i.

Let T ≜ R+ (or more generally a monoid). We assume that each vertex i defines a continuous-time real-value state variable xi : T → R. And we consider a dynamical system defined by ordinary differential equations together with initial conditions:

where t ∈ T and 1 ≤ i ≤ n. In this paper, we use ˙x to denote d dt x (i.e., the velocity) and call Ψi : Rn → R a velocity operator. We assume that each velocity operator is smooth and globally Lipschitz. Also, we restrict Ψi to be autonomous (time-invariant).

Let x denote (x1, x2, . . . , xn)⊤. Then the evolution operators Ψi induce a flow Φ : T × Rn → Rn as

The problem that concerns us in this paper is to recover Ψi from a multivariate time series:

where 0 ≤ t1 < t2 < · · · < tN. Note that we do not assume that the time series is regularly spaced.

Clearly, in its most general form (2), it is an ill-defined problem to recover Φ from a single time series. If we review ∆ti ≜ ti+1 - ti as the value of a roof function r(x(ti)), then Φ is a special flow [Fisher and Hasselblatt, 2019], and there are infinitely many valid special flows under time changes. Therefore, we need to introduce an additional inductive bias. Following [Zang and Wang, 2020], we assume a single operator Ψ : R × XA → R shared by all state variables:

where XA ⊂ 2R. In other words, xi(t) depends on itself and its neighbouring state variables, and Ψ is invariant to the permutation of the neighbouring vertices.

2.2 Encoder-Decoder Based Neural Dynamical Models

How do you model the operator Ψ in Eq (7)? Many neural dynamical models for complex networks (e.g., [Zang and Wang, 2020]) follow a common encoder-decoder architecture in graph neural networks [Kipf and Welling, 2016] and also in general neural ODEs [Lee and Parish, 2021]. In NDCN [Zang and Wang, 2020] (see Fig 1.(a)), the authors defined an encoder fE : R → Rd and a decoder fD : Rd → R, where d ∈ Z+ is the embedding dimension, so that each initial value xi(0) is encoded jointly as d state variables in the vector hi(0) = fE(x(0)). Then they define

where L is the graph Laplacian, H(t) is (h1(t), h2(t), . . . , hn(t))⊤, and fθ : Rd → Rd is a vector-valued function. Clearly, the use of graph Laplacian for aggregation is a convenient choice. It restricts the type of dynamical coupling among neighbouring vertices.

Training an NDCN model involves jointly fitting fE, fθ, fD to minimise the empirical risk on the observations. It is not different from a standard supervised training setup. In particular, no explicit requirement is imposed on the dynamical behaviours of the fitted flow.

3 Issues with Encoder-Decoder Based Neural Dynamical Models

In this section, we show that encoding the scalar state variable x(t) as a latent vector h(t) leads to overfitting and spurious evolution operators. We use a simulated example to demonstrate the problems and then analyse potential causes in Section 3.3. Later in Section 4, we will propose a simple yet effective solution to address these problems.

3.1 Three Tests for Neural Dynamical Models

The universal approximation theorem [Hornik et al., 1989] allows us to fit a neural network to approximate any time series. In NDCN, if we keep increasing the embedding dimension d, we can always fit a given observation time series. But for any fitted neural dynamical model to be useful, we expect the induced flow to have the correct dynamical properties, which we aim to capture with the following three tests.

Extrapolation beyond training data. To fit a neural dynamical model, we rely on observations at finite time points t1, t2, . . . , tN. These observations are often collected during a short-term transient period close to an initial condition. In many studies of NODEs, a test sample from a time period overlapping or close to the training period is used. This is consistent with the standard supervised setting, where data is assumed to be independent and identically distributed (iid) or nearly iid. However, in dynamical modelling, more important is the long-term behaviour of the fitted model. Therefore, the model needs to make accurate predictions on long-term horizons. So our first test is whether the estimated model can make accurate long-term predictions beyond the observation period of the training data.

Same fixed points or invariant sets. One may argue that in dynamic system modelling, what matters is the correct final state, instead of accurate predictions of intermediate state values. Related to the previous test is whether the fitted model needs to have the same collection of fixed points (or invariant sets) as the true model. If the answer is yes, two corresponding fixed points need to have the same dynamical stability. In this paper, we check the stability of a flow by the largest Lyapunov exponent [Wolf et al., 1985].

Neural network as a legal flow. The final test is the most important but is often overlooked: whether the fitted neural dynamical model induces a legal vector field for a flow. In particular, the model should induce a state transition function that is a monoid action. This means that if Φ is a flow, then

for any t1, t2 ∈ T . In other words, the trajectory x(t) cannot branch out at any time t.

3.2 HowNDCN Fares in These Tests

We demonstrate the problems with embedding-based neural dynamical models using an example from [Zang and Wang, 2020]. In this example, NDCN is shown to outperform various RNN+GNN baselines, in short-horizon extrapolation tasks (see Table 3 of [Zang and Wang, 2020]).

Heat diffusion on a frid network. Following [Zang and Wang, 2020], we simulate a vector time series using a heat equation defined on a grid network. We consider a network with 400 vertices arranged in a square 2d lattice. Let {v(k, l) : 1 ≤ k, l ≤ 20} be the vertices. Two vertices v(k, l) and v(k′ , l′ ) are connected if and only if |k - k′ | ≤ 1 and |l - l′ | ≤ 1. In other words, this is a modified version of the 8-regular graph with 2d boundaries. On this network, a flow is defined by the heat equation. For each node i, we have:

where α is a diffusive constant. Also xi(0) is set with a random value from [0, 25]. Note that in theory, the above coupling operator can be perfectly modelled by the graph Laplacian. From the above flow, we created 80 irregularly spaced samples with t ∈ [0, 5], and fitted an NDCN model using the code provided by the authors on Git Hub (https: //github.com/calvin-zcx/ndcn).

Failing to extrapolate. Fig 2 shows the true trajectories of 10 vertices and the corresponding trajectories recovered from the NDCN model. The NDCN model performed well in the neighbourhood of the training observations (to the left of the dashed line). But it has almost poor generalisation capacity outside that region. Fig 2(c) shows that the prediction error grows as we move away from the training data.

Missing the fixed point. The heat flow has a simple fixed point: eventually, every xi converges to the average temperature of the network. Fig 2b shows that the NDCNmodel fails to capture this very stable fixed point, as all state variables actually diverge over time in the fitted model. We measured the largest Lyapunov exponent of the true flow and the estimated flow. While the true value is -35.83, indicating a stable flow, the value of the estimated flow is 9.36, a positive value that implies instability. The same contradiction can be revealed by a 2D visualisation of a vertex-wise phase plot. Fig 8 shows that the estimated flow initially approached the average temperature but eventually diverged from the correct fixed point.

Failing to produce a flow. From the estimated flow ˆΦ, we produce two time series. One from t = 0 with the initial condition xi(0) and another from t = 5 with the initial condition ˆΦ(5, xi(0)). If ˆΦ is a legal flow, then these two time series should coincide from t ≥ 5. Figure 3 shows the time series for two vertices. Clearly resetting the initial condition led to a completely different time series.

3.3 What Went Wrong?

The above example shows that a neural dynamical model can fit an observation time series well in the embedding space, but at the same time fails to produce the correct flow. Embedding has been used to reconstruct flows from regularly spaced time series [Takens, 1981]. Also, it has been shown that any homeomorphism h : X → X,X ⊂ Rp in a pdimensional Euclidean space can be approximated by a Neural ODE ˆΦ : T × R2p → R2p operating on a 2p-dimensional Euclidean space [Zhang et al., 2020]. But then why did the embedding-based model fail to recover the heat dynamics?

Encoder/Decoder or Homeomorphism?

First, we limit our definition of a correct flow. Following [Fisher and Hasselblatt, 2019], a flow Φ = {ϕ(t, x)} on X and another flow Ψ = {ψ(t, y)} on Y are topological equivalent if there is a homeomorphism h : X → Y such that

If we accept the notion of "correctness", then in Fig 1(a) the encoder function fE and the decoder function fD have to satisfy some conditions.

First, fE and fD need to be inverse to each other, as they represent h and h-1 . This implies that fE and fD should not be implemented using two independent neural networks. Secondly, both fE and fD must be one-to-one functions (injections). This implies that the encoder/decoder neural network has to be reversible. In particular, it cannot contain pooling layers or a RELU activation function. However, the encoders and decoders used in most neural dynamical models do not meet these two conditions. In [Zang and Wang, 2020], fD and fE are parametrised using two separate neural networks, just as in a standard GNN.

Curse or Blessing of Dimensionality?

Another problem with embedding is the increased dimensionality of an inverse problem. As in NDCN, the encoder maps a scalar value to a vector of dimension d. So instead of fitting a shared scalar field as in Eq (7), we need to learn the d-dimension vector field. At the same time, we have only one time series of dimension n as the training data. When d is too large relative to n, overfitting is likely to occur. This is illustrated in Fig 4.

One may argue that embedding is routinely used by neural ODEs in classification/regression tasks [Chen et al., 2018; Chamberlain et al., 2021], and there appear to be no issues. In classification/regression tasks, neural ODEs produce a flow Φ and a return time T such that Φ(T, ·) serves as an approximator of the homeomorphism between inputs and label encoding. In such a scenario, it is well established that lifting the dimension of the input features enlarges the set of valid evolution functions, making a "difficult" homeomorphism easier to fit [Dupont et al., 2019]. However, it is unnecessary to consider whether the flow Φ itself is "correct" or not. With that said, superfluous dimensions may still lead to unstable fit [Massaroli et al., 2020].

In NDCN, the additional embedding dimension does provide additional expressive power that compensates for the limitation of the graph Laplacian aggregation function. But we argue that if the graph Laplacian is too restrictive for the task, we should use a more general aggregation function.

4 Dy-Net Neural Dynamics

In the previous section, we have seen that an embeddingbased neural dynamical model can learn a spurious flow that fits the training data well, but is qualitatively different from the true dynamic in terms of long-term properties.

If it is so difficult to learn an encoder/decoder pair that guarantees topological conjugacy, do we need latent embedding in the first place? In this section, we show the possibility of constructing neural networks that directly model the flow in its original state space.

4.1 Replacing Graph Laplacian with Permutation Invariance Aggregation

If we build a neural network to work directly on the original state space, one challenge is the limited expressive power of the graph Laplacian. The graph Laplacian and its variants play a central role in graph convolution and graph neural networks [Kipf and Welling, 2016; Veliˇckovi´c et al., 2017]. However, the graph Laplacian is a linear aggregation operator that cannot handle potentially nonlinear coupling in network dynamics. In other words, to avoid encoders and decoders, we also have to abandon the graph Laplacian and seek a more general aggregation scheme. We propose a model based on the separation of self-dynamic and coupling operators.

4.2 Dy-Net Neural Dynamics

Our first approximation is to replace Ψ(xi(t), {xj(t) : Aij = 1}) with the sum of a node operator and an edge operator:

By Theorem 2 in [Zaheer et al., 2017], ΨE({xj(t) : Aij = 1} can be represented as ρ (∑j:Aij=1 ϕ(xj )), where ρ and ϕ are two scalar transformations.

Inspired by the pairwise coupling feature in [Chen et al., 2014], we further assume that ρ is additive and that each xj can have a non-linear coupling effect on xi governed by a shared termΨe(xi(t), xj(t)). This leads to the following sum decomposition of the operator Ψ:

This turns out to be the same vector-field structure that appears in [Barzel and Barab´asi, 2013; Barzel et al., 2015; Gao et al., 2016].

Note that Eq (14) is more general than the graph Laplacian aggregation in Eq (11). It is also more general than the threeterm structure in [Barzel et al., 2015]:

Figure 5: Loss warm-up schedule for learning F and G. (a) At the early stage of the training process, the contribution of the data points far away from the initial time (t being large) are dampened to stabilise the training. (b) A temperature variable τ controls how the dampening is relaxed with the number of epochs. The temperature increases every 400 epochs.

Eq (14) suggests that we can use two separate neural networks Fθ and Gθ to parametrise ΨN and Ψe (see Fig. 1). This is in contrast to the single-neural-network solution in NDCN. Compared with an encoder-decoder-based model, a model specified by Eq (14) offers important benefits: The neural networks for ΨN and Ψe are independent of A. That means that these two neural networks Fθ,Gθ can be trained even when the adjacency matrix A changes within a time series. Dynamical network rewiring is common in many applications, such as biological and social systems. Because of this decoupling of Ψe from A, we call our two-neural-network model Dy-Net Neural Dynamics (DNND).

Fig 1(b) shows the architecture of DNND, which renders the following flow on networks:

for all 1 ≤ i ≤ n. Clearly xi(t) is always a flow by definition.

Affine regularisation. It is well known that neural ODE training can be numerically challenging due to approximation errors in numerical integration procedures and stiffness in the system [Kim et al., 2021]. To alleviate the challenge of unstable model training, we model Fθ and Gθ as

where cf0, cf1, cg0, cg11, cg12 are scalar parameters to capture the dominating linear terms, and fθ and gθ are MLPs to capture higher-order nonlinear terms. This is similar to the set-up of LassoNet [Lemhadri et al., 2021]. It introduces the inductive bias of preferring affine mappings. It allows us to regularise fθ and gθ to promote sparse and interpretable models. Note that the graph Laplacian is now a special case with only two nonzero terms cg11 = 1 and cg12 = -1.

4.3 Robust Fitting of Two Neural Networks

Given a time series D = x(t1), x(t2), . . . , x(tN), the learning problem is to find Fθ : R → R and Gθ : R2 → R that minimise the following loss:

To reduce the impact of numerical error accumulation from long-range integration [Milani and Nobili, 1987], we devise a warm-up schedule to adapt the loss function dynamically. At the early stage of training, the model parameters are strongly affected by time-series data points closer to the initial value, which results in smaller integration errors. More specifically, we adopt a weighting schedule for the loss function with evolving weight function wk at epoch k as follows:

where τk monotonically increases with the number of training epochs k. So the loss function evolves as follows:

5 Experiments

We use several approaches to recover evolution operators from time series and then compare the quality of the recovered models. We closely follow [Zang and Wang, 2020] to set up the ground-truth models and to simulate the observation data. In this setup, a ground-truth model is a coupled system of ODEs defined on a random network. From initial conditions, we simulate a time series as the observation data. Our goal is to recover a neural-network model for the evolution operators in the ground-truth ODE system.

5.1 Ground-truth Dynamical Models and Observation Data

We considered an extensive set of linear and nonlinear ODEs (Table 1) on different complex networks, including a fixed grid network from [Zang and Wang, 2020] and random networks based on the Erd´os and R´enyi (ER) model [Gilbert, 1959], the Albert-Barab´asi (BA) model [Barab´asi and Albert, 1999], the Watts-Strogatz (WS) model [Watts and Strogatz, 1998], and the LFR model [Fortunato, 2010; Lancichinetti et al., 2008].

Training Time Series. We follow the setup in [Zang and Wang, 2020] and randomly sample 80 (irregularly spaced) times 0 ≤ t1 < t2 < · · · < t80 ≤ 5 for training. The initial conditions x(0) of the dynamical variables on vertices are set with random values from [0, 25]. We use the Runge-Kutta-Fehlberg (RKF45) solver to generate x(t1), x(t2), . . . , x(t80) as the training time series.

5.2 Generalisation beyond the Training Data

Table 2 compares the prediction performance from the shortterm prediction to the long-term prediction. One can see that NDCN performs reasonably well near the training data, but its long-term prediction is unreliable, with the mean absolute percentage error (MAPE) even exceeding 100%. In contrast, DNND produces accurate predictions, even well beyond the region of training data. The pattern persists on other network topologies (see more results in the supplementary material).

Fig. 7 shows time series on 50 vertices. NDCN produces time series with values far exceeding the physically reasonable values.

5.3 Recovery of Dynamic

The overfitting of NDCN is also revealed in the largest Lyapunov exponent of the fitting model. Table 3 shows that all dynamics are stable with the grid network, with the negative Lyapunov exponents. The DNND estimates are similar to the true values. In contrast, NDCN produces positive Lyapunov exponents, suggesting unstable dynamics.

we can also see the fixed point, if any, of the fitted model using a projection into the 2D space spanned by xi and ˙xi. Fig. 8 shows that DNND converges to the correct fixed point. In contrast, NDCN does not converge at all.

Recovering Governing Evolution Function. As shown in [Murphy et al., 2022], predictability and reconstructability are independent and sometimes incompatible problem properties. A time series with less variability is easier to predict but provides less information for recovering the evolution function. Fig. 9 shows that from the two subnetworks in DNND, we can reasonably recover the vector field in terms of self-dynamic and coupling effect: ΨN and Ψe in Eq. (14).

Note that NDCN's vector field is defined in a latent space and it is impossible to recover the evolution function in terms of the original state variables. Unlike NDCN, which models the governing function on latent variables, DNND directly captures the governing function of the dynamic in terms of the original physical measurements. From the Fθ-network and Gθ-network, we can recover the input-output dependency. It is another advantage of our embedding-free approach.

Furthermore, if ΨN and Ψe have analytical forms in terms of basis functions, then symbolic regression can potentially be applied to recover the exact algebraic forms of the two functions. We show some results in the Appendix.

6 Conclusion

We show that blindly applying an encoder-decoder structure, without a sound principle such as topological conjugacy, is harmful to correct recovery of a dynamic flow from observations. Through our proposed Dy-Net Neural Dynamics (DNND) model, we show that state embedding is unnecessary for a broad class of flow-recovery problems. Empirically, DNND reliably recovers a wide range of dynamical systems on different network topologies.

As DNND models the vector field in the original state space, the fitted model can serve as a basis for discovering analytical, hence interpretable, differential equations using techniques such as symbolic regression [Brunton et al., 2016]. We presented some early results in the Appendix.

Theoretically, it would be interesting to characterise the learnability of dynamical flows, up to topological equivalence, from irregularly spaced observations. This is related to finding a suspension flow under a roof function [Fisher and Hasselblatt, 2019]. But as we have access to only samples of the homeomorphism, it is a machine learning problem. One particularly interesting question worth exploring is how to parametrise time changes ( ∂ ∂t ) given a finite sample from an unknown roof function.

Acknowledgments

This work was supported by the National Key R&D Program of China under Grant Nos. 2021ZD0112501 and 2021ZD0112502; the National Natural Science Foundation of China under Grant Nos. U22A2098, 62172185, 62206105 and 62202200; Jilin Province Capital Construction Fund Industry Technology Research and Development Project No.

• 2022C047-1; Changchun Key Scientific and Technological Research and Development Project under Grant No. 21ZGN30; Key Laboratory of Symbolic Computation and Knowledge Engineering of Ministry of Education Open Research Project No. 93K172021K05.