Learning Stochastic Dynamical Systems as an Implicit Regularization with
  Graph Neural Networks
Stochastic Gumbel graph networks are proposed to learn high-dimensional time
series, where the observed dimensions are often spatially correlated. To that
end, the observed randomness and spatial-correlations are captured by learning
the drift and diffusion terms of the stochastic differential equation with a
Gumble matrix embedding, respectively. In particular, this novel framework
enables us to investigate the implicit regularization effect of the noise terms
in S-GGNs. We provide a theoretical guarantee for the proposed S-GGNs by
deriving the difference between the two corresponding loss functions in a small
neighborhood of weight. Then, we employ Kuramoto's model to generate data for
comparing the spectral density from the Hessian Matrix of the two loss
functions. Experimental results on real-world data, demonstrate that S-GGNs
exhibit superior convergence, robustness, and generalization, compared with
state-of-the-arts.
Jin Guo, Ting Gao, Yufu Lan, Peng Zhang, Sikun Yang, Jinqiao Duan
# Learning Stochastic Dynamical Systems as an Implicit Regularization with Graph Neural Networks

###### Abstract

Stochastic Gumbel graph networks are proposed to learn high-dimensional time series, where the observed dimensions are often spatially correlated. To that end, the observed randomness and spatial-correlations are captured by learning the drift and diffusion terms of the stochastic differential equation with a Gumble matrix embedding, respectively. In particular, this novel framework enables us to investigate the implicit regularization effect of the noise terms in S-GGNs. We provide a theoretical guarantee for the proposed S-GGNs by deriving the difference between the two corresponding loss functions in a small neighborhood of weight. Then, we employ Kuramoto's model to generate data for comparing the spectral density from the Hessian Matrix of the two loss functions. Experimental results on real-world data, demonstrate that S-GGNs exhibit superior convergence, robustness, and generalization, compared with state-of-the-arts.

## 1 Introduction

Multivariate time series enable us to thoroughly investigate the statistical pattern among the observed dimensions, and to make predictions. Over last decades, an amount of works have been dedicated to modeling time series data arising in applications including finance, biology, and etc. Real-world time series data often exhibit a certain amount of correlations or causal relationships between the oberved dimensions. For instance, the prices of many financial derivatives, could be simultaneously influenced by the common market signals and interventions. The dynamics of electroencephalogram (EEG) brain signals, implicitly reflect a latent graph structure [1]. Hence, it is crucial to exploit the graph structure in modelling multivariate time series, within a dynamical system.

Graph neural networks (GNN) [2] aim to extract both local and global features, by leveraging available graph structure, using neural networks. Canonical GNNs, including graph convolutional networks (GCN) [3] and graph attention networks (GAT) [4], have demonstrated their strong capacity in capturing graph-structured data. Recently, there has been growing interests in predicting time series with graph neural networks to capture the underlying graph structure. For instance, temporal graph convolutional neural networks (T-GCN) [5] can integrate temporal information with the available graph structure. Nonetheless, the graph information, which may characterize the correlations or causal relations between observed dimensions, are either not immediately available, or contain noise.

Stochastic dynamical system, as a mathematical tool [6] to describe physical system evolving over time, becomes more and more popular nowadays to model various real world complex phenomena. Stochastic dynamical system [7], has been receiving increasing attention specifically in machine learning domain [8, 9], because of its mathematical rigor and strong modelling flexibility. Besides, it is also a bridge connecting real data with deep learning algorithms [10, 11]. On the one hand, many kinds of
neural network has the corresponding continuous version of differential equations, which helps building some convergence guaranteed neural networks such as [12, 13, 14, 15]. On the other hand, mathematical analysis from dynamical system point of view could also provide insights on loss landscape and critical points [16]. Furthermore, investigation on Edge of Stability (EOS) phenomenon are also promising research directions [17, 18].

Moreover, the Gumbel graph neural network (GGN) [19], is advanced to recover the graph structure underlying time series, within a dynamical system. By effectively leveraging the graph structure, the GGN achieves better accuracy in predicting high-dimensional time series. Nevertheless, the GGN suffer from the problems of limited generalization capabilities and excessive smoothness stemming from the inherent properties of graph neural networks themselves. To address these concerns, this paper proposes the stochastic Gumbel Graph Network (S-GGN) model by introducing a diffusion term. The main contributions of this paper can be summarized as follows:

* A stochastic Gumbel graph network (S-GGN) model is proposed to improve the model robustness and generalization ability in capturing high-dimensional time series, by introducing the noise term. In particular, we thoroughly study the convergence of the S-GGN model with theoretical analysis (Sec.2).
* A grouped convolution S-GGN structure, is advanced to capture _noisy_ graph-structured time series. Using convolution operations, the model can effectively reconstruct the dynamics by leveraging the external node features to remove the noise effects.
* Experiments on real-world problems [20], demonstrate the superior generalization capability of the S-GGN model without compromising its accuracy, compared with GCNs (Sec.3).

