Newton–Cotes Graph Neural Networks: On the Time Evolution of Dynamic Systems
Authors: Lingbing Guo
Weiqing Wang
Zhuo Chen
Ningyu Zhang
Zequn Sun
Yixuan Lai
Qiang Zhang
Huajun Chen
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Equivariant Graph Neural Networks
Molecular Dynamics
N-body System
Human Motion

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.

Newton-Cotes Graph Neural Networks: On the Time Evolution of Dynamic Systems
Linghing Guo 12,32 Weiqing Wang'; Zhuo Chen 1.2.3 Ningyu Zhang 1.2 Zequn Sun Yixuan Lai Qiang Zhang and Huajun Chen
'College of Computer Science and Technology, Zhejiang University "Zhejiang University Ant Group Joint Laboratory Knowledge Graph $ZJU- ZZaU-Hangzhou Global Scientific and Technological Innovation Center
'Department of Data Science & AI, Monash University 5 "State Key Laboratory for Novel Software Technology, Nanjing University
Abstract
Reasoning system dynamics is of the most important analytical approaches for many scientific studies. With the initial state system as input, the recent graph neural networks (baee 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 new approach to predict the integration based 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 ofee-----herrr art methods.
1 Introduction
Reasoning the time evolution of dynamic systems has been a long-ter challenge for hundreds of years [1]]]. Despite that the advances in computer manufacturing nowadays make it possible to simulate long trajectories of complicated systems constrained by multiple foree laws, the computational cost still imposes a heavy burden on the scientific communities 3
Recent graph neural networks (GNNs) 2 -bssed methods provide alternative solution to predict the future states directly with only the initial state as input [0420]. Take molecular dynamics example, the atoms in molecule can be regarded nodes with different
(MD) 4 24 labels, and thus can be encoded by GNN. As the input and the target are atomic coordinates, the learning problem becomes 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 x distant to the given initial coordinate x", with x" and the initial velocity v as input. Then, the
"Equal Contribution
Corresponneecce to: (qiang.zhang, es, huajunsir) zju..ed.cen.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

