Graph Spline Networks for Efficient Continuous Simulation of Dynamical Systems
Authors: Chuanbo Hua
Federico Berto
Michael Poli
Stefano Massaroli
Jinkyoo Park
Venue: Submitted to ICLR 2023 (2023)
Keywords: Graph
Spline Collocation Method
Graph Neural Networks
Simulation
Partial Differential Equations
PDEs
Physics
Scientific Computing

Abstract: While complex simulations of physical systems have been widely studied in engineering and scientific computing, lowering their often prohibitive computational requirements has only recently been tackled by deep learning approaches. In this paper, we present GraphSplineNets, a novel deep learning approach to speed up simulation of physical systems with spatio-temporal continuous outputs by exploiting the synergy between graph neural networks (GNN) and orthogonal spline collocation (OSC). Two differentiable time-oriented OSC and spatial-oriented OSC are applied to bridge the gap between discrete GNN outputs and generate continuous solutions at any location in space and time without explicit prior knowledge of underlying differential equations. Moreover, we introduce an adaptive collocation strategy in space to enable the model to sample from the most important regions. Our model improves on widely used graph neural networks for physics simulation on both efficiency and solution accuracy. We demonstrate SplineGraphNets in predicting complex dynamical systems such as the heat equation, damped wave propagation and the Navier-Stokes equations for incompressible flow, where they improve accuracy of more than 25% while providing at least 60% speedup. 

Under review as a conference paper at ICLR 2023
GRAPH SPLINE NETWORKS FOR EFFICIENT CONTINU- OUS SIMULATION OF DYNAMICAL SYSTEMS
A Anonymous authors
Paper under double-blind review
ABSTRACT
While complex simulations of physical systems have been widely studied in engi- neering and scientific computing, lowering their often prohibitive computational requirements has only recently been tackled by deep learning approaches. In this paper, we present GRSS novel deep learning approach to speed up simulation of physical systems with spatio-temporal continuous outputs by exploiting the synergy between graph neural networks (GNN) and orthogo- nal spline collocation (OSC). Two differentiable OSC (time-oriented OSC and spatial- OSC) are applied to bridge the gap between discrete GNN out- puts and generate continuous solutions any location in space and time without explicit prior knowledge of underlying differential equations. Moreover, we in- troduce adaptive collocation strategy in space enable the model sample from the most important regions. Our model improves on widely used graph neu- ral networks for physics simulation on both efficiency and solution accuracy. We demonstrate GRAPHSPLINNNETT in predicting complex dynamical systems such as the heat equation, damped wave propagation and the Navier-Stokes equations for incompressible flows, where they improve accuracy of more than 25% while providing at least 60% speedup.
INTRODUCTION
For a growing variety of fields, simulations of partial differential equations (PDEs) representing physical processes are an essential tool. PDE-based simulators have been widely employed in range of practical issues, spanning from astrophysics (Mücke et al., 2000) to biology (Quarteroni & Veneziani, 2003), engineering (Wu & Porté-Agel, 201 1), finance, (Marriott et al., 2015) or weather forecasting (Bauer et al., 2015). Traditional solvers for phsysics-based simulation oftentimes need significant amount of computational resources (Houska et al., 2012), such as solvers based on first principles and the modified Gauss-Newtor methods. To broaden the scope of applications of dynamics simulation, the scientific machine learning community has put considerable effort into developing computationally simple yet accurate simulation approaches.
Deep learning has been shown be powerful alternative efficiently compute solutions (Raissi et al., 2019) or model dynamical systems directly from data (Mrowca et al., 2018). Among deep learning methods, graph neural networks (GNNs) come with desirable properties such spatial equivariance and translational invariance which allow learning representations of dynamical inter- actions in a generalizable manner (Pfaff al., 2021; Bronstein et al., 2021) and unstructured grids. Despite the benefits of these paradigms, graph-based models have the fundamental drawback of being discrete nature, which makes it challenging implement continuous simulations time and space. While graph models that operate in continuous space continuou time have been intro- duced in the past (Poli et al., 2019), such approaches mainly deal with only one aspect of continuity at once, either in space time, and are hindered by accuracy issues while interpolating in space (Alet et al., 2019) or require a considerable number of iterative evaluations vector field time limiting their performance (Xhonneux et al., 2020).
To bridge the gap between the inherently discrete graphs and the intrinsic continuous nature of the real world, in this work we propose GRAPHSPLINENETS, novel method that exploits the syn- ergy between graph neural networks and the orthogonal spline collocation (OSC) method (Bialecki & Fairweather, 2001; Fairweather & Meade, 2020). By leveraging the OSC, approach can

---