## 2 The S-GGN Frameworks

### S-GGN model

The GGN, consisting of a network generator and a dynamics learner, recovers the underlying dynamical systems from observations such as high-dimensional time series data. The network generator within GGN utilizes the reparameterization technique known as the Gumbel-softmax trick [21], whereby the graph is sampled based on probabilities. The application of the Gumbel-softmax trick allows the GGN to directly apply the backpropagation algorithm to calculate the gradient and optimize the network. Based on the connection between the GGN model and the discrete representation of the dynamical system, we extend its applicability to a formulation that aligns with the discrete form of the stochastic dynamical system, called Stochastic Gumbel Graph Neural Network (S-GGN). Thus, the dynamic learner of S-GGN can be represented as

\[X^{t+1}_{predict}=X^{t}+f(X^{t},A)\Delta t+\sigma(X^{t},A)\xi_{t}\sqrt{\Delta t} \tag{1}\]

where \(X^{t}\) denotes the state vector of all \(N\) nodes at moment \(t\) and \(A\) is the symmetric adjacency matrix constructed by the network generator. Here \(\xi_{t}\sim\mathcal{N}(0,I)\) is an independent standard normal random vector.

The graph neural network module within S-GGN can be depicted as a composition of the following mappings:

\[\begin{split}& H^{t}_{e_{1}}=f_{v\to e}(X^{t}\otimes(X^{t})^{T}), \\ & H^{t}_{e_{2}}=f_{e}(H^{t}_{e_{1}}),\\ & H^{t+1}_{v_{1}}=f_{e\to v}(A*H^{t}_{e_{2}}),\\ & H^{t+1}_{v_{2}}=f_{v}(H^{t+1}_{v_{1}}),\end{split} \tag{2}\]

where \(X^{t}\in\mathbf{R}^{N\times d_{x}}\) denotes the features of \(N\) nodes with each node has feature dimension \(d\). Let \(A\) denotes the adjacency matrix that describes the relationships between the nodes. The \(f_{v\to e},f_{e},f_{e\to v},f_{v}\) consist of a linear layer and an activation function. The operation \(\otimes\) signifies pairwise concatenation, and the symbol \(*\) denotes multiplication by elements. The composition of these four mappings in (2) corresponds to the function \(f\), the same networks' structures as \(\sigma\) in equation (1).

Next, we train the two neural networks for functions \(f\) and \(\sigma\) in (1), denoted as \(f_{NN}\) and \(\sigma_{NN}\) respectively. Their networks structures are the same with different values.
### Spectral Analysis

We denote the parameters of GGN networks as \(\mathbf{\omega}_{t}\) at the \(t\)-th interaction. Consider the division \(0:=t_{0}<t_{1}<\cdots<t_{M}:=T\) of \([0,T]\), we define that \(\delta_{m}:=t_{m+1}-t_{m}\). The discretization of parameters' evolution in GGN network can be expressed as

\[\mathbf{\omega}_{t+1}=\mathbf{\omega}_{t}+f_{NN}(\mathbf{\omega}_{t},X^{t})\delta_{t},\]

where \(f_{NN}:\mathbb{R}^{d_{\mathbf{\omega}}}\times\mathbb{R}^{d_{x}}\). After introducing noise \(\varepsilon\) to our networks, we consider a rescaling of the noise \(\sigma_{NN}\mapsto\varepsilon\sigma_{NN}\), then the following discretization stochastic differential equation (SDE) holds

\[\mathbf{\omega}_{t+1}^{\varepsilon}=\mathbf{\omega}_{t}^{\varepsilon}+f_{NN}(\mathbf{ \omega}_{t}^{\varepsilon},X_{t})\delta_{t}+\varepsilon\sigma_{NN}(\mathbf{\omega} _{t}^{\varepsilon},X^{t})\xi_{t}\sqrt{\delta_{t}}\,, \tag{3}\]

where \(\sigma_{NN}:\mathbb{R}^{d_{\mathbf{\omega}}}\times\mathbb{R}^{d_{x}\times r}\), and \(\xi_{t}\sim\mathcal{N}(0,I)\) is an independent \(r\)-dimensional standard normal random vector. That is the evolution of parameters in our S-GGN network. Besides, we have the \(\mathbf{\omega}_{0}^{\varepsilon}=\mathbf{\omega}_{0}\).

Then we give some conditions on drift and diffusion term. To enhance readability, we denote the network functions \(f_{NN}\) and \(\sigma_{NN}\) as \(f\) and \(\sigma\) respectively.

**Assumption 2.1**: _We assume that the drift term \(f\) and diffusion term \(\sigma\) satisfy_

