Newton-Cotes Graph Neural Networks: On the Time Evolution of Dynamic
  Systems
Reasoning system dynamics is one of the most important analytical approaches
for many scientific studies. With the initial state of a system as input, the
recent graph neural networks (GNNs)-based methods are capable of predicting the
future state distant in time with high accuracy. Although these methods have
diverse designs in modeling the coordinates and interacting forces of the
system, we show that they actually share a common paradigm that learns the
integration of the velocity over the interval between the initial and terminal
coordinates. However, their integrand is constant w.r.t. time. Inspired by this
observation, we propose a new approach to predict the integration based on
several velocity estimations with Newton-Cotes formulas and prove its
effectiveness theoretically. Extensive experiments on several benchmarks
empirically demonstrate consistent and significant improvement compared with
the state-of-the-art methods.
Lingbing Guo, Weiqing Wang, Zhuo Chen, Ningyu Zhang, Zequn Sun, Yixuan Lai, Qiang Zhang, Huajun Chen
# Newton-Cotes Graph Neural Networks:

###### Abstract

Reasoning system dynamics is one of the most important analytical approaches for many scientific studies. With the initial state of a system as input, the recent graph neural networks (GNNs)-based methods are capable of predicting the future state distant in time with high accuracy. Although these methods have diverse designs in modeling the coordinates and interacting forces of the system, we show that they actually share a common paradigm that learns the integration of the velocity over the interval between the initial and terminal coordinates. However, their integrand is constant w.r.t. time. Inspired by this observation, we propose a new approach to predict the integration based on several velocity estimations with Newton-Cotes formulas and prove its effectiveness theoretically. Extensive experiments on several benchmarks empirically demonstrate consistent and significant improvement compared with the state-of-the-art methods.

## 1 Introduction

Reasoning the time evolution of dynamic systems has been a long-term challenge for hundreds of years [1; 2]. Despite that the advances in computer manufacturing nowadays make it possible to simulate long trajectories of complicated systems constrained by multiple force laws, the computational cost still imposes a heavy burden on the scientific communities [3; 4; 5; 6; 7; 8].

Recent graph neural networks (GNNs) [9; 10; 11; 12; 13]-based methods provide an alternative solution to predict the future states directly with only the initial state as input [14; 15; 16; 17; 18; 19; 20]. Take molecular dynamics (MD) [4; 21; 22; 23; 24] as an example, the atoms in a molecule can be regarded as nodes with different labels, and thus can be encoded by a GNN. As the input and the target are atomic coordinates, the learning problem becomes a complicated regression task of predicting the coordinate at each dimension of each atom. Furthermore, some directional information (e.g., velocities or forces) are also important features and cannot be directly leveraged especially considering the rotation and translation of molecules in the system. Therefore, the existing works put great effort into the reservation of physical symmetries, e.g., transformation equivariance and geometric constraints [18; 19]. The experimental results on a variety of benchmarks also demonstrate the advantages of their methods.

Without loss of generality, the main target of the existing works is to predict the future coordinate \(\mathbf{x}^{T}\) distant to the given initial coordinate \(\mathbf{x}^{0}\), with \(\mathbf{x}^{0}\) and the initial velocity \(\mathbf{v}^{0}\) as input. Then, the
predicted coordinate \(\hat{\mathbf{x}}^{T}\) can be written as:

\[\hat{\mathbf{x}}^{T}=\mathbf{x}^{0}+\hat{\mathbf{v}}^{0}T. \tag{1}\]

For a complicated system comprising multiple particles, predicting the velocity term \(\hat{\mathbf{v}}^{0}\) rather than the coordinate \(\hat{\mathbf{x}}^{T}\) improves the robustness and performance. Specifically, by subtracting the initial coordinate \(\mathbf{x}^{0}\), the learning target can be regarded as the normalized future coordinate, i.e., \(\hat{\mathbf{v}}^{0}=(\mathbf{x}^{T}-\mathbf{x}^{0})/(T-0)\), which is why all state-of-the-art methods adopt this strategy [18; 19]. Nevertheless, the current strategy still has a significant drawback from the view of numerical integration.

As illustrated in Figure 1, supposed that the actual velocity of a particle is described by the curve \(v(t)\). To predict the future state \(\mathbf{v}^{T}\) at \(t=T\), the existing methods adopt a constant estimation approach, i.e., \(\hat{\mathbf{x}}^{T}=\mathbf{x}^{0}+\hat{\mathbf{v}}^{0}T=\mathbf{x}^{0}+ \int_{0}^{T}\hat{\mathbf{v}}^{0}dt\). Evidently, the prediction error could be significant as it purely relies on the fitness of neural models. If we alternatively choose a simple two-step estimation (i.e., Trapezoidal rule [25]), the prediction error may be reduced to only the blue area. In other words, the model only needs to _compensate_ for the blue area.

In this paper, we propose Newton-Cotes graph neural networks (abbr. NC) to estimate multiple velocities at different time points and compute the integration with Newton-Cotes formulas. Newton-Cotes formulas are a series of formulas for numerical integration in which the integrand (in our case, \(v(t)\)) are evaluated at equally spaced points. One most important characteristic of Newton-Cotes formulas is that the integration is computed by aggregating the values of \(v(t)\) at these points with a group of weights \(\{w^{0},...,w^{k}\}\), where \(k\) refers to the order of Newton-Cotes formulas and \(\{w^{0},...,w^{k}\}\) is irrelevant to the integrand \(v(t)\) and can be pre-computed by Lagrange interpolating polynomial [26].

We show that the existing works can be naturally derived to the basic version NC (\(k=0\)) and theoretically prove that the prediction error of this estimation as well as the learning difficulty will continually reduce as the increase of estimation step \(k\).

To better train a NC (k), we may also need the intermediate velocities as additional supervised data. Although these data can be readily obtained (as the sampling frequency is much higher than our requirement), we argue that the performance suffers slightly even if we do not use them. The model is capable of learning to generate promising velocities to better match the true integration. By contrast, if we do provide these data to train a NC (k), denoted by NC\({}^{+}\) (k), we will obtain a stronger model that not only achieves high prediction accuracy in conventional tasks, but also produces highly reliable results of long-term consecutive predictions.

We conduct experiments on several datasets ranging from N-body systems to molecular dynamics and human motions [27; 28; 29], with state-of-the-art methods as baselines [17; 18; 19]. The results show that NC improves all the baseline methods with a significant margin on all types of datasets, even without any additional training data. Particularly, the improvement for the method RF [17] is greater than \(30\%\) on almost all datasets.

## 2 Related Works

In this section, we first introduce the related works using GNNs to learn system dynamics and then discuss more complicated methods that possess the equivariance properties.

### Graph Neural Networks for Reasoning Dynamics

Interaction network [30] is perhaps the first GNN model that learns to capture system dynamics. It separates the states and correlations into two different parts and employs GNNs to learn their interactions. This progress is similar to some simulation tools where the coordinate and velocity/force are computed and updated in an alternative fashion. Many followers extend this idea, e.g., using

Figure 1: Illustration of different estimations. The rectangle, trapezoid, and striped areas denote the basic (used in the existing works), two-step, and true estimations, respectively. Blue denotes the error of two-step estimation that a model needs to compensate for, whereas blue \(+\) yellow denotes the error of the existing methods.
hierarchical graph convolution [9; 10; 31; 32] or auto-encoder [33; 34; 35; 36] to encode the objects, or leveraging ordinary differential equations for energy conservation [37].

The above methods usually cannot process the interactions among particles in the original space (e.g., Euclidean space). They instead propose an auxiliary interaction graph to compute the interactions, which is out of our main focus.

### Equivariant Graph Neural Networks

Recently, some methods considering physical symmetry and Euclidean equivariance have been proposed and achieved state-of-the-art performance in modeling system dynamics [14; 15; 16; 17; 18; 19; 38]. Specifically, [14; 15; 16] force the translation equivariance, but the rotation equivariance is overlooked. TFN [39] further leverages spherical filters to achieve rotation equivariance. SE(3) Transformer [29] proposes to use the Transformer [13] architecture to model 3D cloud points directly in Euclidean space. EGNN [18] simplifies the equivariant graph neural network and make it applicable in modeling system dynamics. GMN [19] proposes to consider the physical constraints (e.g., sticks and hinges sub-structures) widely existed in systems, which achieves the state-of-the-art performance while maintaining the same level of computational cost. SEGNN [38] extends EGNN with steerable vectors to model the covariant information among nodes and edges.

Note that, not all EGNNs are designed to reason system dynamics, but the best-performing ones (e.g., EGNN [18] and GMN [19]) in this sub-area usually regard the velocity \(\mathbf{v}\) as an important feature. These methods can be set as the backbone model in NC, and the models themselves belong to the basic NC (0).

## 3 Methodology

We start from preliminaries and then take the standard EGNN [18] as an example to illustrate how the current deep learning methods work. We show that the learning paradigm of the existing methods can be derived to the simplest form of numerical integration. Finally, we propose NC and theoretically prove its effectiveness in predicting future states.

### Preliminaries

We suppose that a system comprises \(N\) particles \(p_{1},p_{2},...,p_{N}\), with the velocities \(\mathbf{v}_{1},\mathbf{v}_{2},...,\mathbf{v}_{N}\) and the coordinates \(\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{N}\) as states. The particles can also carry various scalar features (e.g., mass and charge for atoms) which will be encoded as embeddings. We follow the corresponding existing works [17; 18; 19] to process them and do not discuss the details in this paper.

Reasoning DynamicsReasoning system dynamics is a classical task with a very long history [1]. One basic tool for reasoning or simulating a system is numerical integration. Given the dynamics (e.g., Langevin dynamics, well-used in MD) and initial information, the interaction force and acceleration for each particle in the system can be estimated. However, the force is closely related to the real-time coordinate, which means that one must re-calculate the system states for every very small time step (i.e., \(dt\)) to obtain a more accurate prediction for a long time interval. For example, in MD simulation, the time step \(dt\) is usually set to \(50\) or \(100\) fs (\(1\) fs \(=10^{-15}\) second), whereas the total simulation time can be \(1\) ms (\(10^{-3}\) second). Furthermore, pursuing highly accurate results usually demands more complicated dynamics, imposing heavy burdens even for super-computers. Therefore, leveraging neural models to directly predict the future state of the system has recently gained great attention.

### Egnn

Although the existing EGNN methods have great advantage in computational cost over the conventional simulation tools, the potential prediction error that a neural model needs to compensate for is also huge. Take the standard EGNN [18] as an example, the equivariant graph convolution can be written as follows:

\[\mathbf{m}_{ij}^{0}=\phi_{e}(\mathbf{h}_{i},\mathbf{h}_{j},||\mathbf{x}_{i}^{0 }-\mathbf{x}_{j}^{0}||^{2},e_{ij}), \tag{2}\]

where the aggregation function \(\phi_{e}\) takes three types of information as input: the feature embeddings \(\mathbf{h}_{i}\) and \(\mathbf{h}_{j}\) for the input particle \(p_{i}\) and an arbitrary particle \(p_{j}\), respectively; the Euclidean distance
\(||\mathbf{x}_{i}^{0}-\mathbf{x}_{j}^{0}||^{2}\) at time \(t=0\); and the edge attribute \(e_{ij}\). Then, the predicted velocity and future coordinate can be calculated by the following equation:

\[\hat{\mathbf{v}}_{i}^{0} =\phi_{v}(\mathbf{h}_{i})\mathbf{v}_{i}^{0}+\frac{1}{N-1}\sum_{j \neq i}{(\mathbf{x}_{i}^{0}-\mathbf{x}_{j}^{0})\mathbf{m}_{ij}^{0}}, \tag{3}\] \[\hat{\mathbf{x}}_{i}^{T} =\mathbf{x}_{i}^{0}+\hat{\mathbf{v}}_{i}^{0}T, \tag{4}\]

where \(\phi_{v}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{1}\) is a learnable function. From the above equations, we can find that the predicted velocity \(\hat{\mathbf{v}}_{i}^{0}\) is also correlated with three types of information: 1) the initial velocity \(\mathbf{v}_{i}^{0}\) as a basis; 2) the pair-wise Euclidean distance \(||\mathbf{x}_{i}^{0}-\mathbf{x}_{j}^{0}||^{2}\) (i.e., \(\mathbf{m}_{ij}^{0}\)) to determine the amount of force (analogous to Coulomb's law); and the pair-wise relative coordinate difference \((\mathbf{x}_{i}^{0}-\mathbf{x}_{j}^{0})\) to determine the direction. For multi-layer EGNN and other EGNN models, \(\hat{\mathbf{v}}_{i}^{0}\) is still determined by these three types of information, which we refer interested readers to Appendix A for details.

**Proposition 3.1** (Linear mapping).: _The existing methods learn a linear mapping from the initial velocity \(\mathbf{v}^{0}\) to the average velocity \(\mathbf{v}^{t^{*}}\) over the interval \([0,T]\)._

Proof.: Please see Appendix B.1 

As the model makes prediction only based on the state and velocity at \(t=0\), we assume that \(\mathbf{v}^{t^{*}}\) fluctuates around \(\mathbf{v}^{0}\) and follows a normal distribution denoted by \(\mathcal{N}_{NC(0)}=(\mathbf{v}^{0},\sigma_{NC(0)}^{2})\). Then, the variance term \(\sigma_{NC(0)}^{2}\) reflects how difficult to train a neural model on data sampled from this distribution. In other words, the larger the variance is, the worse the expected results are.

**Proposition 3.2** (Variance of \(\mathcal{N}_{NC(0)}\)).: _Assume that the average velocity \(\mathbf{v}^{t^{*}}\) over \([0,T]\) follows a normal distribution with \(\mu=\mathbf{v}^{0}\), then the variance of the distribution is:_

\[\sigma_{NC(0)}^{2}=\frac{\epsilon_{\text{NC}(0)}^{2}}{T^{2}}=\mathcal{O}(T^{2}) \tag{5}\]

Proof.: According to the definition, \(\sigma_{NC(0)}^{2}\) can be written as follows:

\[\sigma_{NC(0)}^{2} =\frac{\sum_{p}{(\mathbf{v}^{t^{*}}-\mathbf{v}^{0})^{2}}}{M}= \frac{\sum_{p}{((\mathbf{v}^{t^{*}}T-\mathbf{v}^{0}T)/T)^{2}}}{M} \tag{6}\] \[=\frac{\sum_{p}{((\mathbf{x}^{T}-\mathbf{x}^{0})-\mathbf{v}^{0}T )^{2}}}{MT^{2}}=\frac{\sum_{p}{(\int_{0}^{T}(\mathbf{v}(t)-\mathbf{v}^{0}T)dt) ^{2}}}{MT^{2}}, \tag{7}\]

where \(M\gg N\) is the number of all samples. The term \((\mathbf{v}(t)-\mathbf{v}^{0}T)\) is a typical (degree 0) polynomial interpolation error and can be defined as:

\[\epsilon_{\text{IP}(0)}(t) =\mathbf{v}(t)-\mathbf{v}^{0}T=\mathbf{v}(t)-P_{k}(t)\big{|}k=0 \tag{8}\] \[=\frac{\mathbf{v}^{(k+1)}(\xi)}{(k+1)!}(t-t_{0})(t-t_{1})...(t-t_{ k})\big{|}k=0=\mathbf{v}(\xi)t, \tag{9}\]

where \(0\leq\xi\leq T\), and the corresponding integration error is:

\[\epsilon_{\text{NC}}(0)=\int_{0}^{T}\epsilon_{\text{IP}(0)}(t)dt=\mathbf{v}( \xi)\int_{0}^{T}tdt=\frac{1}{2}\mathbf{v}(\xi)T^{2}=\mathcal{O}(T^{2}). \tag{10}\]

Therefore, we obtain the final variance:

\[\sigma_{NC(0)}^{2}=\frac{M\mathcal{O}(T^{4})}{MT^{2}}=\mathcal{O}(T^{2}), \tag{11}\]

concluding the proof.
### Nc (k)

In comparison with the \(0\) degree polynomial approximation NC\((0)\), the high order NC\((k)\) is more accurate and yields only a little extra computational cost.

Suppose that we now have \(K+1\) (usually \(K\leq 8\) due to the catastrophic Runge's phenomenon [40]) points \(\mathbf{x}_{i}^{0},\mathbf{x}_{i}^{1},...,\mathbf{x}_{k}^{K}\) equally spaced on time interval \([0,T]\), then the prediction for the future coordinate \(\hat{\mathbf{x}}_{i}^{T}\) in Equation (4) can be re-written as:

\[\hat{\mathbf{x}}_{i}^{T}=\mathbf{x}_{i}^{0}+\frac{T}{K}\sum_{k=0}^{K}w^{k} \hat{\mathbf{v}}_{i}^{k}, \tag{12}\]

where \(\{w^{k}\}\) is the coefficients of Newton-Cotes formulas and \(\hat{\mathbf{v}}_{i}^{k}\) denotes the predicted velocity at time \(t^{k}\). To obtain \(\hat{\mathbf{v}}_{i}^{k},k\geq 1\), we simply regard EGNN as a recurrent model [41, 42] and re-input the last output velocity and coordinate. Some popular sequence models like LSTM [43] or Transformer [13] may be also leveraged, which we leave to future work. Then, the equation of updating the predicted velocity \(\hat{\mathbf{v}}_{i}^{k}\) at \(t^{k}\) can be written as:

\[\hat{\mathbf{v}}_{i}^{k}=\phi_{v}(\mathbf{h}_{i})\mathbf{v}_{i}^{k-1}+\frac{ \sum_{j\neq i}(\mathbf{x}_{i}^{k-1}-\mathbf{x}_{j}^{k-1})\mathbf{m}_{ij}}{N-1}, \tag{13}\]

where \(\mathbf{v}_{i}^{k-1},\mathbf{x}_{i}^{k-1}\) denote the velocity and coordinate at \(t^{k-1}\), respectively. Note that, Equation (13) is different from the form used in a multi-layer EGNN. The latter always uses the same input velocity and coordinate as constant features in different layers. By contrast, we first predict the next velocity and coordinate and then use them as input to get the new trainable velocity and coordinate. This process is more like training a language model [13, 41, 44, 43, 45].

From the interpolation theorem [26, 46], there exists a unique polynomial \(\sum_{k=0}^{K}C^{k}t^{(k)}\) of degree at most \(K\) that interpolates the \(K+1\) points, where \(C^{k}\) is the coefficients of the polynomial. To avoid ambiguity, we here use \(t^{(k)}\) to denote \(t\) raised to the power of \(k\), whereas \(t^{k}\) to denote the \(k\)-th time point. \(\sum_{k=0}^{K}C^{k}t^{(k)}\) thus can be constructed by the following equation:

\[\sum_{k}^{K}C^{k}t^{(k)}=[\mathbf{v}(t^{0})]+[\mathbf{v}(t^{0}),\mathbf{v}(t^ {1})](t-t^{0})+...+[\mathbf{v}(t^{0}),...,\mathbf{v}(t^{K})](t-t^{0})(t-t^{1} )...(t-t^{K}), \tag{14}\]

where \([\cdot]\) denotes the divided difference. In our case, the \(K+1\) points are equally spaced over the interval \([0,T]\). According to Newton-Cotes formulas, we will have:

\[\int_{0}^{T}\mathbf{v}(t)dt\approx\int_{0}^{T}\sum_{k=0}^{K}C_{k}t^{(k)}dt= \frac{T}{K}\sum_{k=0}^{K}w^{k}\mathbf{v}^{k}. \tag{15}\]

Particularly, \(\{w^{k}\}\) are Newton-Cotes coefficients that are irrelevant to \(\mathbf{v}(t)\). With additional observable points, the estimation for the integration can be much more accurate. The objective thus can be written as follows:

\[\operatorname*{argmin}_{\{\hat{\mathbf{v}}^{k}\}}\sum_{p}\mathbf{x}^{T}-\hat {\mathbf{x}}^{T}=\sum_{p}\mathbf{x}^{T}-\mathbf{x}^{0}-\frac{T}{K}\sum_{k=0}^{ K}w^{k}\hat{\mathbf{v}}^{k}. \tag{16}\]

**Proposition 3.3** (Variance of \(\sigma_{NC(k)}^{2}\)).: \(\forall k\in\{0,1,...\},\epsilon_{NC(k)}\geq\epsilon_{NC(k+1)}\)_, and consequently:_

\[\sigma_{NC(0)}^{2}\geq\sigma_{NC(1)}^{2}\geq...\geq\sigma_{NC(k)}^{2}. \tag{17}\]

Proof.: Please see Appendix B.2 for details. Briefly, we only need to prove that \(\epsilon_{\text{NC}(0)}\geq\epsilon_{\text{NC}(1)}\), where NC (1) is associated with Trapezoidal rule. 

**Proposition 3.4** (Equivariance of NC).: _NC possesses the equivariance property if the backbone model \(\mathcal{M}\) (e.g., EGNN) in NC possesses this property._

Proof.: Please see Appendix B.3.
### NC (k) and NC\({}^{+}\) (k)

In the previous formulation, in addition to the initial and terminal points, we also need the other \(K-1\) points for supervision which yields a regularization loss:

\[\mathcal{L}_{r}=\sum_{p}\sum_{k=0}^{K}\|\mathbf{v}^{k}-\hat{\mathbf{v}}^{k}\|. \tag{18}\]

We denote NC (k) with the above loss by NC\({}^{+}\) (k). Minimizing \(\mathcal{L}_{r}\) is equivalent to minimizing the difference between the true velocity and predicted velocity at each time point \(t^{k}\) for each particle. \(\|\cdot\|\) can be arbitrary distance measure, e.g., L2 distance. Although the \(K-1\) points of data can be easily accessed in practice, we argue that a loose version (without \(\mathcal{L}_{r}\)) NC (k) can also learn promising estimation of the integration \(\int_{0}^{T}\mathbf{v}(t)dt\). Specifically, we still use \(K+1\) points \(\{\hat{\mathbf{v}}^{k}\}\) to calculate the integration over \([0,T]\), except that the intermediate \(K-1\) points are unbounded. The model itself needs to determine the proper values of \(\{\hat{\mathbf{v}}^{k}\}\) to optimize the same objective defined in Equation 16. In Section 4.5, we design an experiment to investigate the difference between the true velocities \(\{\mathbf{v}^{k}\}\) and the predicted velocities \(\{\hat{\mathbf{v}}^{k}\}\).

### Computational Complexity and Limitations

In comparison with EGNN, NC (k) does not involve additional neural layers or parameters. Intuitively, we recurrently use the last output of EGNN to produce new predicted velocities at next time point, the corresponding computational cost should be \(K\) times that of the original EGNN. However, our experiment in Section 4.4 shows that the training time did not actually linearly increase w.r.t. \(K\). One reason may be the gradient in back-propagation is still computed only once rather than \(K\) times.

We present an implementation of NC\({}^{+}\) in Algorithm 1. We first initialize all variables in the model and then recurrently input the last output to the backbone model to collect a series of predicted velocities. We then calculate the final predicted integration and plus the initial coordinate as the final output (i.e., Equation (12)). Finally we compute the losses and update the model by back-propagation.

```
1:Input: the dataset \(\mathcal{D}\), the backbone model \(\mathcal{M}\), number of steps \(k\) for NC;
2:repeat
3:for each batched data in the training set \((\{\boldsymbol{X}^{0},...,\boldsymbol{X}^{k}\},\{\boldsymbol{V}^{0},..., \boldsymbol{V}^{k}\})\)do
4:\(\hat{\boldsymbol{X}}^{0}\leftarrow\boldsymbol{X}^{0},\quad\hat{\boldsymbol{V}} ^{0}\leftarrow\boldsymbol{V}^{0}\);
5:for\(i:=1\)to\(k\)do
6:\(\hat{\boldsymbol{X}}^{i},\hat{\boldsymbol{V}}^{i}\leftarrow\mathcal{M}(\hat{ \boldsymbol{X}}^{i-1},\hat{\boldsymbol{V}}^{i-1})\);
7:endfor
8: Compute the final prediction \(\hat{\boldsymbol{X}}^{k}\) by Equation 12;
9: Compute the main prediction loss \(\mathcal{L}_{main}\) according to Equation (16);
10: Compute the velocity regularization loss \(\mathcal{L}_{r}\) according to Equation (18);
11: Minimize \(\mathcal{L}_{main}\), \(\mathcal{L}_{r}\);
12:endfor
13:until the loss on the validation set converges.
```

**Algorithm 1** Newton-Cotes Graph Neural Network

## 4 Experiment

We conducted experiments on three different tasks towards reasoning system dynamics. The source code and datasets have been uploaded and will be available on GitHub.

### Settings

We selected three state-of-the-art methods as our backbone models: RF [17], EGNN [18], and GMN [19]. To ensure a fair comparison, we followed their best parameter settings (e.g., hidden size and the number of layers) to construct NC. In other words, the main settings of NC and the baselines
were identical. The results reported in the main experiments were based on a model of NC (2) for a balance of training cost and performance.

We reported the performance of NC w/ additional points (denoted by \(\mathbf{NC}^{+}\)) and w/o additional points (the relaxed version, denoted by \(\mathbf{NC}\)).

### Datasets

We leveraged three datasets towards different aspects of reasoning dynamics to exhaust the performance of NC. We used the N-body simulation benchmark proposed in [19]. Compared with the previous ones, this benchmark has more different settings, e.g., number of particles and training data. Furthermore, it also considers some physical constraints (e.g., hinges and sticks) widely existing in the real world. For molecular dynamics, We used MD17 [28] that consists of the molecular dynamics simulations of eight different molecules as a scenario. We followed [19] to construct the training/testing sets. We also considered reasoning human motion. Follow [19], the dataset was constructed based on \(23\) trials in CMU Motion Capture Database [27].

As NC\({}^{+}\) requires additional supervised data, we accordingly extracted more intermediate data points between the input and target from each dataset. It is worth noting that this operation was tractable as the sampling frequency of the datasets was much higher than our requirement.

### Main Results

We reported the main experimental results w.r.t. three datasets in Tables 1, 2, and 3, respectively. Overall, NC\({}^{+}\) significantly improved the performance of all baselines on almost all datasets and across all settings, which empirically demonstrates the generality as well as the effectiveness of the proposed method. Remarkably, NC (without extra data) achieved slightly low results but still had great advantages over performance in comparison with three baselines.

Specifically, The average improvement on N-body datasets was most significant, where we observe a near 20% decrease in prediction error for all three baseline methods. For example, NC (RF) and NC\({}^{+}\) (RF) greatly outperformed the original RF and even had better results than the original EGNN under some settings (e.g., 3,2,1). This phenomenon also happened on the molecular dynamics datasets,

\begin{table}
\begin{tabular}{l r r r r r r r r r r} \hline \hline  & \multicolumn{3}{c}{\(|\begin{subarray}{c}\text{Train}\\ 0.1\end{subarray}=500\end{subarray}\)} & \multicolumn{3}{c}{\(\begin{subarray}{c}\text{Train}\\ 0.1\end{subarray}=1500\end{subarray}\)} \\  & 1,2,0 & 2,1 & 0,10 & 5,3,3 & 1 & 2,0 & 2,1 & 0,10 & 5,3,3 \\ \hline Linear & 8.23\(\pm\)0.75\(\pm\)0.09 & 7.55\(\pm\)0.09 & 7.96\(\pm\)0.10 & 11.36\(\pm\)0.00 & 11.62\(\pm\)0.00 & 8.22\(\pm\)0.07 & 7.55\(\pm\)0.09 & 9.76\(\pm\)0.00 & 11.36\(\pm\)0.00 & 11.62\(\pm\)0.00 \\ GNN [9] & 5.33\(\pm\)0.07\(\pm\)0.01\(\pm\)0.08 & 7.58\(\pm\)0.08 & 9.83\(\pm\)0.04 & 9.77\(\pm\)0.23 & 3.61\(\pm\)0.13 & 3.23\(\pm\)0.34 & 3.42\(\pm\)0.11 & 7.97\(\pm\)0.04 & 7.91\(\pm\)0.31 \\ TFN [39] & 11.54\(\pm\)0.08 & 9.87\(\pm\)0.27 & 11.66\(\pm\)0.08 & 13.43\(\pm\)0.31 & 12.23\(\pm\)0.11 & 5.86\(\pm\)0.34 & 4.97\(\pm\)0.23 & 8.71\(\pm\)0.14 & 11.21\(\pm\)0.12 & 10.75\(\pm\)0.08 \\ SE(3)-Tr. [29] & 5.54\(\pm\)0.08 & 5.14\(\pm\)0.08 & 8.95\(\pm\)0.04 & 11.42\(\pm\)0.00 & 11.59\(\pm\)0.00 & 5.02\(\pm\)0.04 & 4.68\(\pm\)0.08 & 8.39\(\pm\)0.02 & 10.82\(\pm\)0.03 & 10.85\(\pm\)0.02 \\ \hline RF [17] & 3.50\(\pm\)0.17 & 3.07\(\pm\)0.24 & 5.25\(\pm\)0.44 & 7.59\(\pm\)0.25 & 7.73\(\pm\)0.39 & 2.97\(\pm\)0.15 & 2.19\(\pm\)0.11 & 3.80\(\pm\)0.25 & 5.71\(\pm\)0.31 & 5.66\(\pm\)0.27 \\ NC (RF) & **2.84\(\pm\)0.23**\(\pm\)**0.04 & **4.32\(\pm\)**0.08** & **6.67\(\pm\)0.28** & 7.14\(\pm\)0.24 & 2.91\(\pm\)0.11 & **1.76\(\pm\)0.01** & **3.23\(\pm\)0.01** & 5.09\(\pm\)0.05 & 5.34\(\pm\)0.03 \\ NC\({}^{+}\) (RF) & 2.92\(\pm\)0.24 & **3.24\(\pm\)**0.08** & **4.28\(\pm\)**0.05** & 7.02\(\pm\)**0.14** & 7.46\(\pm\)0.23** & **2.87\(\pm\)0.18** & 1.76\(\pm\)0.01 & 3.24\(\pm\)0.02 & **5.06\(\pm\)0.08** & **5.23\(\pm\)0.29** \\ \hline EGNN [18] & 2.81\(\pm\)0.12 & 2.27\(\pm\)0.04 & 4.67\(\pm\)0.47 & 4.75\(\pm\)0.08 & 4.59\(\pm\)0.29 & 2.59\(\pm\)0.10 & 1.86\(\pm\)0.01 & 2.54\(\pm\)0.01 & 2.79\(\pm\)0.34 & 3.25\(\pm\)0.07 \\ NC (EGNN) & 2.41\(\pm\)0.03 & 2.18\(\pm\)0.02 & 3.53\(\pm\)0.01 & 4.26\(\pm\)0.03 & 4.13\(\pm\)0.23 & **2.23\(\pm\)0.13** & 1.91\(\pm\)0.03 & 2.14\(\pm\)0.03 & 2.36\(\pm\)0.05 & **2.86\(\pm\)0.03** \\ NC\({}^{+}\) (EGNN) & **2.25\(\pm\)0.04** & **2.07\(\pm\)**0.06** & **2.01\(\pm\)0.02** & **3.54\(\pm\)0.06** & **3.96\(\pm\)0.04** & 2.32\(\pm\)0.05 & **1.86\(\pm\)0.13** & 2.39\(\pm\)0.01 & **2.35\(\pm\)0.07** & 2.99\(\pm\)0.02 \\ \hline GMN [19] & 1.84\(\pm\)0.02 & 2.02\(\pm\)0.02 & 2.48\(\pm\)0.04 & 2.92\(\pm\)0.04 & 4.08\(\pm\)0.03 & 1.68\(\pm\)0.04 & 1.47\(\pm\)0.03 & 2.10\(\pm\)0.04 & 2.32\(\pm\)0.02 & 2.86\(\pm\)0.01 \\ NC (GMN) & **1.55\(\pm\)0.07** & 1.58\(\pm\)0.02 & 2.07\(\pm\)0.03 & 2.73\(\pm\)0.02 & 2.99\(\pm\)0.03 & 1.59\(\pm\)0.01 & 1.18\(\pm\)0.05 & 1.47\(\pm\)0.01 & **1.66\(\pm\)0.01** & 1.93\(\pm\)0.03 \\ NC (GMN) & 1.57\(\pm\)0.02 & **1.43\(\pm\)0.02** & **2.03\(\pm\)0.04** & **2.57\(\pm\)0.04** & **2.72\(\pm\)0.03** & **1.49\(\pm\)0.12** & **1.17\(\pm\)0.04** & **1.44\(\pm\)0.01** & 1.70\(\pm\)0.06** & **1.91\(\pm\)0.02** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Prediction error (\(\times 10^{-2}\)) on N-body dataset. The header of each column “\(p,s,h\)” denotes the scenario with \(p\) isolated particles, \(s\) sticks and \(h\) hinges. Results averaged across 3 runs.

\begin{table}
\begin{tabular}{l r r r r r r r r r r r} \hline \hline  & \multicolumn{3}{c}{\(\begin{subarray}{c}\text{Aquin}\\ \text{RF}\end{subarray}\)} & \multicolumn{3}{c}{\(\begin{subarray}{c}\text{Run}\\ 0.1\end{subarray}=1500\end{subarray}\)} & \multicolumn{3}{c}{\(\begin{subarray}{c}\text{Mian}\\ 0.1\end{subarray}=1500\)} & \multicolumn{3}{c}{\(\begin{subarray}{c}\text{Train}\\ 0.2\end{subarray}=1500\)} \\  & 1,2,0 & 2,1 & 0,10 & 5,3,3 & 1 & 2,0 & 2,1 & 0,10 & 3,2,1 & 0,10 & 5,3,3 \\ \hline Linear &
where we find that NC\({}^{+}\) (EGNN) not only defeated its original version, but also the original GMN and even NC\({}^{+}\) (GMN) on many settings.

### Impact of estimation step \(k\)

We conducted experiments to explore how the performance changes with respect to the estimation step \(k\). The experimental results are shown in Figure 2.

For all three baselines, the performance improved rapidly as the growth of step \(k\) from \(0\) to \(2\), but this advantage gradually vanished as we set larger \(k\). Interestingly, the curve of training time showed an inverse trend, where we find that the total training time of NC (5) with GMN was almost two times slower than the NC (1) version. Therefore, we set \(k=2\) in previous experiments to get the best balance between performance and training cost.

By comparing the prediction errors of different baselines, we find that GMN that considers the physical constraints (e.g., sticks and hinges) usually had better performance than the other two methods, but the gap significantly narrowed with the increase of step \(k\). For example, on MD17 and motion datasets, EGNN obtained similar prediction errors for \(k\geq 2\) while demanding less training time. This phenomenon demonstrates that predicting the integration with multiple estimations lowers the requirement for model capacity.

Overall, learning with NC significantly reduced the prediction errors of different models. The additional training time was also acceptable in most cases.

### A Comparison between NC and NC\({}^{+}\)

The only difference between NC and NC\({}^{+}\) was the use of the regularization loss in NC\({}^{+}\) to learn the intermediate velocities. Therefore, we designed experiments to compare the intermediate results and the final predictions of NC and NC\({}^{+}\) with GMN as the backbone model.

We first tracked the average intermediate velocity prediction errors on valid datasets in terms of training epochs. The results are shown in Figure 3, from which we observe that the prediction errors were not huge for both two methods, especially NC\({}^{+}\) obtained highly accurate prediction for the intermediate velocities with different \(k\). For NC, we actually find that it implicitly learned to model

\begin{table}
\begin{tabular}{l c c c c c} \hline \hline  & GNN & TFN & SE(3)-Tr. & RF & EGNN & GMN \\ \hline Raw & 67.3\({}_{\pm 1}\) & 66.9\({}_{\pm 2}\) & 60.9\({}_{\pm 0.9}\) & 197.0\({}_{\pm 10}\) & 59.1\({}_{\pm 2}\) & 43.9\({}_{\pm 1}\) \\ NC & N/A & N/A & N/A & 164.8\({}_{\pm 1}\) & 55.8\({}_{\pm 6}\) & 30.0\({}_{\pm 4}\) \\ NC\({}^{+}\) & N/A & N/A & N/A & **162.6\({}_{\pm 18}\)** & **51.3\({}_{\pm 4}\)** & **29.6\({}_{\pm 8}\)** \\ \hline \hline \end{tabular}
\end{table}
Table 3: Prediction error (\(\times 10^{-2}\)) on motion capture. Results averaged across 3 runs.

Figure 2: A detailed comparison of NCGNN without extra data on \(5\) datasets, w.r.t. \(k\), average of \(3\) runs, conducted on a V100. The lines with dot, circle, and triangle denote the results of GMN, EGNN, and RF, respectively. The first, second, and third rows are the results of prediction error, training time of one epoch, and total training time, respectively.
the intermediate velocities more or less, even though its prediction errors gradually improved or dramatically fluctuated during training. This observation empirically demonstrates the effectiveness of the intermediate velocities in predicting the final state of the system. We also reported the results on MD17 and N-body datasets in Appendix C.

We then visualized the intermediate velocities and the final predictions in Figure 4. The intermediate velocities were visualized by adding the velocity with the corresponding coordinate at the intermediate time point. From Figure 4, we can find that NC\({}^{+}\) (green ones) learned very accurate predictions with only a few non-overlapped nodes to the target nodes (red ones). For the loose version NC, it learned good predictions for the backbone part, but failed on the limbs (especially the legs). The motion of the main backbone was usually stable, while that of the limbs involved swing and twisting. We also conducted experiments on the other two types of datasets, please see Figure 8 in Appendix C.

### On the Time Evolution of Dynamic Systems

In previous experiments, we show that NC\({}^{+}\) was capable of learning highly accurate predictions of intermediate velocities. This inspired us to realize the full potential of NC\({}^{+}\) by producing a series of predictions given only the initial state as input. For each time point (including the intermediate time points), we used the most recent \(k+1\) points of averaged data as input to predict the next velocity and coordinate. This strategy produced much more stable predictions than the baseline NC (0) - directly using the last output as input to predict the next state.

The results are shown in Figure 5, from which we can find that both the baseline (blue lines) and NC\({}^{+}\) (green lines) had good predictions at the initial steps. However, due to the increasing acculturated errors, the baseline soon collapsed and only the main backbone was roughly fit to the target. By contrast, NC\({}^{+}\) was still able to produce promising results. The deviations mainly concentrated on the limbs. Therefore, we argue that NC\({}^{+}\) can produce not only the highly accurate one-step prediction, but also the relatively reliable consecutive predictions. We also uploaded the.gif files of examples on all three datasets in Supplementary Materials, which demonstrated the consistent conclusion.

Figure 4: Visualization of the intermediate velocities w.r.t. \(k\). The red, blue, and green lines denote the target, prediction of NC, and prediction of NC\({}^{+}\), respectively.

Figure 5: Consecutive predictions on the Motion dataset. The red, blue, and green lines denote the target, the prediction of NC (0), and the prediction of NC\({}^{+}\) (2), respectively. The interval between two figures is identical to the setting in the main experiment.

Figure 3: The prediction errors of intermediate velocities on valid set, w.r.t. training epoch. The blue and green lines denote the prediction errors of NC and NC\({}^{+}\), respectively.
Conclusion

In this paper, we propose NC to tackle the problem of reasoning system dynamics. We show that the existing state-of-the-art methods can be derived to the basic form of NC, and theoretically/empirically prove the effectiveness of high order NC. Furthermore, with a few additional data provided, NC\({}^{+}\) is capable of producing promising consecutive predictions. In future, we plan to study how to improve the ability of NC\({}^{+}\) on the consecutive prediction task.