predicted coordinate x" can be written as:
x-= + T...
(1)
For complicated system comprising multiple particles, predicting the velocity term v rather than the coordinate x" improves the robustness and performance. Specifically, by subtracting the initial coordinate x", the learning target can be regarded as the normalized future coordinate, i.e., &
(x" x (T-), which is why all state-of- off-the-rr art methods adopt this strategy [8 19] Nevertheless, the current strategy still has significant drawback from the view of numerical integration.
As illustrated in Figure T supposed that the actual velocity of particle is deseribed by the curve v()). To predict the future state x"" at 1 = T, the existing methods adopt a constant
The existing methods
The two-step estimation
v(t)
The true integration over [0,
estimation approach, i.e., x = xx + 8 = x + So dt. Evidently, the prediction error could be significant as it purely relies the fitness of neural models. If we alternatively choose a simple two-step estimation (i.c., Trapezoidal rule (251), the prediction error may be reduced only the blue area. In other words, the model only needs to compensate for the blue area.
2
T Figure 1: Illustration of different estimations. The rectangle, trape- zoid, and striped areas denote the basic (used in the existing works), two-step, and true estimations, re- spectively. 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.
0
In this paper, propose Newton-Cotes graph neural networks (abbr. NC) to estimate multiple velocities at different time points and compute the integration with Newton-Cotes for- mulas, Newton-Cotes formulas are 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-Coke formulas is that the integration is computed by aggregating the values of vtt) at these points with a group of weights {w",... J. where k refers to the order of Newton- Cotes formulas and (w ...., w²) } irrelevant to the integrand ctt) and can be pre-computed by Lagrange interpolating polynomial 26
We show that the existing works theoretically prove that the prediction error of this estimation as well as the learning difficulty will continually reduce the increase of estimation step k.
be naturally derived to the basic version NC (k 0) and
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 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/29), with state-of-the-art methods as baselines [17-19. The results show that NC improves all the baseline methods with 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 Dynamies 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 hierarchical graph convolution auto-encoder (33 to encode the objects, or leveraging ordinary differential equations for energy
99] I0 32) or
2

---

conservation [77. 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 Recent methods considering physical symmetry and Eu- clidean equivariance have achieved state-of-the-art performance in modeling system dynamics [14] Specifically, [1416] force the translation equivariance, but the rotation equivariance is overlooked. TFN [39] further leverages spherical filters to achieve rotation equivariance. SE(3) Trans- former 29 proposes use the Transformer
19
38
[3] architecture to model 3D cloud points directly in Euclidean space. EGNN [88 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 -fate-of-thh-art performance while maintaining the same level of computational cost. SEGNN 3X 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 [8] and GMN (19) in this sub-area usually regard the velocity as an important feature. These methods can be set as the backbone model NC, and the models themselves belong to the basic NC (0).
Neural Ordinary Differential Equations The neural ordinary differential equation (neural ODE) methods 41 parameterize ODEs with neural layers and solve them by numerical solvers such as 4th order Hange-Kutta, These solvers are originally designed for numerical integration. Our work draws inspiration from numerical integration algorithms, where the integrand is known only at certain points. The goal is to efficiently approximate the integral to desired precision. Neural ODE methods repeatedly perform the numerical algorithms on small internals to obtain the numerical solution. Therefore, it may be feasible to iteratively perform our method to solve problems in neural ODEs.
3 Methodology
We start from preliminaries and then take the standard EGNN 18]
as example 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.
3.1 Preliminaries
We suppose that system comprises N particles P1, Pi.P..... PN. with the velocities V1, V V and the coordinates X1. X2, ...; XN 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 to process them and do not discuss the details in this paper.
Reasoning Dynamics Reasoning system dynamics is a classical task with a very long history M. One basic tool for reasoning or simulating system is numerical integration. Given the dynamics (e.g.. Langevin dynamics, well-used MD) and initial information, the interaction force and acceleration for each particle in the system can be estimated. However, the force is closely related the real-time coordinate, which means that one must re-calculate the system states for every very small time step (i.e., dt) 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- second), whereas the total simulation time can be 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.
3.2 EGNN
Although the existing EGNN methods have great advantage in computational cost over the conven- tional simulation tools, the potential prediction error that a neural model needs to compensate for is also huge. Take the standard EGNN [8] as an example, the equivariant graph convolution can be
written follows:
0 m = de (h, h, ||x - x911 eig),
(2)
3

---

where the aggregation function de takes three types of information as input: the feature embeddings h; and h, for the imput particle p; and an arbitrary particle Pj. respectively; the Euclidean distance
||x 22 at time O; and the edge attribute eij. Then, the predicted velocity and future coordinate can be calculated by the following equation:
v.. = (h.)
Σ j
(x] x.)mbi
(3)
N
x = x V.T,
(4)
where 0n RDD 7 R. a learnable function. From the above equations, we can find that the predicted velocity V is also correlated with three types of information: 1) the initial velocity v as a basis; 2) the pair-wise Euclidean distance ||x] - x311? (i.e., m.;) to determine the amount of force (analogous to Coulomb's 's law); and the pair-wise relative coordinate difference (x? - -xy) to determine the direction. For multi-layer EGNN and other EGNN models, v is still determined by these three types of information, which we refer interested readers to Appendix for details.
Proposition 3.1 (Linear mapping). The existing methods learn a linear mapping from the initial velocity v to the average velocity v over the interval [0,
Proof. Please see Appendix B.I
As the model makes prediction only based on the state and velocity at = O, we assume that v² fluctuates around v and follows normal distribution denoted by NNc((o) (v", 1y Nc(o)). Then, variance term (NC()) 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 NNC(o)). Assume that the average velocity v over [0, T" follows normal distribution with μ vo, then the variance the distribution is:
)2 = =
(5)
Proof. According to the definition, xxc(0) can be written as follows:
Ep (v' - voj2 -VVTT-/T)
2 NC(())
(6)
M
M
Ep ((x" - -x²) - oTyz
Ep (o (v(1) - v
MT²
MT2
where M is the number of all samples. The term (v(t) a typical (degree 0) polynomial interpolation error and can be defined as:
(IP( (t) v(t) vT = v()) - Pa(t)]
(8)
(()
(t (t) -t -t)((t)))...( (t - k= v(E)t,
(k + 1)!
where 0 < 5 T, and the corresponding integration error is:
(xC(C) = L." (PP vv(E)
v(E) 1' tdt =
(10)
Therefore, obtain the final variance:
MO(T4) MT = O
NC
(11)
concluding the proof.

---

3.3
NC (k)
In comparison with the degree polynomial approximation NC(0), the high order NC(k) is more accurate and yields only a little extra computational cost.
Suppose that points x,, x.,... x equally spaced on time interval [0, T). then the prediction for the future coordi- nate x in Equation (4) can re-written as:
have K (usually K s 8 due to the catastrophic Runge's phenomenon 42
K
T
v w
(12)
K
k=0
where {w"} is the coefficients of Newton-Coud formulas and v denotes the predicted velocity at time t To obtain vf, k 2 1, we simply regard EGNN as a recurrent model [43 44 and re-input the last output velocity and coordinate. Some popular sequence models like LSTM 45 or Transformer 13 may be also leveraged, which we leave to future work. Then, the equation updating the predicted velocity v at th can be written as:
Eniss (x
x
m
v = 2h (h:lvv
(13)
- N
where vi-1, x denote the velocity and coordinate at th-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 constant features 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 language model
45
From the interpolation theorem 226 48 there exists unique polynomial Σx=o Ctky(h)) of degree at most K that interpolates the K points, where C* the coefficients of the polynomial. To avoid ambiguity, we here use {(k) to denote t raised the power of k, whereas & to denote the k-th time point. Σk=o (k) thus can be constructed by the following equation:
K
=(())==++( +v() +... (v("), v(th M(t- ( - (t....tt tk),
=
(14)
where [-] denotes the divided difference. In case, the K points equally spaced over the interval (0, 7). According to Newton-Cotes formulas, we will have:
[ v(t)tt se J K
E
=
T
K
T
Σ
K
(15)
k=
Particularly, {w*} are Newton- Newton-Cotes coefficients that are irrelevant to v(I). With additional observ- able points, the estimation for the integration can be much more accurate. The objective thus be written as follows:
K
argmin Σ x" x2 (
-xx'
x2
T
w
(16)
K
k=0
Proposition 3.3 (Variance of o
)ow)) ). wk (0, (.,..,,,NCCCL
CNC(++1) and consequently:
2x0 2 2oou 2 oxcc))
(17)
Proof. Please see Appendix B.2 for details. Briefly, we only need to prove that eNC(O) 2 <(C(1))
where NC (1) is associated with Trapezoidal rule.
Proposition 3.4 (Equivariance of NC). NC possesses the equivariance property if the backbone model M (e.g., EGNN) in NC possesses this property.
Proof. Please see Appendix B.3
5

---

Algorithm Newton-Cotes Graph Neural Network
1: Input: the dataset D, the backbone model M, number of steps & for NC;
2: repeat
3: for each batched data in the training set (fx), ...., X*), (Vo V*)) do
4:
X X 0 V
5: for i = 1 to k do
6:
X V M - 1.-1V
7:
end for
8:
Compute the final prediction X by Equation 12
9:
Compute the main prediction loss Lmain according to Equation (16] Compute the velocity regularization loss L, according to Equation 18
10:
11: Minimize Lmain by:
12:
end for
13:
until the loss the validation set converges.
3.4
NC (k) and NC (k)
In the previous formulation, in addition to the initial and terminal points, we also need the other
K points for supervision which yields regularization loss:
C,
4= P
(18)
denote NC (k) with the above loss by NC+ (k). Minimizing C, is equivalent to minimizing the difference between the true velocity and predicted velocity at each time point ( for each particle. ||- can be arbitrary distance measure, e.g., distance. Although the K
points of data can be easily
argue that a loose version (without C,) NC (k) can also learn promising estimation of the integration Io' (f)) Specifically, we still use K points {v*) to calculate the integration over [0, 7], except that the intermediate K points are unbounded. The model itself needs to determine the proper values of (V*} optimize the same objective defined in Equation 16 In Section 4.5 design experiment to investigate the difference between the true velocities {v²] and the predicted velocities {
accessed in practice,
3.5 Computational Complexity and Limitations
In comparison with EGNN, NC (k) does not involve additional neural layers 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 shows that the training time did not actually linearly increase w.r.t. K. One reason may be the gradient in back propaaaaann still computed only once rather than K times.
We present an implementation of + in Algorithm| I and then recurrently input the last output 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 (2)). Finally compute the losses and update the model by Dack-propagation.
We first initialize all variables in the model
Experiment
We conducted experiments three different tasks towards reasoning system dynamics. The source code and datasets are available at https://jitubbcoo/m/uu///NCCCC
4.1 Settings
We selected three state-of-the-art methods our backbone models: RF 17 EGNN 18 and
GMN [19]. To ensure 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

---

Table 1:
Prediction error (x10-2) 3) on N-body dataset. The header of each column "p, 8, h" denotes
the scenario with p isolated particles, sticks and h hinges. Results averaged across 3 runs.
Train] 500
Train| 1500
1,2,0
2,0,1
3,2,1
0,10,0
5,3,3
1,2,0
2,0,1
3,2,1
0,10,0
5,3,3
Lincar
7.55
9.761000
11.36
11.62=000
8
2.02
7.55
9.76
11.36.0m
11.62:0m
GNNN
0.33 2.0
7.58100s
9.83
9.77
3.61
3,23
7.91 10: 7.97204
TFN
11 9.87
11.6620o6
13.43
12.23
5.86
4.97
8,51
11.21
10.75
SB(3) Tr.
20
5.54 5.14
8.95.00
11.42
11.59=0m
5.02-....s
4.682005
8.39
10.82gam
10.85
RF 17
3.5
3.07++024
5.25++0.4
+.59
-.73
2.97 =.19
3.80
=.31 5.66+027
NC (RF)
:.84 2.35
4.32::00
6.67
7.14=025
2.91 1.76:
3.23
5.09:mms
5.34
NC (RF)
1. 2.34
4.28::a.6
7.02+0.14
9.06-0.31
=.87
1.78=om
3.24-002
5,06=006
5.23+029
EGNN [8]
2.81 12
2.27
4.67.
4.75
4.59
2.55
1.86
2.54_nm
2.79
3.25
NC (EGNN)
2.41 2018
2.18=aa
3.53-20.01
4.26::00
4.13=00
2.23
1.91=003
2.14
2.36
2.86
NC* (EGNN)
2.25 2.07.am
2.0
3.54
3.96
2.32.10.6
1.86
2.39
2.35
2.99
GMN 19
:004 2.02#a.:
2.48-10.-
2.92
4.08=0a:
1.68
=.47:003
2.10=0044
2.32=am
2.86
NC (GMN)
gu5 1.58gam
2.07
2.73
2.99
1.59
1.18200s
1047=0m
1.66:um
1.93gam
NC (GMN)
1.57 1.43+a.2
2.03++0M4
2.57
:0:04
2.72
1.49:a4:
1.17=084 1.44-001
1.70=aam
+.91
Table 2: Prediction error (x10-3) on MD17 dataset. Results averaged across 3 runs.
Aspirin
Benzene
Ethanal
Malonnkdeydde
RF
EGNN
GMN
RF
EGNN
GMN
RF
EGNN
GMN
RF
EGNN
GMN
129
48.12:00 25.05
ca
10.87-mm
12.93.mm
12.8 12.81 :M:
12.97
NC
10.87 -M1
9.60 oaa
9.40 -ua
63.04 om
54.76.024
41.26--ax
4.62 mM1
4.62 -011
4.83 1M1
12.93 Am
12.92 PM
Naphthalene
Salicylic
Toloene
Uracil
RF
EGNN
GMN
RF
EGNN
GMN
RF
EGNN
GMN
RF
BONN
GMN
Row
0.50
0.47-
0.40am
1.23
1002au:
10.93aas 10.85
11.78aa
0.64.mm
0.64amm
1.09.ux 1.09 oA1
0.84.am
10.64.mm
10.17 10.14 AM
NC
0.39 .o1
0.37 om
0.37 Am
0.84 oAt
0.86 AA1
10.85 AM1
10.52 om
0.60 MA1
0.57 PM
0.56 IM:
were identical. The results reported in the main experiments were based on model of NC (2) for a balance of training cost and performance.
We reported the performance of NC w/ additional points (denoted by NC+) and w/o additional points (the relaxed version, denoted by NC).
4.2 Datasets
We leveraged three datasets towards different aspects of reasoning dynamics to exhaust the perfor- mance 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.
4.3 Main Results
We reported the main experimental results w.r.t. three datasets in Tables 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 well 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 bascline methods. For example, NC (RF) and NC
7

---

Table Prediction error (x10-2) on motion capture. Results averaged across runs.
GNN
TFN
SF(3)-T.
RF
EGNN
GMN
Raw
67. 1,32
66,9
60.9
197.0
19.1
43.9a11
N/A
N/A
164.8eLa
55.8461
30.0=14
NC
N/A
N/A
N/A
162.6:48
51.3+40
29.6=ua
153
(N6
O
143
123
028
3
NCJK)
NCVK)
NOKK
NDO0)
(9)
12
11.3
33
2.5
50
16.3
MCWK
NC03
NEV0
NOK0
NOW0
3000
4300
5009
3000
500
-100
+-GMN EGNN
-C-RF
1000
NC00)
NC03
(a) N-body (1, 2,
(b) N-body (5,
(c) MD17 (Aspirin)
(d) MD17 (Uracil)
(c) Motion Capture
Figure 2: A detailed comparison of NCGNN without extra data on 5 datasets, w.r.t. average of 3 runs, conducted 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.
V100, The lines with dot, circle, and triangle denote the results of GMN,
(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, where find that NC- (EGNN) not only defeated its original version, but also the original GMN and even NC+ (GMN) on many settings.
4.4 Impact of estimation step k:
We conducted experiments explore how the performance changes with respect to the estimation step k. The experimental results are shown Figure
For all three baselines, the performance improved rapidly as the growth of step k from 0 to 2, but this advantage gradually vanished 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 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 l 2 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.
4.5
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 compare the intermediate results and the final predictions of NC and NC+ with GMN the backbone model.

---

Difference
Differape
Differenee
Epoch 104
Epoch 194
Epoch 1903
(4) NC(2)
(b) NC(3)
(c) NC(5)
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.
NC(2) intermediate
(b) NC(3) intermediate
(c) NC(S) intermediate
NC(2) final
Figure 4: Visualization the intermediate velocities w.r.t. k. The red, blue, and green lines denote the target, prediction of NC, and prediction of NC respectivelly
t-0
[-]
t-2
t-3
t-5
t-6
1-7
[-8
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 NC+ (2), respectively. The interval between two figures is identical to the setting in the main experiment.
We first tracked the average intermediate velocity prediction errors on valid datasets in terms of training epochs. The results shown in Figure 3 from which 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. actually find that it implicitly learned model 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 17 and N-body datasets in Appendix C
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 the intermediate time point. From Figure [ we can find that NC (green ones) learned very accurate predictions with only 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 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 in Appendix C
4.6
On the Time Evolution of Dynamic Systems
In previous experiments, we show that NC+ was capable learning highly accurate predictions of intermediate velocities. This inspired us realize the full potential of NC by producing series of predictions given only the initial state as input. For each time point (including the intermediate time points), used the most recent &s + points of averaged data input to predict the next velocity and

---

Table Multi-step prediction results (x10-2) on the N-body and MD17 datasets. ADE and FDE
denote average displacement error and final displacement enror, respectively. The results of four
baseline methods are from [49].
N-body (Springs)
MD17 (Aspirin)
MD17 (Benzene)
MD17 (Ethanol)
MD17 (Malonaldehyde)
ADE
Accuracy
ADE
FDE
ADE
FDE
ADE
FDE
ADE
FDE
LSTM
45
53.5
17.59
24.79
6.06
9.46
7.73
9.88
15.14
22.90
NRI
93.0
12.60
18.50
1.89
2.58
6.69
8.78
12.79
19.86
GroupNet
51
10.62
14.00
2.02
2.95
6.00
7.88
7.99
12.49
EgMotion
49
16.8
97.6
5.95
8.38
1.18
1.73
5.05
7.02
5.85
9.02
NC (LgMotion)
16.2
98.9
4.74
6.76
1.15
1,63
5.02
6,87
5,19
8.07
coordinate. This strategy produced much more stable predictions than the baseline NC (0) direetly 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.
We also compared our method with the EGNN methods for multi-step prediction 52). Specifically, multi-step prediction takes sequence of states with fixed interval input and predicts sequence of future states with same length and interval. As recent multi-step prediction methods are still based on sequence encoder-decoder architectures, it would be interesting examine the applicability of NC to them. To this end, we adapted NC (k=2) to the source code the best-performing method EqMotion 49 and used identical settings for the hyper-parameters.
The experimental results are shown in Table 4 where we can observe that NC (EqMotion) out- performed the base model Eqtion all datasets and across all metrics. It is worth noting that the archiecture and parameter-setting of the EGNN module were identical to EqMotion and NC (Eqiotion), which clearly demonstrates the effectiveness of the proposed method. We have also released the source code of NC (EqMotion) our GitHub repository.
4.7 Ablation Study
We conducted ablation studies to investigate the effectiveness of the proposed Newton-Cotes in-
Table 5: Ablation study results (x10-2).
tegration. We implemented two alternative meth-
N-body (1,2.0)
MD17 (Aspirin)
Motion
ods: NC (coeffs= 1), where we set all the coef- ficients for summing the predicted velocities 1; and depth-varying ODE, where we employed common neural ODE solvers from 44 41 to
GMN
1.84
10.14
43.9
(cuellss= depth-varying ODE
1.82
10.08
42.2
NC
address our problem. We used an identical GMN as the backend EGNN model and conducted a simple hyper-paramede search for the learning rate and ODE solver.
The results are presented in Table 5 The two alternative methods methods generally performed better than the baseline GMN, The performance of depth-varying GDE was comparable to that of NC (coeffs= 1), but both were lower than that of the original NC. This observation further supports the effectiveness of method.
5 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 theoreticallyymmpirraalll 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 study how to improve the ability of NC+ on the consecutive prediction task.
10

---

Acknowledgement
We would like to thank all anonymous reviewers for their insightful and invaluable comments. This work is funded S
National Key R&D Program of China (Funding No.S022018FF00000O Zhejiang Provincial Natural Science Foundation of China (No.LGG22FF3001) and sponsored by CCF-Tencent Open Fund (CCF- RAGR20230122).
References
[I] Georges Louis Leclere comte de Buffon. Histoire naturelle, générale et particuliére: De la manière d' étuudier & de traiter l'histoire naturelle. Histoire & théorie de la terre. 1749, volume 1. L'Imprimerie Royale, 1774.
[2] Loup Verlet. Computer" experiments" on classical fluids. i. thermodynamical propertics of lennard-jones molecules. Physical review, 159(1):98, 1967.
[3] Glenn M Torrie and John P Valleau. Nonphysical sampling distributions in monte carlo free- energy estimation: Umbrella sampling. Journal of Computational Physics, 23(2):187-199, 1977.
[4] Daan Frenkel and Berend Smit. Understanding molecular simulation; from algorithms to applications, volume Elsevier, 2001.
[5] Frank Noé, Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science, 365(645)) 147, 2019.
[6] Ilya Vakser. Protein- Protein-protein docking: From interaction interactome. Biophysical journal, 10758):1785-1793, 2014.
[7] Han Wang, Linfeng Zhang, Jiequn Han, and E Weinan. Deepmd-kit: A deep learning package for many-body potential energy representation and molecular dynamies. Computer Physics Communications, 228:178: 178-184, 2018.
[8] Mieczyslaw Torchala, Jain H Moal, Raphael AG Chaleil, Juan Fernandez-Recio, and Paul A Bates. Swarmdock: a server for flexible protein-protein docking. Bioinformatics, 29(6): 807-809, 2013.
[9] Thomas N. Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. In ICLR, 2017.
[10] Petar Velickovic, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Liò, and Yoshua Bengio. Graph attention networks. In ICLR, 2018.
[11] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in neural information processing systems, 30, 2017.
[12] Jiale Liu and Xinqi Gong. Attention mechanism enhanced Istm with residual architecture and its application for protein-protein interaction residue pairs prediction. BMC bioinformatics, 20 (1):1-11,201 2019.
[13] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. In NeurIPS, 2017.
[14] Benjamin Ummenhofer. Lukas Prantl, Nils Thuerey, and Vladlen Koltun. Lagrangian fluid simulation with continuous convolutions. In ICLR, 2019.
[15] Alvaro Sanzhez-Gonzalez. Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning simulate complex physics with graph networks. In ICML, pages 8459-8468, PMLR, 2020.
[16] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gonzalez, and Peter W Battaglia. Learning mesh-based simulation with graph networks. arXiv preprint 2rXiv:2010.03409, 2020.
11

---

[17] Jonas Köhler, Leon Klein, and Frank Noé. Equivariant flows: sampling configurations for multi-body systems with symmetric energies. arXiv preprint arXiv:1910.00753, 2019.
[18] Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E(n) equivariant graph neural networks. In ICML, volume 139 of Proceedings of Machine Learning Research, pages 9323- 9332, 2021.
[19] Wenbing Huang, Jiaqi Han, Yu Rong, Tingyang Xu, Fuchun Sun, and Junzhou Huang. Equiv- ariant graph mechanics networks with constraints. In ICLR, 2022.
[20] Jiaqi Han, Yu Rong. Tingyang Xu, Fuchun Sun, and Wenbing Huang. Equivariant graph hicrarchy-basee neural networks. In ICLR, 2022. [tttlttt///ppeeeeee wtwpe/nrrviieeeeeefooreeee net/forum?
id=ywxtmGlnU_6
[21] Bowen Dai and Chris Bailey- Protein interaction interface region prediction by geo- metric deep learning. Bioinformatiex, 2021.
[22] Alex Fout, Jonathon Byrd, Basir Shariat, and Asa Ben-Hur. Protein interface prediction using graph convolutional networks. In Proceedings of the s1 International Conference on Neural Information Processing Systems, pages 6533-6542, 2017.
[23] Min Zeng, Fuhao Zhang, Fang-Xiang Wu, Yaohang Li, Jianxin Wang, and Min Li. Protein- protein interaction site prediction through combining local and global features with deep neural networks.
Bioinformatics, 36(4): 14-1120, 2020.
[24] Helen M Berman, John Westbrook, Zukang Feng, Gary Gilliland, Talapady N Bhat, Helge Weissig, Ilya Shindyalov, and Philip Bourne. The protein data bank, Nucleic acids research, (8(1):235-243 2000.
[25] Kendall Atkinson. An introduction to numerical analysis. John wiley & sons, 1991.
[26] Serge Bernstein. Sur I' l'ordre de la meilleure approximation des fonctions continues par des polynômes degré donné, volume 4. Hayez, imprimeur des académies royales, 1912.
[27] CMU. Carnegie-mellon motion capture database. 2003. URL httthtttt//maappa ca.ccap..cc.... cmu
cdu edu
[28] Stefan Chmiela, Alexandre Tkatchenko, Huziel Sauceda, Igor Poltavsky, Kristof T Schütt, and Klaus-Robert Robert Müller. Machine learning of accurate energy-conserving molecular force fields. Science advances, 3(5):e6630015, 2017.
[29] Fabian Fuchs, Daniel E. Worrall, Volker Fischer, and Max Welling. transtrrrassss:mrs::: 3d troo-translation equivariant attention networks. In NeurIPS, 2020.
[30] Peter Battaglia, Razvan Pascanu, Matthew Lai, Danilo Jimenez Rezende, et al. Interaction networks for learning about objects, relations and physics. In NeurlPS, volume 29, 2016.
[31] Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haher, Li F Fei-Fei, Josh Tenenbaum, and Daniel L Yamins. Flexible neural representation for physics prediction. In NeurlPS, volume 31, 2018.
[32] Tobias Pfaff, Meire Fortunato, Alvaro Sanchez-Gionzalez, and Peter Battaglia. Learning mesh- based simulation with graph networks. In ICLR, 2021. URL httpps://openreriex net/ Forum?id=rol@yYLo. XP
[33] Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathematical statistics. 22(1):79-866 1951.
[34] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow. NeurIPS, 29, 2016.
[35] Casper Kaae Sønderby, Tapani Raiko, Lars Maaløe, Søren Kaae Sønderby, and Ole Winther. Ladder variational autoencoders. NeurlPS, 29, 2016.
[36] Thomas Kipf, Ethan Fetaya, Kuan-Chieh Wang, Max Welling, and Richard Zemel. Neural relational inference for interacting systems. In ICML, pages (688-2697, 2018.
12

---

[37] Alvaro Sanchez- Gonzalez, Victor Bapst, Kyle Cranmer, and Peter Battaglia. Hamiltonian graph networks with ode integrators, arXiv preprint 19X91:1909.120 12790, 2019.
[38]
Johannes Brandstetter, Rob Hesselink, Elise van der Pol, Erik J. Bekkers, and Max Welling. Geometric and physical quantities improve E(3) equivariant message passing. In ICLR, 2022.
URL httpps::/ppperreeiee wtw.://www..ee/oomomiiim - 1
[39] Nathaniel Thomas, Tess Smidt, Steven Kearnes, Lusann Yang, Li Li, Kai Kohlhoff, and Patrick Riley. Tensor field networks: Rotation-and translationeequivarian neural networks for 3d point clouds. arXiv preprint arXiv:1802.082219 2018.
[40] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In NewrIPS, pages 6572-6583, 2018.
[41] Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. arXiv preprint arXiv: arXiv: 1.07522, 2019.
[42] Carl Runge. Über empirische funktionen und die interpolation zwischen äquidistanten ordinaten. Zeitschrift für Mathematik und Physik, 46(224-243):20, 1901.
[43] Ronald J Williams and David Zipser. learning algorithm for continually running fully recurrent neural networks. Neural computation, 1, 1989.
[44] Lingbing Guo, Zequn Sun, and Wei Hu. Learning to exploit long-term relational dependencies in knowledge graphs. In ICML, 2019.
[45] Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural Computation, 1997.
[46] Rafal Józefowicz, Oriol Vinyals, Mike Schuster, Noam Shazeer, and Yonghui Wu. Exploring the limits of language modeling. In ICLR, 2016.
[47] Oriol Vinyals, Lukasz Kaiser, Terry Koo, Slav Petrov, Ilya Sutskever, and Geoffrey E. Hinton. Grammar as foreign language. In NeurIPS, 2015.
[48] Isaac Newton. The Principia: mathematical principles of natural philosophy. Univ of California Press, 1999.
[49] Chenxin Xu, Robby T. Tan, Yuhong Tan, Siheng Chen, Yu Guang Wang, Xinchao Wang, and Yanfeng Wang. Eqmotion: Equivariant multi-agent motion prediction with invariant interaction reasoning. In CVPR, pages 1410-1420. IEBE, 2023.
[50] Thomas N. Kipf, Ethan Fetaya, Kuan-Chiel Wang. Max Welling, and Richard S. Zemel. Neural relational inference for interacting systems. In ICML, volume 80 of Proceedings of Machine Learning Research, pages (693-2702. PMLR, 2018
[51] Chenxin Xu, Maosen Li, Zhenyang Ni, Ya Zhang, and Siheng Chen. Groupnet: Multiscale hypergraph neural networks for trajectory prediction with relational reasoning. In CVPR, pages 6488 6497. IEEE, 2022.
[52] Miltiadis Kofinas, Naveen Shankar Nagaraja, and Efstratios Gavves. Roto-translated local coordinate frames for interacting dynamical systems. In NeurlPS, pages 6417-6429, 2021.
[53] Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
[54] Lei Jimmy Ba, Jamie Ryan Kiros, and Geoffrey E. Hinton. Layer normalization. CORR, abs/1607.06650. 2016.
[55] Vinod Nair and Geoffrey E Hinton. Rectified linear units improve restricted boltzmann machines. In ICML, 2010.
13

---

A The Derivations of NC (0) for Other EGNN Models
In this section, first illustrate how the other two baseline models (i.e. RF 17 and GMN 19 be derived to NC (0), and then discuss multi-layer EGNN.
A.1 RF
The overall equivariant convolution process of RF is similar that of EGNN. The main difference is that RF does not use the node feature. Instead, it leverages the L2 norm of the velocity as an additional feature to update the predicted velocity:
v. = Norn(v))@
1 + N
Σ (xf x)mm?; j
(19)
x x! + TTT,
(20)
where Norm(v denotes the L2 norm of the input velocity vp. Therefore, RF can be also viewed as a NC (0) model.
A.2 GMN
The main difference between GMN and EGNN is that GMN detects the sub-structures in the system and process the particles in each special sub-structure independently. Specifically, GMN re-formulates the original EGNN (i.e., Equations
follows:
m (h; hj, ||x x911² exj):
(21)
v e h.vvo + Σ
dn (x mly),
(22)
iES
iES
v = FK(V
(23)
X = x + T.TT
(24)
where vi the predicted velocity of the sub-structure. GMN then uses calculate the velocity of each particle by a function FK which can be either learnable or based on the angles and relative positions in the sub-structure [
A.3 Multi-layer EGNN
In current EGNN methods, the input coordinate x and velocity v are regarded constant feature and used in different layers. For example, in the multi-layer version EGNN, my, will be formulated as follows:
0,1 m = dhe (h .b- .|xx x9I1.en))
(25)
0./ where m.; is my, in layer 1 and h denotes the hidden feature in layer / - 1. Evidently, only the hidden features are updated within different layers. The overall formulation multi-layer layer akin to the single-layer EGNN.
14

---

B
1 Proofs of Things
B.1 Proof of Proposition 3.1
3.1
Proof. The objective of the existing methods for single system can be defined as:
argmin
Σ
{x
x
(26)
P1
Σ
[x]
x
T)
(27)
P1
T x.
T(
(28)
T
Pi
Ejzi Σ (vf -doo(h,) +
-x)mm (29)
=T
N
T Σ (v. - (w +b)),
(30)
where wD E R! and bd E R denote the learnable variables irrelevant to vi and t, concluding the proof.
B.2 Proof of Proposition 3.3
Proof. As the higher order cases (k: 2 1) have already been proved, we only need to show that
exc(O) 2 Exc(1). The first order of Newton-Cotes formula NC (1) is also known Trapezoidal rule, ie.:
T (vo + v²)). [- v(t)at
(31)
As aforementioned, the actual integration x" x for different training examples is different, and we assume that it fluctuates around the base estimation (v + +v") and follows a normal distribution NNC(1), where the variance of ()C))) is positively correlated with the difficulty of optimizing the overall objective. The variance of NNC(I) is:
E, Ch
Ei_ (kt()))dt)2
NC(L)
(32)
NT
where the integration term general form of polynomial interpolation error. According to Equa- tion 14 it can be derived to:
[" (1 (1-
(- (E)
)al,
(33)
2!
where v" denote the second derivative of v. Let 8 Tds. the above equation can re- re-written as:
t-to 7
then t = 10 + sh and di = d(v" sh)
8(8 1)T2v" (E)
T
ds
3 (E) = 0
(34)
12
0
Therefore, the final o
NC(1) is:
No = Ep 17 T (())
(35)
NT2
O(T*) < O(T2) = = :NC)))
(36)
concluding the proof.
15

---

Differsse (N))
Diffbranne
FPach
Epoch
NC(2)
(b) NCC3)
NC(5)
Figure 6: The prediction errors of intermediate velocities MD17 dataset, w.r.t training epoch. The blue and green lines denote the prediction errors of NC and NC respectively.
Difference (5)
Epoch
Epoch
NC(2)
(6) NC(5)
Figure 7: The prediction errors of intermediate velocities N-body dataset, w.r.. training epoch.
The blue and green lines denote the prediction errors NC and NC+, respectively.
NC12) intermediate
(b) NC13) intermediato
NCX2) final
NC(5) intermediate
Figure 8: Visualization of the intermediate velocities w.r.t. k of MD17 (Aspirin). The red, blue, and green lines denote the target, prediction of NC, and prediction of NC+, respectively.
B.3 Proof of Proposition 3.4
Proof. The GNN models possessing equivariance property are equivariant to the translation, rotation, and permutation input. NC directly feeds the input into these backbone models and naturally possesses this property.
Formally, let Tg equivariant then we will have:
be group of transformation operations. If the backbone model M is
(x)((x)) = (M))))).
(37)
where S, is an equivalent transformation to T, on the output space. NC can be regarded as weighted combination of the outputs of M:
Σ (x) Σ (wss((((x)),
(38)
where the Newton-Cotes weights w; constant and irrelevant to the input, the output, and the model M itself. Therefore, the above equation will always holds.
C Additional Experimental Results
The average intermediate velocity prediction errors on MD17 and Motion datasets are shown in Figure 6 and Figure 7 respectively. NC still learned the intermediate velocities did not feed
16

---

Table 6: Hyper-parameter settings in the main experiments.
Datascts
k
velocity regularization
velocity regularization decay
parameter regularization
parameter regularization decay
criterion
normalization
feature
normaiization
intermediste
velocity
N-body
0.001
0.999
1.0
0.99
MSE MSE MSE
False
True
MD17
0.1
0.999
1.0
0.95
True
True
Motion
0.01
0.999
1.0
0.95
True
True
clip
iraining examples
learning rale
# opoch
batch-size
activation
layers
optimizer
gradient norm
N body
1,500
200
500
ReLUU
0,00015 0,0005
Adam
1.0
ReL.U
Motion
1,500
100
200
ReLU
Adam
1.0
the intermediate into it. Particularly, the error N-body dataset was small and stable, which may demonstrate the effectiveness estimating intermediate velocities even without supervised data.
We also provide some visualized examples on MD17 dataset in Figure 8 from which we can observe similar results compared with Figure in the main paper.
D Hyper-parameter Setting
We list the main hyper-parameter setting of NC with different EGNN models different datasets Table [ We used the Adam optimizer [3] and adopted layer normalization [54] and ReLU activation (55) for all settings. For fair comparison, the parameter settings for the backbone EGNN models were identical to these in the existing implementation [1
17


Reviews:
[{"review_id": "H0h9zErYAS", "reviewer": "Reviewer_P4zh", "summary": "The paper introduces a new method for predicting system dynamics. Unlike past works that predict the average velocity between start and end time, the proposed approach makes several velocity estimation at different time points and then uses Newton-Cotes formula for computing the integration. This method reduces the prediction error and improve accuracy in predicting future state of the system. The paper improves this method further by incorporating loss for each predicted velocity estimate, and shows better performance than current SOTA systems.  ", "questions": "I would like authors to address the points raised in weakness section. What is the intuition behind accuracy not increasing past K=2. It seems that higher K should result in even better performance, but that's not the case. ", "limitations": "The authors have not addressed limitations in their work. ", "rating": 6, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "1. The paper is clearly written with good figures. Supplementary material provides additional GIFs and proofs to help with understanding.\n2. The paper proposes an interesting method to improve the current dynamics prediction systems by predicting velocity estimates at different time points and then using Newton-Cotes method to integrate them. Paper proofs that this method leads to reduction in variance, and therefore the error. NC+ approach further reduces the error by incorporating additional loss.\n3. The paper shows that even their NC method outperforms the previous SOTA approaches significantly.\n4. Despite predicting velocities at different time points, this method doesn't degrade the training time a lotl.", "weaknesses": "I feel the novelty of the method is limited. The method replaces the average velocity estimate with multiple velocity estimates that are then integrated over, something that seems pretty obvious. Also, something that is not clear is why the performance improvement saturates at only K=2, as shown in Fig 2? Another weakness is that the inference time scales linearly as the value of K is increased. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your detailed and insightful comments. We have addressed your concerns below and hope our responses provide clarity:\n\n**Weaknesses:**\n\n**1. The novelty of the method may be limited. The method seems intuitive.**\n\n  The novelty of our paper lies not just in the design of neural models but also in the theoretical exploration of EGNNs. We reformulate the existing EGNN methods in the form of numerical integration, which offers a more comprehensive view of reasoning system dynamics.\n \n  NC (k) is straightforward, yet powerful and fully supported by the theoretical proof. It significantly improves the performance of EGNNs on all datasets. Thanks to its simplicity, the computational complexity is also acceptable. The training time increases slightly as it estimates more points (as shown in Fig. 2).\n\n**2. Why does the performance improvement saturate at K=2, as shown in Fig 2?**\n\n  Sorry for the confusion and we will improve the presentation. As it may not be very clear from Fig. 2, we provide the detailed results in Table 1 in the rebuttal PDF. The prediction error is still reduced (slightly) as we set a larger k. This behavior could be attributed to the following reasons: (1) The neural model learns to compensate for the prediction error (as illustrated in Fig. 1 and Lines 36-44). Setting a larger k may reduce the variance of the prediction error and make the optimization easier, but this strategy does not always work, especially when the model is robust. (2) The accuracy gain exponentially decreases as we estimate more points in numerical integration. Take Fig. 1 as an example, a shift from 1e-1 to 1e-2 yields a more significant gain compared to a shift from 1e-3 to 1e-4. (3) Runge’s oscillation phenomenon (Lines 152-154) tells us that more interpolation points may not always lead to more accurate results.\n\n**3. The inference time scales linearly with the increase in k. The authors haven't addressed this limitation.**\n \n  There seems to be a misunderstanding. As shown in Fig. 2 (rows 2 and 3), the training time only increases slightly. We provide the detailed training time on N-body (1, 2, 0) below:\n\n| k  | GMN   | EGNN   | RF   |\n|---|:---:|:---:|:---:|\n| 0  | 1817.0   | 1355.4   |  1328.2  |\n| 1  |  1800.9 | 1357.4  |  1333.2 |\n| 2  |  2057.9 | 1395.6  |  1328.9 |\n| 5  |  3749.0 | 1841.2  | 1946.9  |\n\n  Please note that k is in {0, 1, 2, 5}. \n\n**Questions:**\n\n**1. I would like authors to address the points raised in the weakness section.**\n\n  Thank you again for your feedback. We hope our responses to **Weaknesses** have addressed your concern.\n\n\n\n**Limitations:**\n\n**1. The authors have not addressed limitations of their work.**\n\n  We acknowledge this concern and have indeed discussed this issue in Sec. 3.5. As elaborated in our response to Weakness 3, the increased computational cost with higher k values is affordable. We believe this trade-off is worth the advantages offered by our method.\n"}, {"review_id": "v1JEn4Sq1m", "reviewer": "Reviewer_MNf4", "summary": "This paper integrates Newton-Cotes formulas from numerical analysis in (equivariant) graph networks. The authors introduce an approach that makes several velocity estimations to predict the integration towards the final state of a system.\nFurther, they show that standard equivariant graph networks can be interpreted as a trivial case of their formalism. The proposed method (NC) has provably lower errors/variance, and becomes better with more velocity estimations. NC can be used with supervision from intermediate velocities, but it can also work well without supervision.\nThe authors demonstrate the effectiveness of their method on several benchmarks, and using multiple graph network architectures as a backbone.", "questions": "1. Following the weaknesses listed above, in figures 3 and 6, why do the errors fluctuate or keep increasing with epochs? What is its impact on performance?\n\n2. The Newton-Cotes formulas exhibit a reflection symmetry. Is it possible that the model is learning to predict the reflected intermediate velocities? An ablation study on that would be useful.\n\n3. In subsection 4.6, the authors perform an ablation study on the method's ability to make consecutive predictions. There are related works that propose equivariant graph networks for multiple step prediction, including LoCS [2] and EqMotion [1]. It would be interesting to examine the applicability and added value of NC to such methods.\nAs an example, it would be interesting to see the integration of NC in LoCS, since this work treats velocities in a principled manner and makes predictions in local frames.\n\n\n#### Minor \n\n4. In lines 161-162, the authors mention that \"EGNN always uses the same input velocity and coordinate as constant features in different layers\". However, EGNN does use the updated coordinates in different layers, at least in the n-body setting.\n\n#### References\n[1] Xu, Chenxin et al. EqMotion: Equivariant Multi-agent Motion Prediction with Invariant Interaction Reasoning. CVPR 2023.\n\n[2] Kofinas, Miltiadis et al. Roto-translated Local Coordinate Frames for Interacting Dynamical Systems. NeurIPS 2021.\n", "limitations": "The authors have adequately addressed the limitations of their work.", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is very well-written and easy to follow. The proposed method is well-motivated from numerical analysis, and it is simple and clear.\nThe proposed method shows improvements on 3 different benchmarks, and using 3 established equivariant graph networks as a backbone.\nThe non-supervised version of the method (NC) can be directly integrated in many existing graph network architectures.", "weaknesses": "The analysis on the prediction errors of intermediate velocities (figures 3, 6, 7) could be more in-depth, especially for the unsupervised version of the method (NC). It is unclear why the errors often increase as training progresses, or why they tend to fluctuate erratically, and whether this has an impact on performance, or if the model can make up for it.\n\nThe qualitative visualizations, especially for MD17 and n-body (figure 8) are very hard to interpret.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are grateful for your encouraging comments and insightful observations. We hope the following responses address your concerns:\n\n**Weaknesses:**\n\n**1. The analysis on the prediction errors of intermediate velocities (figures 3, 6, 7) could be more in-depth, why the error of the unsupervised NC tends to fluctuate erratically, and whether this has an impact on performance, or if the model can make up for it.**\n\n  Thanks very much for your suggestion. To ensure visual clarity, we have magnified the prediction errors (%), which may make the unsupervised errors seem more erratic than they truly are; otherwise, the error of the supervised version would look like a zero line. \n\n  We provide a new figure (Figure 1 in the rebuttal PDF) depicting both the intermediate error and final prediction error. Despite the fluctuations, the intermediate error’s impact on the final prediction error is minimal, as the neural model compensates accordingly. \n\n  The unsupervised NC does not have an explicit loss to minimize the intermediate error (as stated in Lines 261-264). We want to use this experiment to illustrate that even without an explicit supervised loss, the model automatically adjusts to fit the intermediate points to obtain a better prediction for the final states.\n\n**2. The qualitative visualizations for MD17 and n-body (figure 8) could be more readable.**\n\n  Thanks for your suggestion. The input objects are molecules in MD17, but the base model GMN considers only a subset of chemical bonds during training and inference. This is why the objects in the figures do not look like molecules. We have updated the illustrations and included a new figure (Figure 2 in the rebuttal PDF). For the N-body figures, presenting a more readable version is challenging as the input object is a mixture of points, sticks, and hinges.\n\n**Questions:**\n\n**1. Following Weakness 1, why do the errors fluctuate or keep increasing with epoch and what is its impact on performance?**\n\nPlease see our response to **Weakness 1**. \n\n**2. Is it possible that the model is learning to predict the reflected intermediate velocities?**\n\n  Thanks for your suggestion. It is indeed possible to predict the intermediate velocity given the initial and final states. We perform a comparison with a GMN variant where we concatenate the outputs of the initial and final states to predict the intermediate velocity. The results on the motion dataset are shown in the table below:\n\n|  Method | Intermediate Velocity Prediction  | Final State Prediction  |\n|---|---:|---:|\n| GMN  | 5.4  |  43.9 |\n| NC  | 2.8  | 29.6  |\n\nWe can observe that our method employing NC formulas to inversely infer intermediate velocities significantly outperforms the GMN variant. A more interesting and challenging task may be predicting the continual trajectories between the initial and final states, which we intend to explore in the future.\n \n**3. There are EGNN methods for multiple step prediction, including LoCS [2] and EqMotion [1]. It would be interesting to examine the applicability and added value of NC to such methods.**\n\nWe have applied our method to EqMotion and obtained some preliminary results. To ensure a fair comparison, we use identical hyper-parameter settings and we do not employ any additional data. The results are shown in the following table (ADE and FDE denote average displacement error and final displacement error, respectively): \n\n| Method  |  N-body Prediction ADE | N-body Reasoning Accuracy  | MD17 (Aspirin) ADE/FDE  | MD17 (Benzene) ADE/FDE  | MD17 (Ethanol) ADE/FDE  | MD17 (Malonaldehyde) ADE/FDE | \n|---|---|---|---|---|---|---|\n| EqMotion | 16.8  | 97.8  | 5.95/8.38  | 1.18/1.73 | 5.05/7.02 | 5.85/9.02 | \n| NC (EqMotion)  | **16.2**  | **98.9**  | **4.74/6.76**  | **1.15/1.63** | **5.02/6.87** | **5.19/8.07** | \n\n\nClearly, our NC (EqMotion) outperforms the base model EqMotion on all datasets and across all metrics. We will also release the source code of NC (EqMotion) if the paper gets accepted. Thank you again for your insightful comment.\n\n\n**4. Minor - In lines 161-162, the introduction to EGNN should be revised. It uses the updated coordinates in different layers at least in the n-body setting.**\n\nMany thanks. We have corrected it."}, {"review_id": "uxh3rYab7L", "reviewer": "Reviewer_Gu3W", "summary": "This paper builds on recent advances in deep learning for physical simulation to produce reliable results of long-term consecutive predictions. The authors provide a proof that equivariant graph neural networks learn a linear map from initial velocity to a target velocity, and propose to inject inductive bias based on Newton-Cotes (NC) formulas as alternative of standard equivariant graph convolution. The paper introduces techniques for training, such as a new objective and a new regularization loss for higher order features, and also gives theoretical guarantee that variance can be reduced by increasing the degree of NC formula. Experiments in various multiple body interaction systems demonstrate that the proposed model outperforms  state-of-the-art baselines .", "questions": "Please have a look at the weakness mentioned in the strength and weakness section and address these. Overall the idea is interesting, the authors need to substantiate their claims in light of possibly a few more experiments.\n\n**Minor question**\n\nWhat does “steps” of Table 4 in Appendix C mean? ", "limitations": "Yes, the limitations are discussed in the appendix. ", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "The paper is appropriately placed in the current growing literature on scientific machine learning to create a novel fast model for tackling complex problems in the realm of particle systems. Although the proposed architecture is simple in appearance, the motivation of introducing NC formulas is convincing giving both concise visualization and theoretical claims on the limitation of the existing models. Results of the experiments show that introducing NC formulas is effective to perform long-term accurate prediction. \n\nIn my opinion this is a good paper and above the acceptance threshold.", "weaknesses": "Although the proposed methods NC(k) and NC(k)+ outperforms state-of-the-art baselines, it is still somewhat unclear which components of NC(k)+ help to improve the performance of NC(k)+. What if setting all the coefficients of Newton-Cotes formula to be 1 (or some constant) while remaining the other parameters same? Without some additional ablation study, it is a bit difficult to get a clear picture of how new components of the proposed method would help improving its prediction performance. \n\nIn the forward process, how is the coefficients of Newton-Cotes formula obtained? Is the coefficients trainable or derived from the equation (15)? It is also unclear how to derive coordinate at intermediate time step $t^{k}$ once velocity is updated by the formula (13), while the formula (12) is used only when computing the final state at T. \n\nFor training time shown in Figure 2, why does the training cost of NC with GMN grow much faster than NC with the other backbones? It will be helpful to get more insight of the trends. Also, the order of NC(k) seems to be set K<=8, due to the catastrophic Runge’s phenomenon. Are there any limitation on the length of time interval [0, T] to which new NC points are spaced? How does the length of the interval affects the performance?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Many thanks for your helpful feedback and suggestions. We will carefully incorporate them into our paper.\n\n**Weaknesses:**\n\n**1. Performing a comparison with a variant NC model with constant coefficients may be helpful to better understand the proposed method.**\n\n  Great idea. We implement a new variant where the coefficients are set to 1 and other hyper-parameters remain unchanged. We report the results as follows:\n\n| Method  |  Intermediate Velocity Error | Final Prediction Error  |\n|---|---:|---:|\n| GMN  | -  | 10.14  |\n| NC (coefficients=1)  | 100.50  | 9.64  |\n| NC   | 6.37  | 9.50  |\n|  NC+ (coefficients=1)  |  10.53 | 9.91  |\n|  NC+   | 0.73  | 9.40  |\n\n  For NC, the model does not implicitly fit the NC formulas when we set all coefficients to 1.  As we can observe, it obtains a large intermediate velocity error. In this case, the model may find a new intermediate value to fit the new coefficients.  It still performs better than the base model GMN, but worse than the original NC. According to NC formulas, if the initial and final points are set, the optimal intermediate points and their coefficients are also determined. Setting all coefficients to 1 may be a worse hint to the model. \n\n  For NC+, we set all coefficients to 1 and compel it to learn the original intermediate velocity. We can see that there is a contradiction between these two objectives. The model performed worse on either prediction. \nWe also report the results on the other two datasets in Table 2 in the rebuttal PDF. The proposed method with constant coefficients performs slightly better than GMN on the N-body (1, 2, 0) dataset, but fails to reason dynamics on the more complicated human motion dataset.\n\n  We also report the results on the other two datasets in Table 2 in the rebuttal PDF. The proposed method with constant coefficients performs slightly better than GMN on the N-body (1, 2, 0) dataset, but fails to reason dynamics on the more complicated human motion dataset.\n\n**2. In the forward process, how to obtain the coefficients of Newton-Cotes formula,  are the coefficients trainable or derived from Equation (15)? How to update the coordinate at intermediate time step once velocity is updated by Equation (13).**\n\n  The coefficients are derived from Equation (15) and are hence not trainable. We calculate the intermediate coordinate by adding the predicted velocity to the last coordinate during the inference phase. In Equations (12-13), we can directly feed NC+ with true preceding coordinates during the training phase (a common trick in training sequence models).\n\n**3. Why does the training cost of NC with GMN grow much faster than NC with the other backbones in Figure 2? Also, the order of NC(k) seems to be set K<=8, due to the catastrophic Runge’s phenomenon. Are there any limitations on the length of time interval [0, T] to which new NC points are spaced? How does the length of the interval affect the performance?**\n\n  Regarding the training cost, the GMN model generally requires more time compared to others (e.g., 0.8s per update v.s. 0.6 per update on N-body (1, 2, 0)). It considers the substructure (e.g., hinges and chains) within each object, which yields more cost when we consider more points. In fact, the training cost of GMN doesn’t grow so much faster than that of others. Please note that the y-ticks don’t always start from 0 in Figure 2 (for better readability).\n\n  Regarding the order of NC (k), the neural model learns to compensate for the prediction error (as illustrated in Fig. 1 and Lines 36-44). Setting a larger k may reduce the variance of the prediction error and make the optimization easier, but this strategy does not always work especially when the model is robust. Furthermore, the accuracy gain exponentially decreases as we estimate more points in numerical integration. Take Fig. 1 as an example, a shift from 1e-1 to 1e-2 yields a more significant gain compared to a shift from 1e-3 to 1e-4. \n\nWe will clarify these aspects further in the revision.\n\n**Questions:**\n\n**1. Overall the idea is interesting, the authors need to substantiate their claims in light of possibly a few more experiments.**\n\nWe appreciate your constructive feedback and hope that our responses have satisfactorily addressed your concerns.\n\n**2. Minor question -What does “steps” of Table 4 in Appendix C mean?**\n\nSorry for the confusion. “# steps” is a shorthand for “the number of estimation steps”, i.e., k in our main paper. We have corrected it for better clarity. \n\n"}, {"review_id": "59CAbqUE5D", "reviewer": "Reviewer_9Zk2", "summary": "This paper proposes Newton-Cotes Graph Neural Networks (NC), which are a time-continuous extension of Equivariant Graph Neural Networks (EGNNs). The authors also conducted a mathematical analysis on the predictive variance of NC under the Gaussian assumption. Through numerical experiments, the authors have demonstrated that the NC-formulated EGNNs outperform the baseline EGNNs.", "questions": "- Regarding the selection of the numerical integrator, it would be beneficial to clarify the advantages of using NC over more well-known and commonly used integrators such as Runge-Kutta or Dormand–Prince.\n- When it comes to training NC+, is it necessary to collect the K-1 intermediate points?", "limitations": "- The current manuscript lacks a comparison with the most competitive baseline methods, as well as a comprehensive review of relevant literature.", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "- The method is clear and can be generalized to other GNN classes beyond Equivariant GNNs.", "weaknesses": "**The proposed method is a special case of a larger GNN class, but there is no comparison provided with this GNN class in the submitted paper.**\n\nAlthough the authors have demonstrated the effectiveness of the NC integrator in conjunction with GNN through numerical experiments, the paper fails to compare the proposed Newton-Cotes Graph Neural Networks (NC) with the most competitive baseline class, Graph Neural Ordinary Differential Equations (GDEs) [1].\n\nIt would be highly beneficial to compare the proposed method with the more well-established GDE models. In my opinion, the proposed NC can be considered a special case of GDE where the integrator is the Newton-Cotes method and the GNN corresponds to the proposed EGNN. Furthermore, based on the advancements in Neural ODE research, the limitation mentioned in section 3.5 could already be addressed using adjoint sensitivity methods.\n\n[1] Graph Neural Ordinary Differential Equations", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your insightful comments. We have taken them into account and provided clarifications below:\n\n**Weaknesses:**\n\n**1. The proposed method can be considered a special case of GDEs. It would be highly beneficial to compare the proposed method with GDEs.**\n\nThanks very much for your insightful and detailed suggestion. We have added GDEs to our related works in the revision. GDEs propose a set of strategies to incorporate GNNs with neural ODEs:\n\nThe depth-varying GDE, which can be regarded as “continuous” ResNet + GNNs. The GNN module acts as the derivative and can be stacked arbitrary times. In contrast, our paper uses the Newton-Cotes coefficients as a kind of prompt to encourage the model to fit the intermediate values and minimize the training loss. We have observed that NC outperforms the baselines significantly (Tables 1, 2, and 3), even without the additional data. We have also observed that NC implicitly learns the intermediate velocities in Figure 3. In our response to Reviewer Gu3W’s comment (Weakness 1), we implement a new variant of NC with coefficients set to 1. This variant is similar to the depth-varying GDE, but it performs worse than NC. To thoroughly investigate this issue, we have modified the source code of GDEs for our task. We use an identical GMN as the backend EGNN model and perform a simple hyper-parameter search for the learning rate and ODE solver. We report the best results as follows:\n\n| Method  |  N-body (1,2,0) | MD17 (Aspirin)  | Motion  |\n|---|---:|---:|---:|\n| GMN  | 1.84  | 10.14  | 43.9  |\n| NC (coefficients=1)  | 1.78  | 9.64  | 48.4  |\n| depth-varying GDE  | 1.82  | 10.08  | 42.2  |\n| NC   | 1.55  | 9.50  | 30.0  |\n\nThe performance of depth-varying GDE is comparable to that of NC (coefficients=1), but both are lower than that of the original NC. All the methods generally perform better than the baseline GMN. This observation further supports the effectiveness of our method.\n  \nAnother related GDE model is the spatio-temporal GDE, but we are unable to adapt it to our task due to the unavailability of the source code. The spatio-temporal GDE, like many neural ODE methods, aims to build a continuous prediction and capture observations with time steps. However, it is not appropriate to consider our method as a special case of GDEs. To draw an analogy, it is acceptable to state that we use numerical integration algorithms to solve ODE problems, but it may sound peculiar to claim that we use ODE solvers to solve numerical integration. The numerical solvers (e.g., 4th order Runge-Kutta, rk4) in neural ODEs are originally designed for numerical integration. Regarding our method as a special case of GDEs is akin to treating numerical integration as an ODE problem and solving it using numerical algorithms. In contrast, it may be feasible to iteratively perform our method to solve problems in GDEs. Our work draws inspiration from observations on numerical integration algorithms, where the integrand is known only at certain points. The goal is to efficiently approximate the integral to the desired precision. Neural ODE methods repeatedly perform the numerical algorithms on small internals to obtain the numerical solution. \n\nFinally, we would like to emphasize that our contributions extend beyond proposing a new method. The theoretical analysis is equally important, as it provides a novel and comprehensive perspective on reasoning dynamics with neural models. \n \n**2. The limitation mentioned in section 3.5 could be addressed using adjoint sensitivity methods.**\n\nThe adjoint sensitivity methods are generally used to reduce memory usage (at a cost of time). They cannot reduce the training time of our method. \n\n**Questions:**\n\n**1. It would be beneficial to clarify the advantages of using NC over other algorithms.**\n\nThanks for your suggestion. We use NC because it aligns closely with our scenarios. The intervals between states are always constant in trajectories. This characteristic allows us to easily obtain the intermediate points that satisfy NC (but not other algorithms, e.g., Gaussian quadrature and Clenshaw–Curtis quadrature). We will include the above explanations in the revision.\n\n**2. When it comes to training NC+, is it necessary to collect the K-1 intermediate points?**\n\nWe acknowledge the need for the intermediate points to train NC+ in a supervised fashion. Fortunately, the additional data can be readily obtained since the sampling frequency used in constructing the datasets is much higher than what we require. \n\nOnce again, thank you for your valuable feedback. We hope our responses have adequately addressed your concerns, and we kindly request you to reconsider the evaluation in light of these points."}]

Author rebuttal:
Dear all reviewers:

We sincerely appreciate the time and effort you have dedicated to reviewing our paper. Your comments are invaluable in helping us enhance the quality of our paper.

We would like to thank Reviewers P4zh and MNF4 for their suggestions regarding the readability of some figures. We have addressed their concerns by updating the figure and including some examples in the rebuttal PDF.

We are also grateful to Reviewers Gu3W and 9Zk2 for their suggestions regarding new baselines. We have presented the new experimental results in the rebuttal PDF and our replies.

We especially thank Reviewer MNf4 for the suggestion about applying our method to the recent paper EqMotion [1]. We have conducted experiments and we are pleased to report that some preliminary results are highly promising:

| Method  |  N-body Prediction ADE | N-body Reasoning Accuracy  | MD17 (Aspirin) ADE/FDE  | MD17 (Benzene) ADE/FDE  | MD17 (Ethanol) ADE/FDE  | MD17 (Malonaldehyde) ADE/FDE | 
|---|---|---|---|---|---|---|
| EqMotion | 16.8  | 97.8  | 5.95/8.38  | 1.18/1.73 | 5.05/7.02 | 5.85/9.02 | 
| NC (EqMotion)  | **16.2**  | **98.9**  | **4.74/6.76**  | **1.15/1.63** | **5.02/6.87** | **5.19/8.07** | 

The empirical evidence from these new tasks further validates the effectiveness of our method. We are also committed to releasing the source code of NC (EqMotion) if our paper is accepted.


Best Regards,

Authors

[1] Xu, Chenxin et al. EqMotion: Equivariant Multi-agent Motion Prediction with Invariant Interaction Reasoning. CVPR 2023.


Meta-review: The paper proposes a new method called Newton-Cotes Graph Neural Networks (NC) to improve prediction of system dynamics over long time horizons. It uses multiple velocity estimations based on Newton-Cotes numerical integration formulas instead of a single velocity estimate used in prior work.
Several reviewers recommended acceptance, citing the clear writing, novel integration of numerical methods with graph neural networks, strong experimental results, and theoretical analysis. Their main criticisms were around novelty, limitations of the method, and some minor presentation issues.

Reviewer 9Zk2 initially recommended rejection, arguing the method is a special case of a broader class of models (Graph Neural ODEs) that was not compared to.
In their rebuttal, the authors added comparisons to depth-varying Graph Neural ODEs, showing their proposed NC method achieves better performance. They argued their method is complementary to, not a subset of, Neural ODE methods. Reviewer 9Zk2 was partially convinced but still questioned the necessity of using Newton-Cotes formulas specifically.

Through the author rebuttal, more justification for using Newton-Cotes was provided:
- Their data has fixed intervals between states, which makes it straightforward to select intermediate points satisfying Newton-Cotes conditions. This would be more difficult for methods like Runge-Kutta that use adaptive step sizes.
- Newton-Cotes allows them to easily incorporate intermediate supervision through NC+, which helps improve performance.
- The simplicity of Newton-Cotes makes it easy to integrate into existing graph network architectures.

Overall, I believe the paper does make solid contributions experimentally and theoretically. So I recommend acceptance, while encouraging the authors to expand on the rationale for Newton-Cotes in the camera-ready version.
Decision: Accept (spotlight)