_(i) For all \(t\in[0,T]\) and \(X\in\mathbb{R}^{d_{x}}\), the maps \(\mathbf{\omega}\mapsto f(\mathbf{\omega},X)\) and \(\mathbf{\omega}\mapsto\sigma(\mathbf{\omega},X)\) have Lipschitz continuous partial derivatives in each coordinate up to order three (inclusive)._

_(ii) For any \(\mathbf{\omega}\in\mathbb{R}^{d_{\mathbf{\omega}}},t\mapsto f(\mathbf{\omega},X)\) and \(t\mapsto\sigma(\mathbf{\omega},X)\) are bounded and Borel measurable on \([0,T]\)._

Under the above assumption, we define the distinct between the loss of S-GGN and GGN networks as

\[\mathbf{\mathcal{D}(\mathbf{\omega})}:=\mathbb{E}[l_{S-GGN}(\mathbf{\omega}_{M}^{ \varepsilon})-l_{GGN}(\mathbf{\omega}_{M})], \tag{4}\]

where \(l\) denotes the loss function.

**Proposition 2.2**: _(Comparison of the noise induced loss and the deterministic loss) Under Assumption 2.1, the following holds_

\[\mathbf{\mathcal{D}(\mathbf{\omega})}=\frac{\varepsilon^{2}}{2}[\hat{R}(\mathbf{\omega})- \hat{S}(\mathbf{\omega})]+\mathcal{O}(\varepsilon^{3}), \tag{5}\]

_as \(\varepsilon\to 0\), where the \(\hat{R}\) and \(\hat{S}\) represent_

\[\hat{R}(\mathbf{\omega}) =(\nabla l(\mathbf{\omega}_{M}))^{T}\sum_{k=1}^{M}\delta_{k-1}\hat{ \Phi}_{M-1,k}\sum_{m=1}^{M}\delta_{m-1}\mathbf{v}_{m},\] \[\hat{S}(\mathbf{\omega}) =\sum_{m=1}^{M}\delta_{m-1}\text{tr}(\sigma_{m-1}^{T}\hat{\Phi}_{ M-1,m}^{T}H_{\mathbf{\omega}_{M}l}\hat{\Phi}_{M-1,k}\sigma_{m-1}),\]

_with \(\hat{\Phi}_{m,k}:=\hat{J}_{m}\hat{J}_{m-1}\cdots\hat{J}_{k}\), the state-to-state Jacobians \(\hat{J}_{m}=I+\delta_{m}\frac{\partial f}{\partial\mathbf{\omega}}(\mathbf{\omega}_{ m},X_{m})\) and the \(\mathbf{v}_{m}\) is a vector with the \(p\)-th component (\(p=1,\cdots,d_{\mathbf{\omega}}\)):_

\[[\mathbf{v}_{m}]^{p}=\text{tr}(\sigma_{m-1}^{T}\hat{\Phi}_{M-2,m}^{T}H_{\mathbf{ \omega}}[f_{M}]^{p}\hat{\Phi}_{M-2,m}\sigma_{m-1}).\]

_Moreover, we have_

\[|\hat{R}(\mathbf{\omega})|\leq C_{R}\Delta^{2},|\hat{S}(\mathbf{\omega})|\leq C_{S}\Delta \tag{6}\]

_for \(C_{R},C_{S}>0\) independent of \(\Delta\), where \(\Delta:=\text{max}_{m\in\{0,1,\cdots,M-1\}}\delta_{m}\)._

**Proof.** We refer a proof similar to [22]. First, we can apply a Taylor expansion to \(\omega_{t}\), the drift and the diffusion coefficients at a small neighbourhood of \(\omega_{0}\). With Ito formula and comparing the corresponding terms of \(\varepsilon\) in the two sides of equation (3), we can obtain the result in 3. Next, in conjunction with Lemma 1, 2 and 3 in the literature [22], it can be demonstrated that (6) holds. \(\blacksquare\)

**Remark 2.3**: _Proposition 2.2 indicates that introducing noise to the state of a deterministic graph can be considered, on average, as an approximation of a regularized objective functional._
Experiments

We conduct two experiments to verify the performance of our S-GGN.

### Kuramoto Model

The Kuramoto model is a nonlinear model describing the interaction and synchronization of oscillator groups:

\[\frac{d\theta_{i}}{dt}=\omega_{i}+K\sum_{j\neq i}A_{ij}sin(\theta_{j}-\theta_{i}),i=1,2,...,N,\]

where \(\theta_{i}\) denotes the \(i\)-th oscillator phase, \(\omega_{i}\) denotes its natural frequency, \(N\) denotes the number of vibrator, \(K\) denotes coupling strength which measures the strength of the interaction between the oscillator. Here \(A_{ij}\in\{0,1\}\) are the elements of \(N\times N\) adjacency matrix. The model takes into account the phase differences of the oscillators and the interactions between them to explain the synchronization phenomenon.

Here, we need to initialize the corresponding initial phase and natural frequency for each oscillator, and then calculate the phase of the oscillator attractively according to the above formula.

* Data preparation: the numerical solution of the Kuramoto model by the fourth-order Longe-Kutta method.
* Data pre-processing: the sin value and frequency of the phase value at the corresponding time, and these two characteristics are taken as the characteristics of the node at the corresponding time. Set the window length to 20.
* Experiment settings: the optimizer of network generator and dynamic learner, the number of iteration steps is 3 and 7 respectively.

For both models, the Hessian Matrix of the empirical loss with respect to weights can be obtained every 10 epochs, and the corresponding eigenvalues over epochs is shown in the Figure 1. Observing that S-GGN has to smaller largest eigenvalue, which indicates that S-GGN can find flatter optimal weights, allowing its better performance from sharpness awareness point of view. Figure 2 shows the

Figure 1: Eigenvalues of the Hessian matrix for GGN (left) and S-GGN (right).

Figure 2: The spectral distribution for GGN (left) and S-GGN (right).
distribution of eigenvalues of the two models after the first, 50-th and 100-th epochs. We can see that the eigenvalues' concentration of S-GGN is much stronger than that of GGN, suggesting different convergence of the two models.

### Wireless communication data

In this experiment, the data is obtained from channel measurements in real-world scenarios. Consider a typical line-of-sight (LOS) scenario with a 24-millisecond spacing between points, and an 8x1 uniform linear array (ULA) at the transmitter end.

* Data standardization: the wireless communication signal data from each base station is decomposed into real and imaginary parts for normalisation.
* Data preparation: time window is chosen as 72 points. Within each time window, features are extracted using group convolution applied to data with a window length of 36.
* Experimental settings: The Adam optimizer is selected to optimize both the network generator and the dynamic learner, with 3 and 12 iteration steps respectively.

Due to the highly noisy, non-linear and non-smooth nature of Wireless communication data, direct modelling of the raw signal is challenging.Convolutional neural networks offer significant advantages in feature extraction, automatically learning local features and retaining spatial structure information. We take a rolling prediction approach wherein a dataset of length \(W\) is employed to forecast sequences ranging from \(W+1\) to \(W+p\). Notably, the actual data ranging from \(W+1\) to \(W+p\) is not incorporated into the model during this process. The construct is illustrated in Figure 3.

The S-GNN model is employed to predict the Wireless communication data. Figure 4 presents the mean square error and mean absolute error of both the S-GGN and GGN models on the test set. Notably, the S-GGN model demonstrates a smaller error and exhibits superior generalization performance compared to the GGN model.

Figure 4: The MSE (left) and MAE (right) for GGN and S-GGN in test set.

Figure 3: The process of S-GGN Prediction.
Figure 5: Predictions for 8 base stations signals.
Figure 5 exhibits the comparative prediction outcomes of the GGN and the S-GGN models, utilizing a training sample on the test set. Both the real and imaginary components of the signals originating from the eight node base stations are plotted. The solid blue line represents the actual data, while the solid yellow line corresponds to the prediction results of the S-GGN model, and the solid green line represents the prediction results of the GGN model. We can see the S-GGN outperforms and its prediction result is more close to the true data.

## 4 Conclusion

Considering the complexity of the noise and underlying dynamics of the data, we bring Stochastic Dynamical Systems as a tool to address this problem. First, by comparing the loss of S-GGN and GGN, we can see the term of the multiplicative noise can be treated as a regularization term for the perturbation in a small neighborhood of neural network's weights. As a result, the S-GGN model could achieve better generalization capabilities and robustness on noisy data. Second, we aim to explore the spectral density over iteration steps for Hessian Matrix eigenvalues of the empirical loss w.r.t. weights. We conduct experiments on data from Kuramoto model to verify the effectiveness of S-GGN. Finally, in real-world applications such as wireless communication data, we introduce group convolution techniques as our data preprocessing, which helps us to get better long-term prediction results. Despite this, there are still some problems that need to be solved which we would like to further notice, such as loss analysis through a sharpness awareness point of view, more applications in complex spatial-temporal data like EEG signals in the brain, financial data, molecular dynamics, climate forecasting and so on.

## Acknowledgments

This work was supported by the National Key Research and Development Program of China (No. 2021ZD0201300), the National Natural Science Foundation of China (No. 12141107), the Fundamental Research Funds for the Central Universities (5003011053).

## Data availability

The data that support the findings of this study are available in GitHub at [https://github.com/xiaolangege/sggn](https://github.com/xiaolangege/sggn).