Under review as a conference paper at ICLR 2023
produce predictions at any location in space and time without explicit prior knowledge of the under- lying differential equation. GRAPHSPLINENETS achieve significant speedups compared to GNN baselines by making use of efficient sparse linear solvers (de Boor & Weiss, 1980) for the OSC problem and training the model end-to-end with larger temporal resolutions. Moreover, thanks to the super-convergent approximations at nodes of the OSC partition (Qiao et al., 2021) and an adap- tive sampling strategy of collocation points, GRAPHSPLINENETS improve the solution accuracy of predictions in continuous space and time.
We summarize our contributions as follows:
We introduce GRAPHSPLLINNNTTT continuous time and space leveraging the OSC method.
learning framework for complex dynamical system in
We introduce an adaptive collocation sampling strategy to improve accuracy and a differentiable algorithm for fast inference of the OSC that allows for end-to-end training.
We demonstrate that GRAPHSPLINENETS outperform or competitive against baselines in predicting continuous complex dynamics in terms of both accuracy and speed.
2 RELATED WORKS
Graph Neural Networks for Dynamics Predictions Deep neural networks have recently been successfully employed in a variety of different tasks, ranging from simulated (Long et al., 2018; Li al., 2020) and real datasets (Pathak et al., 2022; Li et al., 2022a; Poli et al., 2022) demonstrating their capabilities in predicting complex dynamics often orders of magnitude faster than traditional numerical solvers. We aim finding efficient and accurate surrogate models: unlike data-driven approaches for solving PDEs such as PINNs (Raissi et al., 2019), that aim at finding solutions to set of equations, our methods does not need to know the exact equations of dynamical system and can directly learn mappings from data. One major line of work for dynamics prediction involves the use of graph neural networks (GNNss): these models provide several benefits compared to other deep learning methods based on regular grids such as convolutional networks. In particular, they make it possible to learn on irregular grids and varying connectivity and inherit physical properties derived from geometric deep learning, such as permutation and spatial equivariance (Bronstein et al., 10. Alet et al. (2019) represent adaptively sampled points in a graph architecture to simulate continuous underlying physical processes without any a priori graph structure. Sanchez-Gonzalez et al. (2020) introduce particle-based graph simulators with dynamically changing connectivity simulating inter- actions through message passing; Pfaff et al. (2021) extend particle-based simulations to base- ones. GNN-based approaches have also been shown to represent some parts of classical numerical solvers, such as finite differences and volumes (Brandstetter al., 2022; Lienen & Günnemann, 2022). Graph neural networks have also recently been applied to large-scale weather predictions (Keisler, 2022). Another deep learning direction on irregular grids that avoids graphs altogether is to convert the input domain into a regular grid via learnable deformations to make usage of neural operators possible (Li et al., 2022b). Compared to other deep learning methods for physics predic- tions on unstructured grids, we do not need to learn transforms and predict directly in the target domains; moreover, we place emphasis on bridging the gap between discrete graph nodes in space and time by allowing for fast and accurate continuous predictions.
Collocation Methods and Graphs Collocation and interpolation methods' are used to estimate unknown data values from known ones (Bourke, 1999). GNNs for predicting dynamics inherently lack important aspect characterizing physical systems: continuity. The concept of continuity can be separated into two categories: continuity in time and continuity in space. The former has been investigated by using continuous ODE models (Poli al., 2019; Xhonneux et al., 2020; Fang et al., 2021) that can, in theory, represent a system evolving in time continuously. However, such methods employ numerical solvers that introduce a considerable number of function evaluations and do not consider the fact that deep learning models can be capable of overcoming time discretization
CCollocation and interpolation are terms that are oftentimes used interchangeably. While interpolation is defined as obtaining unknown values from known ones, collocation is usually defined as a finite solution space satisfying equations dictated by known (collocation) points. Thus, collocation can be considered flexible subset of interpolation methods that satisfics certain conditions, such as C¹ class continuity.
2

---

Under review as a conference paper at ICLR 2023
errors by learning residual values (Poli et al., 2020; Berto et al., 2022). In approach, aim at learning state updates directly rather than vector fields to allow for fast predictions. In the case of space continuity, the graph structure itself is inherently discrete. Alet et al. (2019) employ linear interpolation in graph learning; however, this interpolation method is known to be inaccurate and does not respect C² continuity. Unlike earlier methods, we employ the orthogonal spline collocation (OSC) (Bialecki & Fairweather, 2001) method to quickly find C² continuous solutions to differential equations given small number of partition points, producing space-time continuous simulation. Unlike methods requiring full matrix inversion, OSC has a complexity of only O(n² log n)) thanks its sparse structure, allowing for fast inference. Another strain of work involves directly learning collocation weights (Guo et al., 2019; Brink et al., 2021); ); in our case, however, we use graphs to predict spatio-temporal locations of nodes and employ the OSC approach an end-to-end manner. A further benefit of the OSC choice are its theoretical guarantees on convergence (Bialecki, 1998).
SPLINEGRAPHNETTS
Message Passing Neural Networks (MPNN)
XX, fo
Encoder
Processor
Decoder
MPNN
f A
[Noo Vil
tAr
ts
(xen al2
y
12
y
(xc 1181
Space-oriented OSC
Time-oriented OSC
Linese Equation with ABD matiix
y
Linerr Equation with ADD matrix
Sinvalation Equation
Simulition Fquation ( + At)
Figure The overall scheme GRAPHSPLINENETS.
3.1 PROBLEM SETUP
We first introduce necessary notation that we will use throughout the paper. Let superscripts de- note time and subseripts denote space indexes well others with abuse of notation. We denote the state of a PDEs process at physical space location X = i I.... N) at time as Y {g/ 1.... NN} where N represents the number of sample points and is the physical domain. More specifically,
have x; Ω RDD. yf I can be described by solution of the PDES, i.c. yf = u(x,,t), t). The final objective of a physics process simulator is to estimate future states y given the history states (Ý', 1.... ,T). Long-term predictions can be obtained via autoregressive rollouts of the model. our method, we also aim infer a spatio-temporal continu- ous prediction TT+A! aT+A((x)) (x), S2 At E R,, where the aT+At the simulated functions of the domain at time T + At.
3.2 METHOD OVERVIEW
Fig. depicts the entire architecture our model that can divided into three main components: message passing neural networks (MPNN), space-oriented collocation and time-oriented colloca- tion. The MPNN takes history state observations as input and infers a sequence discrete future predictions via autoregressive rollouts. On these discrete predictions, then use the time oriented and space-oriented collocation methods to obtain simulation functions that can provide both time and space continuous simulations.
time-oriented

---

Under review as a conference paper at ICLR 2023
3.3 MESSAGE PASSING NEURAL NETWORKS
We employ graph G' (V. E' G to represent the nodes in the physics domain, where N V = {v:] == and E' = (ef LN attribute of the directed edge between node and j. respectively. Each node has the attribute encoded structure to predict
j- v and ely denote the attribute of sample node i and the
from its state information. The MPNN employs an encoddr-prrocesor-dcccooo the states of sample points the next timestep:
y D
Ca.....(P, YOJ)))
decoder
processor
encoder
where E-) is the encoder, P.(-) is the i-th message passing layer, and (()) is the decoder.
3.4 ORTHOGONAL SPLINE COLLOCATION METHOD
The orthogonal spline collocation (OSC) method consists of four steps in total: (1) partitioning and selection of collocation points (2) generating the equation of simulator polynomial parameters (3) solving equations and (4) simulating in the physical domain. In the following part of this section, introduce the OSC method for the 1-D and 2-D cases, as both are relevant to the proposed model. GRAPHSPLINENETTS continuously predicts in time domain using the 1-D OSC and in space domain using its 2-D counterpart.
Time-oriented OSC For a process of a specific sample point x, in the physics domain with its state changing over time, we can consider it as an ordinary differential equation (ODE) process [(),(()) = 0, E [0, 7), where (.).) is the solution with boundary conditions (,( = y ((T)
y/ The target of time-oriented OSC is find a series of polynomials under order r and satisfy C continuity to simulate the solution.
To find these polynomials, fp.N, = T. Note that these partitions can be not isometrie. Then, we initialize one polynomial of order r in each partition. These polynomials have (r - 1) degrees of freedom in total, which is the number of variables to be specified to uniquely determine these polynomials. To decide these variables, we need to select r collocation points in each partition; in total N = Np 1). In our model, we consider each message passing layers prediction time step collocation point Ne k= which means each r- 1 rollout prediction belongs to one partition, i.e. 4p.n < (c,a(r-1) < (c((+))((r-1) < ( nl = 0,...,.. ,N,
select Np partitions in the time domain 0 00= 'p.u < (pp11 <
(le,
N. By substituting the state of these collocation points fu," k=o to polynomials, we can transfer this simulation problem an algebraic equation. We emphasize that the coefficient matrix of this alge- braic equation is almost block diugonal (ABD) (De Boor & De Boor, 1978). This kind of system allows for efficient computational routines (Amodio et al., 2000), that we will introduce By solving the equation, value for this sample point at any time أ²² At e [0, 2). All these steps can be calculated by matrix operations, so that the model is fully end-to-end differentiable.
$ 3.5. obtain the simulation polynomial (-(-) that can be used to simulate the
Space-oriented OSC For a specific time frame t.. the states of sample points in the physical domain (u!* [ Ji=o be described by PDEs f("* (x,)) 0,æ; E Ω, where 2²* (-) is the solution for the state at this time frame. The target of the space-oriented OSC is to find one polynomial of order each partition for every dimension and make these N D polynomials C2 continuous in the
domain.
Similar with the time-oriented OSC, we select Np.a partitions in the dth dimension Na Biof a Pa.o < Pd. < ..... < Pd,N,.a Bsup.d, d = 1, ... D, where Bint,d, Baup.d are the lower and upper boundary of dimension d. Note that these partitions can be isometric. Then we initialize one polynomial under order in each partition for each dimension. These polynomials have the degree of freedom (r))r in total. To determine these polynomials, we need select - collocation points in each partition, in total N. = NP P) In our model, we consider the prediction states of sample points as collocation points {x1,* [
4

---

Under review as a conference paper at ICLR 2023
1.N. By substituting the state of these collocation points {u, T== to polynomials, simulation problem to an algebraic problem with an ABD coefficient matrix. By solving this prob- lem we obtain the simulation polynomial â'*, that can be used for simulating the result in the whole domain 8 (*(()), (x), vx E Ω. More details about deriving the degree of freedom, vi- sualization of the ABD coefficient matrix, and further techniques of applying OSC shown in Appendix A.
transfer this
3.5
ALGORITHM FOR EFFICIENTLY SOLVING THE ABD MATRIX
Most interpolation methods need to solve linear equa- tions. Gaussian elimination is one of the most widely used methods to solve a (dense) linear equation, which however has ( complexity (Strassen 1969). Even with the best algorithms known to date, the lower bound of time complexity to solve such equations is (on))) (Golub & Van Loan, 2013). In the OSC method, the coefficient matrix for the lin- ear equation follows the ABD structure, which we can efficiently solve with a time complexity Onn²) by the COLROW algorithm (Diaz et al., 1983) as shown in The core idea for this method is that by using the pivotal strategy and elimination multipliers, we can decompose the coefficient matrix into a set of permu- tation matrix and upper lower triangular matrix that can be solved in (n²) time each. The most recent package providing this algorithm is in the FORTRAN programming language: our re- implementation in PyTorch (Paszke et al., 2019) allows for optimized calculation, GPU support and enabling the use of automatic differentiation.
al.,
Time
10-
0-
OSC
OSC
OSC ROW Solver
Inverse Solver
COLROW
Solver
Fig. 2.
Figure 2: The COLROW algorithm allows for fast solutions the OSC problem. GPU imple- mentation further improves on inference times.
3.6 ADAPTIVE COLLOCATION SAMPLING
To allow for prioritized sampling of important locations regions, we optimize the positions of collocation points via gradient descent on the data of states history, projecting them back to the partition if the gradient step moves them outside to make sure there is a sufficient num- ber of collocation point in each partition cell. We calculate the gradient of each col- location point along each dimension from the continuity predictions and then use this gra- dient vector optimize the position of the point by weighted sum illustrated in Fig. 3. We use the states at optimized positions adapted from history rollouts as the next rollouts input. By adapt- ing the collocation points position our model is able to place greater emphasis harder parts of the space to get a more accurate prediction.
Static
Adaptive
3.7 TRAINING STRATEGY AND LOSS FUNCTION
To train the model, the collocation points at the initial- ized state are input to the MPNN to propagate roll- outs on these points as shown in Fig. 4. Then two OSC methods allow for apatio-temporal continuous outputs, so that we can use any states in the predic- tion range as target to train our model end-to-end. More specifically, given an input (X, Y"), the model outputs the simulated polynomial ú(x, t), x Ω, t [0, T) for T seconds prediction. By making use of higher resolution sample points along time and space in the train set as the target I=NN. where N, is the number of target sample point and N, is the number of target
Figure 3: Adaptive collocation strategy: mesh points converge towards areas with higher infor- mation density.
(xi, ytt )==0,k=0 sample time frames, we can calculate the reconstruction loss of sample points Lg. Moreover, we also calculate the reconstruction loss of predicted collocation points Tec. The complete loss is
5

---

Under review as a conference paper at ICLR 2023
N₂
No l
N,
N1 1
L
ú(x;, 12
(2)
L.
points reconstruction
semple points reconstruction
where the N₂ is the number of MPNN rollout steps. The whole model is then trained end-to-end with automatic differentiation through the OSC.
EXPERIMENTS
4.1 DATASETS AND TRAINING
We test our GRAPHSPLINENETS on three settings of increasing challenge:
to
tk
N
Heat Equation: this PDE describes the behavior of heat through a domain and characterized by diffusivity leading to stable solutions over time.
Damped Wave: this system describes the propaga- tion of waves through medium characterized by their velocity; an additional damping term smooths time.
their amplitude
Collocation Points
Sample Points
Navier-Stokes:
these equations describe the be- havior of 2D incompressible fluid flows; turbu- lence leading to chaoticity makes them notoriously hard problem to solve.
Figure 4: Data points used for model training.
All the models we test employ the same structure of encoderprroessoo-dddoeer for fair comparisons well the same amount of training data in each testing domain. While inputs of baseline models are directly all the available data points, inputs of our OSC-based models are only an initialized 12 12 collocation point unitary mesh at the initial state and fewer in-time sample points.
4.2 EVALUATION METRICS AND BASELINES
We evaluated our model by calculating the average mean square error (MSE) of 1(s), 2(s), 5s)) rollout predictions steps with respectively ground truth. We employ baselines relevant works in the field of discrete-step graph models for dynamical system predictions. Graph convolution networks (GCNs) (Kipf & Welling, 2016) and GCN with a hybrid multilayer perceptron (MLP) model employed as baselines in the ablation study. We also compare our approach with one widely used baseline that employs linear interpolation for physics simulations allowing for continuous predic- tions in space, i.e., GEN (Alet et al., 2019). A similar setup considering the inherent graph structure and employing message passing neural networks in mesh space is employed by (Pfaff et al., 2021) (MPNN in our table). We employ the latter as the first building block for our GRAPHSPLINENETS.
4.3 QUANTITATIVE ANALYSIS
Empirical quantitative results on the three dataset are shown Table In the heat equation dataset, approach reduces long range prediction errors by 64% error while only using 20% of the running time compared with the best baseline model. In the damped wave dataset, our approach reduces emror by 42% with 48%. In the Navier-Stokes dataset, our approach reduces 31% long-range prediction errors while requiring 37% less time to infer solutions compared to the strongest baseline.
4.4 ABLATION STUDY
We consider an ablation study on three datasets demonstrate the effectiveness of our model com- ponents in our approach in multiple aspects. Quantitative results of ablation study models are shown Table [right]. The ablated models are:
MPNN: base graph model with 3 message passing layers.

---

Under review as a conference paper at ICLR 2023
Table 1: Mean square error (MSE) propagation at different time stamps in seconds. Runtimes consider model inference for the full rollouts. Smaller is better (J). Best in bold; second underlined.
Disaset
Metric
Duselines
GRAPESPLINENETS
GCN
GUN+MLP
GEN
MENN
MPNN+OSCIPPOI
MPNN+OSC
MPNN+OSC+AAdpph
166)
11.2
0.200
0:48 = 0.05
0.28
0.02
0.37 = 1.05
0,38 2: 2:04
0.38 a 0.00
11.2 =: 0.00
266
1.05 0.12
0.87 + 0.06
0.53
0.06
0.60 + 0.06
0.68 +..0-
0.56 + 0.00
0.46+-0. +.00+
(XC00
3ts)
2.58
2.05
1 0.10
1.55
0:09
1.85 L 0.13
1.28 1.0.0
0.9 1 1.008
0.87 0.00
Healt Equition
4.12 0.44
8.87
0.23
= 0.36
2.09
+.111
1.02
0.38
6.87
0.25
3.01
1.14 0.31
0.22
Rollout
50
50
12
12
3.26
1:61
1.41
0.81
0.00
0.27
1.15
1.09
0.14
(X10 ))
5.12
0.80
0.25
0.27 =: 0:23
2.57 0.18
2,48
Damped Wave
7.77
6.01
0.82
0:41
6.27 D.1
5.88 0:25
3:41 0.22
0.25
56s)
10.5
0.00
1.52
0.00
0.62
7.82 0.88
4.98 0.20
0.27 4.60
4.51 0.31
Rollost stops Runime [s]
10
10
0.05 0.08
0.82 1: 0.07
1.15
0.00
1.38 1L 0.10
0:45 0.05
0.39 1 0.04
0.42 L: 0.00
16s)
1.47 0.10
1.22
0.11
0.42 0.007
0.66 L. 0.10
0.72 0.00
0.70 0.08
0.54 0:09
266)
2.01 + 0:2.2
1.76 + 0.20
0.30 0.08%
1.18 11:11
1:20 + 0.11
1,022 0.11
0.80 0.14
MSE (x (x:10-1)
3(5)
2.81 +: 0.33
2.45 + 0.96
1.05
+:5==000 0.12
1.64 0.24
1:66 0.18
1.44 0.20
1.23 0.20
Navier Stokes
4is)
1.51
1: 0.64
2.94 a: 0.82
2.38
d: 0.16
2.57 L: 0.28
1.98 1: 0.24
1.72 a 0.27
1.30 0.29
5(6)
4.24 + 0.95
3.91 + 0.30
3.46
0.24
3.06 = 0.53
2:58 + 0.28
2.21 + 0.27
2.02 0.30
Rollowt steps Ruttime
10
10
10
10
5 0.51
5
5
0.31 +
0.88 0.07
1.01
101+0.000 0.00
1.21 = 0.10
+.0.5
0.47 + 0.04
0.40 0.00
(PNN+ model with 3 message passing layers and only post-processing with the OSC method, i.e., we firstly train a MPNN model, then we use the OSC method to collocate the prediction as final result without end-to-end training.
MPNN++SC: MPNN with OSC- OSC-in-the-loop that allows for end-to-end training.
MPNN+OSC+AAaatiiv: MPNN+OSC model that additionally employs our adaptive col-
location point sampling strategy.
Interpolation and collocation methods We demon- strate the efficiency of the OSC method by compar- ing the combination of MPNN with different interpo- lation and collocation methods, including linear inter- polation, cubic interpolation, and B-spline collocation methods. These models are implemented in the end- to-end training loop and we use these methods in both the time and space dimensions. Experiment results are shown in Fig. 5 where we measured the mean square error and the running time of 3 second rollouts pre- dictions. We also test each method with different number of collocation points, i.c. from (2 x 2) to (16 x 16) in Fig. 6. The model MPNN+OSC shows the best performance obtaining the highest accuracy prediction among these approaches with shorter run- ning time. Even though the linear interpolation can be faster than the OSC, it shows a considerable error n the prediction and does not satisfy basic assumptions such as Lipschitz continuity in space.
10
10-1
MSE
2
10-1
10
10-1
10"
Prediction Running Time [s]
MPNN
Linear
MPNN B Spline
MPNN-
Cubic
MPNN OSC
Figure 5: Pareto front different interpolation and collocation points in terms of accuracy and speed with MPNNs elearly exhibits the advan- tages of the OSC method.
Number of collocation points We study the effect of the number of collocation points on the sec- ond rollout prediction eтror by testing the MPNN, MPNN+OSC and MPNN+OSC+Adapin mod- els. The MPNN will always be directly trained with the whole domain data. We use different number of collocation points (from (2 2) to (28 x 28)) into the MPNN process in the rest two models and then compare the output of the OSC with the whole domain to train. With the increase in num- ber of collocation points, Fig. 6 shows that the MPNN+OSC and MPNN+OSC+Adapiiv achieve significant improvements in prediction accuracy over the MPNN. The MPNN+OSC+Adaptive has stable better performance compare with the MPNN+OSC and the improvement is larger when there are fewer collocation points. The reason is the with fewer collocation points, the MPNN+OSC has insufficient ability to learn the whole domain. With the help of the adaptive collocation point,

---

Under review as a conference paper at ICLR 2023
the MPNN+OSC+Adaativee can focus on hard-to-learning regions during training obtain overall better predictions.
COLROW solver and GPU acceleration We show the effectiveness of the COLROW algorithm in accel- erating the OSC speed by comparing the OSC method with one the most commonly used algorithms for efficiently solving linear systems² and the OSC with the COLROW solver. Fig. shows the experimental results 1000 solutions of ABD linear equations with the size of 256 x 256. Our package speeds up the OSC method by 32%. By making use of the GPU to speed up the OSC simulation, we can further lower the run- ning time by 60%.
10-14
MSE
10
[8 8)
[16 x 16)
(34 x 14)
Number of Collocation Points
MPNN
Adaptive
MPNN OSC
Number of rollout steps We show the ef-
fectiveness of the OSC method improving long-range prediction accuracy by comparing the MPNN and MPNN+OSC model. Fig. shows the MPNN+OSC can keep stable in long-range rollouts compare with the MPNN. The reason is that with the OSC, we can use fewer neural network rollout steps to obtain a longer- range predictions, which avoids the error accumulation during the multi-step rollouts and implicitly learns for compensating integration residual errors. End-to-end learning lets the neural networks in MPNN+OSC learn the states between rollout steps, which will make the prediction stable and accurate.
Figure 6: Number of collocation points and mean squared error.
Post processing vs end-to-end learning We show the effectiveness of end-to-end learning architecture by comparing the MPNN+OSC(Post) and MPNN+OSC models. Table 1 show that MPNN+OSC has a more accurate prediction than MPNN+OSC(Post) by more than 8% percent across datasets. This can be ex- plained by the fact that, since the OSC applied end- to-end, the error between MPNN prediction steps is backpropagated to the message passing layers, while in the post processing steps the model has no way of considering such error.
3-
MSE
2-
23
40
Rollout Steps
MPNN
MPNN OSC
Adaptive collocation point We further show the effectiveness of the adaptive collocation strategy by comparing the MPNN+OSC and MPNN+OSC+
Figure 7: MPNN+OSC contains long rollout emrors better than standard MPNN.
Fig. 6 shows that the MPNN+OSC+Adaptive has a better performance than the MPNN+OSC in all collocation point setups. And Table 1 shows that MPNN+OSC+Adaptiv has a more accurate prediction than MPNN+OSC, i.e. around 10% improvement on long rollouts the Navier-Stokes dataset. Adaptive collocation points encourage those points to move to the most dynamic regions in the domain, which is not only able to place greater emphasis on hard-to-learn parts in space, but can let the OSC method develop better implicit representation of the domain.
4.5 QUALITATIVE ANALYSIS
We visualize the damped wave equation propagation and the Navier-Stokes evolution results of the GCN, GEN and MPNN+OSC(Our) eтror distribution and more stable long-range prediction. Thanks to GRAPHSPLINENETS continuous predictions, we can simulate high resolutions without needing additional expensive model inference routines, while the other two models can only achieve lower resolution predictions. For long-range
Fig. 8, Fig. and Fig. 10. Our model has smoother
We used for our experiments torch. forch.linale 1g. solve, which uses LU decomposition with partial piv- oting and row interchanges. This is faster and numerically stable than matrix inversion; however, it has still (n time complexity.
8

---

Under review as a conference paper at ICLR 2023
Gmann Trush
GCN
GFN
Colleration Points
GCN
GEM Our Modell
Figure 8: The visualization of results and the enror the wave dataset. Black dots in the collocation points figure the position of sample points for models training.
t-3
t-44
5
t-
t=7
t-8
Ground Truth
GCN
0.6
GEN
0.4
0.2
Our
Figure 9: Wave dataset prediction results. GRAPHSPLINENETS manages to obtain more stable and smoother predictions compared to baselines.
predictions, while baselines visibly accumulate error, our model can lower the error with smoother and more accurate predictions in space and time.
CONCLUSION
We introduce GRAPIISPLINENETS, novel method that aims briding the gap between inherent discrete graph predictions in space and time and the continuous essence of natural processes. Our approach integrates the theory of Orthogonal Spline Collocation (OSC) methods achieve space and time continuous simulations without the need for computationally expensive numerical routines. We introduce an effective adaptive collocation strategy to prioritize sampling of points the space domain and implement the OSC end-to-end for achieving continuous predictions. We demonstrate how GRAPHSPLINENETS are robust in predicting processes which are characterized by several different PDEs arising directly from the differential equations. We believe this work represents an important step forward new direction the research area at the intersection deep learning and dynamical systems that aims at finding fast and accurate learned surrogate models.
Ground Truch
GCN
GEN
Our
Collocation Points
CCN Errrr
CEN Ernrr
Our Model Erro
3
Figure 10: Qualitative results emoor on the Navier- Navier-Stokes dataset. Black in the collocation points figure are the positions of sample points for our models training.

---

Under review as a conference paper at ICLR 2023
REFERENCES
Ferran Alet, Adarsh Keshav Jeewajee, Maria Bauza Villalonga, Alberto Rodriguez, Tomas Lozano- Perez, and Leslie Kaelbling. Graph element networks: adaptive, structured computation and memory. In International Conference Machine Learning, pp. 212-222. PMLR, 2019.
Pierluigi Amodio, JR Cash, G Roussos, RW Wright, Graeme Fairweather, lan Gladwell. GL Kraut, and Marcin Paprzycki. Almost block diagonal linear systems: sequential and parallel solution techniques, and applications. Numerical linear algebra with applications, 7(5):275-317, 2000.
Peter Bauer, Alan Thorpe. and Gilbert Brunet. The quiet revolution of numerical weather prediction. Nature, 525(7567):47-55, 2015.
Federico Berto, Stefano Massaroli, Michael Poli, and Jinkyoo Park. Neural solvers for fast and ac- curate numerical optimal control. International Conferences on Learning Representations, 2022.
B. Bialecki and G. Fairweather. Orthogonal spline collocation methods for partial differential equations. Journal of Computational and Applied Mathematics, 128(1):55-82, 2001. ISSN 0377-0427. doi: https://doi.org/ http://dbi.org100011//////00 00000/00160007770200000002000000 sciencedirect com/ wttp//ciecncceeeeeeeeeeeeeeeeeee /ttp:///rr..eee///////7770000000 Analysis 2000, Vol. VII: Partial Differential Equations.
URL https::/weet Numerical
Bernard Bialecki. Convergence analysis of orthogonal spline collocation for elliptic boundary value problems. SIAM journal on numerical analysis, 35(2):617-631, 1998.
Paul Bourke. Interpolation methods. Miscellaneous: projection, modelling, rendering, 1(10), 1999.
Johannes Brandstetter, Daniel Worrall, and Max Welling. Message passing neural pde solvers. International Conference on Learning Representation, 2022.
Adam R Brink, David A Najera-Flors and Cari Martinez. The neural network collocation method for solving partial differential equations. Neural Computing and Applications, 33(1 (33(11):591111 15551155608, 2021.
Michael M Bronstein, Joan Bruna, Taco Cohen, and Petar Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2 13478, 2021.
Carl De Boor and Carl De Boor. A practical guide to splines, volume 27. springer-verlag New York, 1978.
Carl de Boor and Richard Weiss. Solveblok: A package for solving almost block diagonal linear systems. ACM Transactions on Mathematical Software (TOMS), 6(1):80-87, 1980.
JC Diaz, Graeme Fairweather, and Patrick Keast. Fortran packages for solving certain almost block diagonal linear systems by modified alternate row and column elimination. ACM Transactions on Mathematical Software (TOMS), 93)):558-377 1983.
Graeme Fairweather and Daniel Meade. A survey of spline collocation methods for the numerical solution of differential equations. In Mathematics for large scale computing, pp. 297-341. CRC
Press, 2020.
William Falcon et al. Pytorch lightning. compPyyooolLiaannns/ooohhhiiiiie 3(6), 2019.
GitHub.
Note:
https://github.
Zheng Fang, Qingqing Long, Guojie Song, and Kunqing Xie. Spatial-temporal graph ode networks for traffic flow forecasting. In Proceedings of the 27th ACM SIGKDD Conference on Knowledge Discovery & Data Mining, pp. 364-373, 2021.
Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013.
Hongwei Guo, Xiaoying Zhuang, and Timon Rabezuk. A deep collocation method for the bending
analysis of kirchhoff plate. Computers, Materials & Continua, 59(2):433-455 456, 2019. ISSN 1546-2226, doi: 20w.040/mc..00.00000.m URL http: /xx, doi. ooo.oorg/100 /02.0//...0.000.... cmc. 201 9. 06660.
10

---

Under review as a conference paper at ICLR 2023
Boris Houska, Filip Logist, Moritz Diehl, and Jan Van Impe. A tutorial numerical methods for state and parameter estimation nonlinear dynamic systems. Identification for Automotive Systems, pp. 67-88, 2012.
Ryan Keisler. Forecasting global weather with graph neural networks. arXiv:2202.07575. 2022.
arXiv preprint
Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional net- 2016.
works. arXiv preprint
Matthieu Kirchmeyer, Yuan Yin, Jérémie Donå, Nicolas Baskiotis, Alain Rakotomamonjy, and Patrick Gallinari. Generalizing new physical systems via context-informed dynamics model. arXiv preprint arXiv:2202.0188 2022.
Jiachen Li, Chuanbo Hua, Jinkyoo Park, Hengbo Ma. Victoria Dax, and Mykel J Kochenderfer. Evolvehypergraph: Group-aware dynamic relational reasoning for trajectory prediction. arXiv preprint arXiv:2208.054422 2022a.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, An- drew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint arXiv:2010.08895, 2020.
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. arXiv preprint arXiv:2207.05209, 2022b.
Marten Lienen and Stephan Günnemann. Learning the dynamics of physical systems from sparse observations with finite element networks. International Conference on Learning Representation, 2022.
Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012.
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong. Pde-net: Learning pdes from data. In International Conference Machine Learning, pp. 2216-3216. PMLR, 2018.
Pru Marriott, Siew Min Tan, and Neil Marriott. Experiential learning-a a case study of the use of computerised stock market trading simulation in finance education. Accounting Education, 24 (6):480-499 2015.
Damian Mrowca, Chengxu Zhuang, Elias Wang, Nick Haber, Li F Fei-Fei, Josh Tenenbaum, and Daniel Yamins. Flexible neural representation for physics prediction. Advances in neural information processing systems, 31, 2018.
Anita Mücke, Ralph Engel, Jörg Rachen, Raymond Protheroe, and Todor Stanev. Monte carlo simulations of photohadronic processes in astrophysics. Computer Physics Communications, 124 (2-3):290-314, 2000.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high- performance deep learning library. Advances neural information processing systems, 32, 2019.
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay. Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al. Four- castnet: A global data-driven high-resolution weather model using adaptive fourier neural opera- tors. arXiv preprint arXiv: 12214, 2022.
Tobias Pfaff, Meire Fortunato, Alvaro Gonzllzz and Peter W Battaglia. Learning mesh- based simulation with graph networks. International Conference on Learning Representations, 2021.
11

---

Under review as a conference paper at ICLR 2023
Michael Poli, Stefano Massaroli, Junyoung Park, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Graph neural ordinary differential equations. arXiv preprint arXiv:1911.07532, 2019.
Michael Poli, Stefano Massaroli, Atsushi Yamashita, Hajime Asama, and Jinkyoo Park. Hyper- solvers: Toward fast continuous-depth models. Advances in Neural Information Processing Sys- tems, 33:21 105-21117, 2020.
Michael Poli, Stefano Massaroli, Federico Berto, Jinkyoo Park, Tri Dao, Christopher Re, and Ste- fano Ermon. Transform once: Efficient operator learning in frequency domain. 2022.
Leijie Qiao, Da Xu, and Zhibo Wang. Orthogonal spline collocation method for the -wo-dimensional time fractional mobile-immobile equation. Journal of Applied Mathematics and Computing, pp- 1-19, 2021.
Alfio Quarteroni and Alessandro Veneziani. Analysis of a geometrical multiscale model based on the coupling of ode and pde for blood flow simulations. Multiscale Modeling Simulation, 1(2): 173-195, 2003.
Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686-707, 2019.
Alvaro Sanchez-Gonzaler Jonathan Godwin, Tobias Pfaff, Rex Ying, Jure Leskovec, and Peter Battaglia. Learning to simulate complex physics with graph networks. In International Confer- ence on Machine Learning, pp. 8459-8468. PMLR, 2020.
Volker Strassen al. Gaussian elimination not optimal. Numerische mathematik, 13(4):354-356, 1969.
Minjie Wang, Da Zheng, Zihao Ye, Quan Gan, Mufei Li, Xiang Song, Jinjing Zhou, Chao Ma, Lingfan Yu, Gai, Tianjun Xiao, Tong He, George Karypis, Jinyang Li, and Zheng Zhang. Deep graph library: A graph-centric,
package for graph neural networks, 2020.
Yu-Ting Wu and Fernando Porté-Agel. Large-eddy simulation of wind-turbine wakes: evaluation of turbine parametrisations. Boundary-layer meteorology. 138(3):345-366, 2011.
Louis-Pascal Xhonneux, Meng Qu, and Jian Tang. Continuous graph neural networks. In Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learn- ing. volume 119 of Proceedings of Machine Learning Research, pp. 10432-10441, PMLR, 13-18 Jul 2020. URL wwt.ttttp//cccc httpss:/pproceedinggg
www.mmrrreess....m.mm press / v119/ /tt11//1h11/////2//./.. wwwnneux20a.....mmmmmoo......... html.
Yuan Yin, Vincent Le Guen, Jérémic Dona, Emmanuel de Bézenac, Ibrahim Ayed, Nicolas Thome, and Patrick Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment, 2021(12):124012, 2021.
12

---

Under review as a conference paper at ICLR 2023
A ADDITIONAL OSC MATERIAL
We further illustrate the OSC method by providing numerical examples in this section.
A.1 1-D OSC EXAMPLE
For simplicity and without loss of generality, we consider the function domain unit domain [0, 1] and we set N 3, = 2. which means we will use a threc-order three-piece function to simulate the I-D ODE problem as shown in Equation. 2. We firstly choose the partition points 24, 0,... 3, To = 0, 33 1. The number of partition points is N + partition points can fixed or not fixed. Then we based on Gauss- Legendre quadrature rule choose collocation points. The number of collocation point within one partition is r -1 1 = 1, so we have in total x (r - 1) collocation points E,, i = 0, ,3.
4. Distance between
After getting partition points and collocation points, we will construct the simulator. Here we have three partitions, in each partition, we assign a order polynomial
00,0 00 + a0,20² x
[20 x1]
(3a)
01,0 + 01,1 ,1,2
[21, 322
(3b)
02,0 02, + 022 x [22 T3
(3c)
Notice that these three polynomials should be C' continuous the connecting points, i.e. parti-
tion points within the domain. For example, Equation Eq. (3a) and Equation Eq. (3b) should be continuous 2i, then we can get two equations
00.00 + 00,1 + 00,201
01,0 + @1,1211 + 01,20r
(4)
00,1 + 200
+ 01,1 + 201,2
For boundary condition
u(x) = (bi I = 30
(5)
12, x 23
we can also get two equations
00.0 + 0 0
b
(6)
a11o + 01,1 + 01,2
= b2
Sum up equations get far. Firstly our undefined polynomials have N (r+1) = 9 parameters. The C² continuous condition will create (N - 1) x = 4 equations and the boundary condition will create 2 equations. Then have N x (r 1) collocation points. For each collocation point, we substitute it to polynomials to get an equation. For example, if the ODE is
a(x) + (x) f(z). x [0, 1]
By substituting collocation point &o into the equation, we can get
((Eo) + úl (Eo) = f(Eo)
40,0 + Co + 2 + 40,1 + 200 = f(Eo)
(8)
40.0 + 1 + 1) + 2 + 25o) = f(Eo)
Now we can know that the number of equations can meet with the degree of freedom of polynomials
Parameters
Boundary
C¹ Contiimouss
Collocation
+
2
(N 1) N)
At this example, generated equations will be constructed to an algebra problem Aa = f where the weight matrix is ABD matrix.
13

---

Under review as a conference paper at ICLR 2023
Not
Figure 11: Visualization of an ABD matrix.
1
0
0
0
0
Eo
+
1
+
250
0
0
0
A
-1
(10a)
2x
O
-1
-2
0
E1
+1
E
2
0
0
1
1
d0.0
bi
00,1
f(Eu)
00,2
f
(10b)
01,0
01,1
f(G)
01,2
b2
To solve this problem, we can get the simulation results.
A.2 2-D OSC EXAMPLE
For simplicity and without loss of generality, consider the function domain unit domain [0, x [0, and we set N₂, = N, 2,r 3. Partition points and collocation points selection are similar with 1-D OSC method, we have N² (r 1)² 16 collocation points in total. For simplicity, we note the partition points at two dimensions to be the same, i.e. zi, i = 0, 1,2. Unlike the 1-D OSC method, we choose Hermit bases to deseribe as the simulator, which keeps C continuous. As case, the base function at point 21 would be
H1(x) = (x) + (x)
(2-2a)(2; =x)²
(21 (#1-20)² E (Jon v1]
f(())= = (x-22)(z -1 E (I1, F2
(#2-21)²
(11)
(2) (2)) 2(21 - 2n)] (ar so--aao x
E (ro, m1]
9: (x =
(2) xo
(2g 2(1 [1ng x)²
E (I1,
(0--21))
We separately assign parameters to basis functions, i.e. H1 (x) = (x)) (x)++b. b1.99h (1:) for variable
(ro, 21] [ Vil partition. Then the polynomial in a partition is the multiple combinations of bases functions of two dimensions. For example, the polynomial in the partition [zo, x1 x [Yo- V1]
fo(x) 16 + (z) (x)) x fo(v) bg (v) (v + (v)]
(12)
Now we consider the freedom degree of these polynomials. From definition, we have (n 1n 1) 24 parameter. Consider boundary conditions, we have 24 4 x N 16 parameters. The number is equal with collocation points 2 x (r ()2,, which means we can get an algebra equation by substituting collocation points. Solving this equation can get the simulator parameters.
We can similarly multiple basis functions and set parameters to the simulation result for the higher dimension OSC method. And then select partition points and collocation points by the same strategy with 2-D OSC method. The rest algebra equation generating and solving equations parts will not be different.
14

---

Under review as a conference paper at ICLR 2023
A.3 SIMPLE NUMERICAL EXAMPLE
We set N = 3, r = 3 to simulate the problem
{
u v' sin() 2500os u(0) = 0 u(1) = 0
(13)
we can get simulation solution as following, with is visualized in Fig. 12.
6.22 0.4x² - 31.4x x E [0,1/3)
a(x)=
1.5 + 1.6x 13.8.² + +9x2, x [1/3, 2/3)
(14)
28.5 100x 108.5x² - 37x E [2/3,
0.5
0.0
0.5
-1.0
0.0
0.2
0.4
0.6
0.8
1.0
Simulation
Ground Truth
Partition Points
Collocation Points
Figure 12: Visualization of OSC solution.
A.4 NUMERICAL ANALYSIS FOR INTERPOLATION AND COLLOCATION METHODS
We compared the OSC with linear, bilinear, O-D cubic, and 2-D cubic interpolation methods on four types problem: 1-D linear, 1-D non-linear, 2-D linear, and 2-D non-linear problems. At these experiments, tested different simulator orders of the OSC method. For example, we set the order of the simulator to 4 for 1-D linear problem and for 2-D linear problem. When the order of the simulator matches the polynomial order of the real solution, OSC directly find the real solution. For non-linear problems, increasing the order of the simulator would be an ideal way to get lower loss. For example, we the order of the simulator to 4 for 1-D non-linear problem and for 2-D non-linear problem. Thanks to the efficient calculation of OSC, even though we use higher-order polynomials simulate, we use less running time to get results.
Table 2: Error of OSC and four interpolation methods different PDEs problems: u(x) x 3x + 1.16x 0.16x (1-D linear), u(x) = (3 (1-D non-linear), u(z, y) = x² 2 - x y xy² + zy (2-D linear), u(æ, y) = sin (3 (inn)nnnnnnn))) (2-D non non-linear).
MODEL
1.D LINEAR
I.D NON LINEAR
2.D LINEAR
2.D
NEAREST
2.3670x 10
1.7558x10
1.9882x 10-3
3.8695 10
LINEAR INTERPOLATION
1.8928x 10
8.7731 10
3,4317x 10
1.1934x 10
CUDIC INTERFOLATION
3.5232x 10
2,2654x 10
2.9117x10
4,5441 10
OSC
3,4153x10
4,1948x10
1.723xx10
3.4462x10
B ADDITIONAL EXPERIMENTAL DETAILS
HEAT EQUATION
The heat equation describes the diffusive process of heat conveyance and be defined by
15

---

Under review as a conference paper at ICLR 2023
Du
Au
(15)
where u denotes the solution to the equation and the Laplacian operator over the domain. In n-dimenniinal space, it can be written as:
Ou u dz.
Au
(16)
Dataset generation We employ FEniCS (Logg et al., 2012) to generate a mesh from the domain and solve the heat equation on these points. The mesh then used by the graph neural network for training.
B.2
WAVE EQUATION
The damped wave equation can be defined by
82w dw
Aw
ot2
Öt
where e is the wave speed and & is the damping coefficient. The state is X (w,
Data generation We consider spatial domain Ω represented 64 x 64 grid and discretize the Laplacian operator. A is implemented using a 5 x discrete Laplace operator in simulation; null Neumann boundary condition are imposed for generation. We 330 and k 50 similarly the original implementation in Yin et al. (2021).
B.3 NAVIER-STOKEE
The Navier-Stoker equations describe the dynamics of incompressible flows with 2-dimensional PDE. They can be described in vorticity form as:
dw
ow vA +
Øt
(17)
Ve
0
= V X
where v is the velocity field and w is the vorticity, is the viscosity and f is a forcing term. The domain is subject to periodic boundary conditions.
Data generation We generate trajectories with a temporal resolution At and time horizon of = 10. We use similar settings as in Yin et al. (2021) and Kirchmeyer et al. (2022): the space is discretized on a 64 x grid and we set [x, y) 0.1 (sin((2x< + v)) + cos(2 + y))), where are coordinates on the discretized domain. We use viscosity value = 10-3
B.4 MODELS AND IMPLEMENTION DETAILS
For all experiments, batch size of 32 used and the models trained for up to 5000 epoches with early stopping. We used the Adam optimizer (Kingma & Ba, 2014) with an initial learning rate of 0.001 and step scheduler with 0.85 decay rate every 500 epochs. For all datasets, used the split of 5 for training, validating and testing for fair comparison.
Specific details of model components introduced below:
MPNN encoder: three-layer MLP with hidden size= 64.
MPNN processor: in total 3 processors, each with three-layer MLP with hidden size= 64.
MPNN decoder; three-layer MI.P with hidden size= 64.
16

---

Under review as a conference paper at ICLR 2023
All MLPs have ReLU x max(0, z) nonlinearities between layers.
Specific details of applying the OSC method are introduced as follows:
Time-oriented OSC: polynomial order= 3, number of collocation points in one partition=
2.
Space-oriented OSC: x dimension polynomial order= 3. polynomial order in the dimension= 3, number of collocation point in partition dimension= 2; and num- ber of collocation point in one partition on y dimension= 2.
B.5 HARDWARE AND SOFTWARE
Experiments were carried out on a machine equipped with an INTEL CORE 19 7900X CPU with 20 threads and NVIDIA RTX A5000 graphic card with 24 GB of VRAM. Software- we used FEniCS (Logg et al., 2012) for Finite Element simulations for the heat equation experiments and PyTorch (Paszke et al., 2019) for simulating the damped wave and Navier-Stokes equations. We employed the Deep Graph Library (DGL) (Wang et al., 2020) for graph neural networks and the PyTorch Lightning library (Falcon et al., 2019) for training.
Ground Truth
GCN
GCP-
GEN
MPNN
GNN OSC(Ot)(
HNCC
GNN OSC(Adaptive)
Figure 13: Visualization of the dataset.
Ground Truth
GCN
GCN-ME MLP
GEN
MPNN
CNN-O5C(PPsE)
GNN-OSC
(NNddS((A
Figure 14: Visualization of the Navier Stokes dataset.
17


Reviews:
[{"review_id": "LITh_XU2-V", "reviewer": "Reviewer_ESAs", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposes to combine message-passing NNs with the orthogonal spline collocation (OSC) to build a method capable of predicting PDE solutions that are continuous in time and space and which can be efficiently evaluated to generate a value at arbitrary points in space and time within the simulated domain. The method is tested on the heat equation, damped wave, and 2d Navier-Stokes (incompressible fluid flow), and shown to be more efficient than classical interpolation schemes, to produce smoother solutions than baselines, and generate more stable long-range predictions.\n", "summary_of_the_review": "The paper introduces an interesting combination of existing techniques, with experiments illustrating gains in accuracy and performance in a number of simple settings. This is interesting, but would be strengthened by more extensive experiments, or by more focus on reproducibility (ideally in the form of open source code).", "strength_and_weaknesses": "The strengths of the paper include the experimental results showing the utility of the proposed method, and the inclusion of ablations illustrating the impact of OSC as a post-processing step, or when trained end-to-end with and without adaptive colocation points. It is also very helpful that the authors include error images in Figs 8 and 10.\n\nA limitation is that all the presented testcases are ones with relatively simple structure and dynamics. It would be interesting how the method performs in a more complex situation. For instance, in the 2D NSE case, one could consider analyzing a Karman vortex street or unstable flow at a higher Reynolds numbers. There is also little novelty on the ML front -- the proposed method is more of a (trainable) postprocessing step for GNN-based learned surrogate models.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "KlxW6BEjzB1", "reviewer": "Reviewer_5HfZ", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper combines graph network simulators with orthogonal spline collocation (OSC) methods to model complex dynamic systems at continuous spatial and temporal locations. The technical method first applies standard graph network models to predict a dynamic system’s state at discrete spatial locations and time steps. Next, it applies OSC to augment the discrete results into continuous solutions. The whole model is end-to-end differentiable so that standard optimization algorithms can apply. The paper presents evaluation experiments on a few common PDE problems, e.g., heat equations and Navier-Stokes equations.", "summary_of_the_review": "My major concern is mostly about the technical approach itself. While I appreciate the paper’s effort, I don’t think the proposed approach is significant or novel enough to reach the ICLR standard. However, I am open to different opinions from the authors or other reviewers.", "strength_and_weaknesses": "The paper is a respectable submission, but I am not sure I am fully convinced by the proposed technical method. I have listed my questions below, and it would be good if the authors could address them during the rebuttal:\n\nThe whole method reminds me of the Physics-Informed Neural Network (PINN) [Raissi et al.], which also attempts to reconstruct a continuous spatial-temporal solution to PDE problems. The related work section seems to overlook PINN, but I think it deserves a more careful discussion to clarify the difference between this paper and PINNs.\n\nActually, I am not sure using OSC to construct continuous spatial-temporal solutions is as significant as the paper claims. With the nodal prediction from graph networks only, one can still generate continuous solutions with standard interpolation schemes (e.g., using the basis functions in the Galerkin method from FEM). I don’t think it is super challenging to obtain a continuous spatial-temporal solution once we have the finite-dimensional discretization, so I find it hard to see the necessity of fusing OSC into a learning model in the first place.\n\nI am also not sure I agree with a few statements that this paper makes about solving linear systems of equations:\n- “...with the most widely used linear equation solver (i.e., use the inverse matrix to solve the equation)”. Quite oppositely, directly inverting a matrix is perhaps the last thing people should do to solve a large linear system of equations.\n- The time complexity analysis in Sec. 3.5 is also questionable to me: The ABD matrix seems quite sparse, so the time complexity for solving dense linear systems of equations does really apply here.\n- I am also not sure COLROW is the state-of-the-art solver here since it was published nearly 40 years ago. Would it make more sense to consider more recent methods and packages, e.g., SuiteSparse, Eigen, or AMGCL? On a related note, if the ABD matrix is SPD, more dedicated solvers (e.g., Pardiso or CHOLMOD) can apply.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "ICyCtOSbcm", "reviewer": "Reviewer_mydW", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The aurhors provide GRAPHSPLINENETS, a novel deep learning approach to speed up simulation of physical systems with spatio-temporal continuous outputs by exploiting the synergy between graph neural networks (GNN) and orthogonal spline collocation (OSC).\n\nTwo differentiable OSC (time-oriented OSC and spatial-oriented OSC) are applied to bridge the gap between discrete GNN outputs and generate continuous solutions at any location in space and time without explicit prior knowledge of underlying differential equations.", "summary_of_the_review": "I have some concerns on \"experiment design\", \"motivations\", \"evaluations\"", "strength_and_weaknesses": "Strength\nThe physics application is quite interesting\nCombining deel learning with physics are interesting\n\n\n\n\nWeakness\n\n\n1 experiment designs: \nHeat Equation:  Damped Wave:  Navier-Stokes are simulated? not real world data? \nsimulated data will be much easier to fit. \n\n2 motivation and clarity\nnot fully board with the motivation of Time–oriented OSC/Space–oriented OSC:  Is Time–oriented OSC/Space–oriented OSC for introducing continuous solution on mesh grid data?\n\nI am not fully on board on why using GNN is helpful, because mesh-grid data?\n\n3 evaluation\nTable 1: existing methods seem not for this tasks specific (i.e. GNN in general), so performance on this tasks may not be good. \ncompare with methods that  without or without  “space and time continuous simulations” are more relevant. GEN is one of the example. \n Table is shows larger gain on 5s versus 1s, are there some discussion on this? \n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject