Learning Efficient Surrogate Dynamic Models with Graph Spline Networks
Authors: Chuanbo Hua
Federico Berto
Michael Poli
Stefano Massaroli
Jinkyoo Park
Venue: NeurIPS 2023 poster (2023)
Keywords: Graph
Spline Collocation Method
Graph Neural Networks
Simulation
Partial Differential Equations
PDEs
Physics
Scientific Computing
Surrogate Models
Weather Forecasting
TL;DR: We propose a novel model to exploit the synergy between graph neural networks and orthogonal spline collocation to accelerate learned simulations of physical systems by interpolating solutions of graph neural networks.

Abstract: While complex simulations of physical systems have been widely used 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 method to speed up the forecasting of physical systems by reducing the grid size and number of iteration steps of deep surrogate models. Our method uses two differentiable orthogonal spline collocation methods to efficiently predict response at any location in time and space. Additionally, we introduce an adaptive collocation strategy in space to prioritize sampling from the most important regions. GraphSplineNets improve the accuracy-speedup tradeoff in forecasting various dynamical systems with increasing complexity, including the heat equation, damped wave propagation, Navier-Stokes equations, and real-world ocean currents in both regular and irregular domains.

Learning Efficient Surrogate Dynamic Models with Graph Spline Networks
Chuanbo Hua" KAIST DIf AI4CO
Federico Berto" KAIST DIffeqML, AI4CO
Michael Poli
Stefano Massaroli
Jinkyoo Park KAIST OMELET
Stanford University DiffeqML
Mila and University of Montreal Di
Abstract
While complex simulations of physical systems have been widely used 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 GRAPHSPLINENETT up the forecasting of physical systems by reducing the grid size and number iteration steps of deep surrogate models. Our method uses two differentiable or- thogonal spline collocation methods to efficiently predict response at any loca- tion in time and space. Additionally, introduce adaptive collocation strat- egy in space to prioritize sampling from the most important regions. GRAPH- SPLINENETS improve the accuracy-speedup tradeoff in forecasting various dy- namical systems with increasing complexity, including the heat equation, damped wave propagation, Navier-Stoker equations, and real-world ocean currents in both regular and irregular domains.
novel deep-learning method to speed
Introduction
For growing variety of fields, simulations of partial differential equations (PDEs) rep- resenting physical processes are an essential tool. PDE-based simulators have been widely employed in range of practical applications, spanning from astrophysics (Mücke et al., 2000) to biology (Quarteroni and Veneziani, 2003), engineering (Wu and Porté- Agel, 2011). finance, (Marriott al., 2015) or weather forecasting (Bauer et al., 2015). However, traditional solvers for physics- basedd simula- tion often need significant amount of computational resources (Houska et al., 2012), such solvers based on first principles and the modified Gauss-Newton meth- ods. Traditional physics- ahysics-based simulation heavily relies on knowledge of underlying physics and parameters re- quiring ad-hoc modeling, which is sensitive to design choices. Finally, even the best traditional simulators are often inaccurate due to the difficulty in approximating real dynamics (Kremer and Hancock, 2006; Oberkampf, 2019; Sanzhez-Gonzalez al., 2020).
Figure 1.1: Wind speed forecasts on the Black Sca. GRAPHSPLINENETS can leam fast and accurate surrogate models for com- plex dynamics on irregular domains.
attractive alternative traditional simulators use deep learning to train surrogate models directly from ob- served data. Data-driven surrogate models can generate
"Equal contribution authors.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

predictions based on the knowledge learned automatically from training without the knowledge of the underlying differential equations. Among deep learning methods, graph neural networks (GNNss have desirable properties such as spatial equivariance and translational invariance, allowing learning representations of dynamical interactions in a generalizable manner (Pfaff et al., 2021; Bronstein et al., 2021) and on unstructured grids. Despite the benefits of these paradigms, by increasing res- olution, graph-baset based models require significantly heavier calculations. Graph models as Poli et al. (2019); Lienen and Günnemann (2022) predict in continuous time or continuous space. However, the substantial computational overhead due to the need for iterative evaluations a vector field over time requires considerable computation and limits their scalability (Xhonneux et al., 2020).
In this work, we propose GRAPHSPLINENETS a novel approach that exploits the synergy between GNNs and the orthogonal spline collocation (OSC) method (Bialecki and Fairweather, 2001; Fair- weather and Meade, 2020) to both efficiently and accurately forecast continuous responses of a target system. We use the GNNs to efficiently get the future states on coarse spatial and temporal grids and apply OSC with these coarse predictions to predict states at any location in space and time (i.e., continuous). As a result, our approach can generate high-resolution predictions faster than previous approaches, even without explicit prior knowledge of the underlying differential equation. Moreover, the end-to-end training and an adaptive sampling strategy of collocation points help the GRAPHSPLINENETTS to improve prediction accuracy in continuous space and time.
We summarize our contributions as follows:
We introduce GRAPHSPLINENETS, a novel learning framework to generate fast and accurate predictions of complex dynamical systems by employing coarse grids to predict time and space continuous responses via the differentiable orthogonal spline collocation.
We propose an adaptive collocation sampling strategy that adjusts the locations of collocation points based fast-changing regions, enabling dynamic improvements in forecasting accuracy.
We show GRAPHSPLINENETS are competitive against existing methods in improving both ac- curacy and speed in predicting continuous complex dynamics on both simulated and real data.
2 Related Works
Modeling dynamical systems with deep learning Deep neural networks have recently been suc- cessfully employed in accelerating dynamics forecasting, demonstrating their capabilities in pre- dicting complex dynamics often orders of magnitude faster than traditional numerical solvers (Long al., 2018; Li et al., 2021; Berto et al., 2022; Pathak et al., 2022; Park et al., 2022; Poli et al., 2022; Lin et al., 2023a,b; Boussif et al., 2022, 2023). However, these approaches are typically lim- ited to structured grids and lack the ability to handle irregular grids or varying connectivity. To address this limitation, researchers have turned to graph neural networks (GNNs) as promising alternative: GNNs enable learning directly on irregular grids and varying connectivity, making them well-suited for predicting system responses with complex geometric structures or interactions (Li al., 2022a). Additionally, GNNs inherit physical properties derived from geometric deep learning, such as permutation and spatial equivariance (Bronstein et al., 2021), providing further advantages for modeling complex dynamics. Alet et al. (2019), one of the first approaches to model dynamics with GNNs, represent adaptively sampled points in graph architecture to forecast continuous un- derlying physical processes without any a priori graph structure.
Gonzalez al. (2020) extend the GNN-based models to particle-based graph surrogate models with dynamically changing connectivity, modeling interactions through message passing. Graph neural networks have also re- cently been applied to large-scale weather predictions (Keisler, 2022). However, GNNs are limited in scaling in both temporal and spatial resolutions with finer grids as the number of nodes grows.
Accelerating graph-based surrogate models An active research trend focuses on improving the scalability of graph-based surrogate models. Pfaff et al. (2021) extend the particle-base deep mod- els of (onnalz al., 2020) to mesh-based ones. By employing the existing connectivity of mesh edges, Pfaff et al. (2021) demonstrate better scalability than proximity- based online graph creation modeling particle interactions. Li et al. (2022b) extend convolutional operators to irregular grids by learning a mapping from an irregular input domain to a regular grid: under some assump- tions on the domain, this allows for faster inference time compared to message passing. Fortunato et al. (2022) build on the mesh-based approach in (Sanchez-Goonaler, et al., 2020) to create multi-

---

scale GNNs that, by operating on both coarse and fine grids, demonstrate better scalability com- pared to single-scale GNN models. GRAPHSPLINENETS do not need to learn additional mappings between different domains or need additional parameters in terms of different scales thus reducing computational overheads by operating directly on unstructured coarse meshes. Our method bridges the gaps between coarse and fine grids in time and space via the orthogonal spline collocation (OSC), thus improving GNN-based methods' scalability.
Collocation methods Collocation and interpolation methods are used to estimate unknown data values from known ones (Bourke, 1999). A fast and accurate collocation method with desir- able properties in modeling dynamics is the orthogonal spline collocation (OSC) (Bialecki, 1998; Bourke, 1999). Compared to other accurate interpolators with O time complexity. the OSC has complexity of only On log n.) thanks to its sparse structure, allowing for fast inference; more- over, is extremely accurate, respects C2 continuity, and has theoretical guarantees on convergence (Bialecki, 1998). Due to such benefits, few research works have aimed at combining OSC with deep learning approaches. Guo et al. (2019); Brink et al. (2021) combine deep learning and collo- cation methods by directly learning collocation weights; this provides advantages over traditional methods since the OSC can ensure C²l continuity compared to mesh-based methods. A recent work Boussif et al. (2022) has studied learnable interpolators and demonstrated that it is possible to learn an implicit neural representation of a collocation function. GRAPHSPLINENETS differ from previ- neural collocation methods in two ways: first, decrease the computational overhead by not needing to learn interpolators or collocation weights; secondly, by interpolating not only in space but also in time, our method enables considerable speedups by employing coarser temporal grids.
3 Background
3.1 Problem Setup and Notation
We first introduce the necessary notation that we will use throughout the paper. Let superscripts denote time, and subscripts denote space indexes as well as other notations. We represent the state of physical process space location X = {xt, i = I,... ,N} and time ( R as Y {uf 1.... , N) where N represents the number of sample points and Ω c TD is D-dimension physical domain. A physical process in this domain can be described by a solution of PDEs, i.e., of = u(x,,t), t). The objective of dynamics surrogate model is to estimate future states (Y', t E R+} given the initial states Y
3.2 Message Passing Neural Networks
We employ graph G [vf, efy) to encode the status of physical process at sample points X at t. v and e/, denote the attribute of sample node and the attribute of the directed edge between node i and j. respectively. Node attributes are encoded from sample points state information while edge attributes are the distance between every two nodes. Message passing neural networks (MPNNs) employ an corr--pooosssodecodeeeee structure to predict the states of sample points at the next timestep by giving the previous states:
D
Pm (E(---(P(E(X,)))) Y.))
(1)
decodar
processor
ancodar
where E(- is the encoder, P.(.) is the i-th message passing layer, and (--) is the decoder.
3.3 Orthogonal Spline Collocation Method
The OSC method consists of three steps in total: (1) partitioning the input domain and selecting collocation points, (2) generating polynomials with parameters to be determined with observations, and (3) solving equations to determine the parameters.
*Collocation and interpolation are terms that are frequently used interchangeabll. While interpolation detined 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 as a flexible subset of interpolation methods that satisfies certain conditions, such as C² class continuity.
3

---

Ground Truct
Polysomill 1 Simrarion
Polysomal 2 Simuration
Polysomil 3 Simulation
0.0
Partiion Points
Collocation
0.2
Figure 3.1: Visualization of an 1-D OSC example with the order and partition number N 3.
Partitioning and initializing polynomials The OSC aims find series of polynomials under order r while satisfying C² continuity to approximate the solution of PDEs. To find these polyno- mials, we split each dimension of the domain into N partitions* Then, we initialize one r-order polynomial with r ++1 1 unknown coefficients in each partition. These polynomials have N x (r (r+1) degrees of freedom, i.e., the variables be specified to determine these polynomials uniquely.
Collocation points selection and equation generation We firstly have 2 equations from the boundary condition and (N 1 1) equations from the C continuity restriction. Then need N x (r-1) 1) more equations to uniquely define the polynomials. So we select r collocation points where the observation obtained to determine the parameters each partition. By substituting the state of collocation points to the polynomials, we can get N x (r 1) equations. Now we transfer the OSC problem to an algehraic problem. With properly selected base functions, this algebraic equation will be full rank and then has unique solution (Lee and Micchelli, 2013).
Solving the equations The coefficient matrix of the generated algebraic problem is almost block diagonal (ABD) (De Boor and De Boor, 1978). This kind of system allows for efficient computa- tional routines (Amodio al., 2000), that we introduce in $ 3.4. By solving the equation, we obtain the parameters for the polynomials that can be used to predict the value of any point in the domain. We illustrate an example of 1-D OSC problem in Fig. 3.1. More details and examples of the OSC are provided in Appendix A. The OSC method can also be extended to 2-D and nigher-dimensional do- mains (Bialecki and Fairweather, 2001). In our work, we employ the 1-D and 2-D OSC approaches in the time and space domains, respectively.
3.4 Efficiently solving ABD matrices
Most interpolation methods need to solve linear equations. Gaussian elimination is one of the most widely used methods to solve a dense linear equation, which, however, has a O(n²) complexity (Strassen et al., 1969). Even with the best algorithms known to date, the lower bound of time complexity to solve such equations is On?logg) (Golub and Van Loan, 2013).
In the OSC method, the coefficient matrix for the lin- ear equation follows the ABD structure, which we efficiently solve with time complexity ((n?)
(b)
103 Iteration Running Time
[Gaanke))
Calculation Time
by the COLROW algorithm (Diaz et al., 1983) as
shown in Fig. 3.2 (a). The core idea for this method
10-
is that by using the pivotal strategy and elimina- tion multipliers, we can decompose the coefficient matrix into a set of a permutation matrix and up- per or lower triangular matrix that can be solved O(n²) time each. The latest package providing this algorithm in the FORTRAN programming lan- guage: our re-implementation in Py PyCCh (Paszke al., 2019) allows for optimized calculation, GPU support, and enabling the use of automatic differen- tiation, as shown in Fig. 3.2 (b).
Digh
Number of Elements
Figure 3.2: (a) The COLROW algorithm has
lowest computational complexity with O(n") compared to the Gaussian elimination algorithm with O(re2) or lower bounds of the generic sparse solvers with O Orrrlgg)). (b) Running time for 100 iterations of OSC+LU decomposition, B: OSC+COLROW solver, and OSC+COLROW solver with GPU accelcration.
P Note that these partitions do not need to be isometric, thus allowing flexibility in their choice. "We employ the Gauss- LegendLe quadrature for choosing the collocation points in regular boundary prob- lems. We further discuss the choice of collocation points in Appendix B.I.

---

Methodology
@ Space-oriented OSC For each time frame, we take l example, the neural network gen- erates prediction values of sample points which can be considered 2-D OSC problem. In this space-oriented OSC problem, collocation points are positions of sample points {x;). their values at this time frame [ ), and the output is polynomials / (x), E Ω which provides continuous prediction the domain. Now we get the prediction of every time frame at any position. Note that the space-oriented OSC can use the time-oriented OSC as the input, which means that we can use the space-oriented OSC at any time frame and then get a spatio-temporal continuous predictions a(x, t), (x, Ω x [to, tx]. 4.1 GRAPHSPLINENETS
Fig. 4.1 depicts the entire architecture of our model. GRAPHISPLINENETS can be divided into three main components: O message passing neural networks (MPNN), time-oriented collocation and space-oriented collocation. The MPNN takes observations at collocation points as input and infers a sequence of discrete future predictions via autoregressive rollouts. We then use the time-oriented and space-oriented collocation methods on these discrete predictions to obtain functions that can provide both time and space-continuous predictions.
NN
NN
TL
fx
Biy At
1/1
Space oriented OSC
Tome-oriented
Eqution ABO mants
lines Equaiin
Stepktion Equatic g1 (x44)
Sriveinn Equatir Rilfe + Af)
Thodeer
Bloonder
(X, Y
cor G
G
P
Y
Constracs Graph
Puo Montge
Get Upskaed Stao
Dly..
olja - Ploea)
Figure 4.1: O Message passing neural networks take as inputs current states and employ an encoder-pricessor- decoder architecture for ohtaining the next states autoregressively. time-oriented OSC obtains continuous predictions in time while Space-oriented OSC is employed obtain continuous predictions in space. The model trained end-to-end efficiently by leveraging sparsity almost block diagonal (ABD) matrices of the collocation method.
Assign collocation points We follow the same rule to select collocation points in 2-D domain described in $ 3.3 to select sample points (x,) for 2-D domain. We have the states of these sample points at the initialization time frame {u/"}.
O Message passing neural networks The neural networks take the states of sample points at the initial time frame ((x;, y²")) }} input and employ prcoder-proessso-dccodee architecture for obtaining the next steps states 0."² 1, autoregressively.
Time-oriented OSC For each sample point, taking x; as an example, the neural networks gen- erate sequence of predictions along time which be considered time-oriented OSC problem, collocation points are time steps (to, t1, ... .tk) and values at these time steps {u."" 5/ .... of J. and the output is the polynomials à; (t), [to, tx] which provides continuous prediction along time; thus, can obtain predictions at any time.
1-D OSC problem. In this
5

---

4.2 Training strategy and loss function
GRAPIISPLINENETS predicts the fine resolution states in hierarchical manner. MPNN uses the sample points at the initialized state as input to propagate rollouts on these points as shown in Fig. 4.2. Based on MPNN outputs, two OSC problems are solved to find the polynomial parameters = OSC(MPNN(X, ()) which allow for spatiotemporal continuous outputs ú(x,, t&: Q). Note that solving OSC problems is equivalent to predicting the future space and time continuous states, given the sample inputs.
Inherently, training the parameters, A of the MPNN is cast as bi-level optimization problem as 0 argming where e OSC(MPNN(X, Y 0)). The loss function L is defined as:
I. sample points arror
10tt
0/113
Σ
lu th tx)
(2)
L, = Interpolation points error
The left term L, of Eq. (2) is the loss of the MPNN output y where NN. is the number of spatial sample points, and K is the number of MPNN rollout steps. The left term L, solely depends on the MPNN
parameters 0. The right term I2 is the pre-
to
tk
diction loss of the high-resolution interpola- N.,N, tion points (x,, vs === along time and space, where N; is the number of spatial interpolation points at the one-time frame, and N, is the number of interpolation time frames. Note that the space and time con- tinuous prediction i(x,, 44i Φ) is constructed based upon the optimized parameters & = OSC(MPNN(X, Y 8)), the parameters that are obtained by solving OSC problems with MPNN predictions, MPNN(X, (). Thus, both loss terms can be optimized with respect to the MPNN parameters 9, and the whole model is trained end-to-end with automatic differentiation through the OSC. This optimization scheme is akin to Bialas and Karwan (1984); Tuy et al. (1993) that solve bi-level optimization by flattening the loss for the upper- and lower-level problems.
Somple Points
Interpolation Points
Illustration of data points used for model training. Interpolation points include ? time-oriented OSC and space-oriented OSC interpolation points. Here, i can be any time frame in (lo, (xo.
Figure 4.2:
4.3 Adaptive collocation sampling
To allow for prioritized sampling of important location regions, we dynamically optimize the po- sitions of collocation points via their space gradients of the partial derivative values over time. More specifically, after getting the continuous sim- ulation function ú, the partial derivative equation of time provides the changes of the domain with time. Then, we calculate the maximum gradient vector for each collocation point to move them to more dynamic regions, i.e., with the largest change over time, which causes higher errors in the static collo- cation points condition. We formulate this process by
Static
Adaptive
Figure 4.3: Adaptive collocation strategy: mesh points converge towards areas with faster chang- ing regions. this example, collocation points move towards the moving wave crests.
du (x, (x)
X = x; + Bv, argmax, VV
(3)
where 3 is the coefficient of the maximum gradient vector. We project collocation points back to the partition if the gradient step moves them outside to ensure a sufficient number collocation points in each partition cell. We illustrate the adaptive collocation points in Fig. 4.3. We use the states at optimized positions adapted from history rollouts as the subsequent rollouts input. By adjusting the collocation points position, our model can focus on more dynamic parts of the space to get more accurate prediction.

---

Experiments
5.1 Datasets and Training Settings
five dynamical systems of increasing challenge. Simulated datasets consist of three PDEs (partial differential equations), including Heat, Damped Wave, and Navier-Stokes Equations. The two empirical datasets include the Ocean and Black Sea datasets; the Black Sea introduces significant challenges in scalability in the number of nodes, complex dynamics, irregular boundary, and non-uniform meshing. All the models employ the same structure of encoder- processor-decodees for fair comparisons and the same amount of training data in each testing domain. While inputs of baseline models are directly all of the available data points, inputs of our OSC-based models are only an initialized 12 x 12 collocation point unitary mesh at the initial state and fewer in- time sample points. For the Black Sea dataset, we construct a mesh via Delaunay Triangulation with 5000 data points for baselines and 1000 for GRAPHSPLINEN (2022). More details on datasets available in Appendix B. We value open reproducibility and make our code publicly available".
We evaluate GRAPHSPLINENT
as done in Lienen and Günnemann
5.2 Evaluation metries and baselines
We evaluated our model by calculating the mean square error (MSE) rollout prediction steps with respectively ground truth for f E [0,5] 5] seconds. We employ relevant baselines in the field of discrete-step graph models for dynamical system predictions. Graph convolution networks (GCNs) (Kipf and Welling, 2016) and GCN with hybrid multilayer perceptron (MLP) model are 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 predictions in space, iee., GEN (Alet et al., 2019). A similar setup considering the inherent graph structure and employing MPNNs in mesh space is used by (Pfaff et al., 2021) (MeshGraphNee,, MGN in our comparisons). We utilize MGN as the MPNN building block for our GRAPHSPLINENETS.
5.3 Quantitative Analysis
We consider series of ablation models on four datasets to demonstrate the effectiveness of our model components in our approach in multiple aspects. Quantitative results of ablation study models are shown Table 5.1 [right]. The ablated models are:
MGN: MeshGraphNet model with 3 message passing layers from Pfaff et al. (2021).
MGN+OSC(Post): model with 3 message passing layers and only pest-processing with the OSC method, i.e., we firstly train a MGN model, then we use the OSC method collocate the prediction as a final result without end-to-end training.
MGN+OSC: MGN with OSC-in-the-loop that allows for end-to-end training.
MGN+OSC+Adaptive: MGN+OSC model that additionally employs our adaptive colloca- tion point sampling strategy.
Post processing end-to-end learning We show the effectiveness of end-to-end learning ar- chitecture by comparing the MGN+OSC(Post) and MGN+OSC models. Table 5.1 shows that MGN+OSC has more accurate prediction than MGN+OSC(Post) by more than 8% percent across datasets. This can explained by fact that, since the OSC is applied end-to-end, the error between MGN prediction steps is backpropagated to the message passing layers. In contrast, the model has no way of considering such errors in the post proc-psressiing steps.
Adaptive collocation points We further show the effectiveness of the adaptive colloca- tion strategy by comparing the MGN+OSC and MGN+OSC+Adaptive Table 5.1 shows that MGN+OSC+Adaptive has a more accurate prediction than MGN+OSC, i... more than 5% im- provement on long rollouts in the Black Sea dataset. Adaptive collocation points encourage those points to move the most dynamic regions the domain, which is not only able to place
7

---

Table 5.1: Mean Squared Error propagation and runtime for 5 second rollouts, GRAPHSPLINENETS consis-
tently outperform baselines both accuracy and runtime. Smaller better (.). Best in bold; second underlined.
Dataset
Metric
Haselines
GRAPESPLINENETS
GCN
GCN+MLP
GEN
MGN
MGN++SC((PsS)
MGN+OSC
MGN-OSC-AÁpiix
Heat Equation
6.87 4: 1:00
5.02 4: 0.59
2.52 à: 0.23
3.01 4: 0.38
1,8 0.18
1
1.07 0.258
Rumtime [s]
3.26 1 0.12
3,02 1 0.100
6.87 1. 0.10
6.99 0.12
1.52 0.09
1.38 0.30
1.41 1 0.122
Damped Wave
MSE (x10-1)
10.5 1.65
9.90 1.52
0.62
7.80 0.88
4,96
4.60
4.51
Runtime [s]
0.95 1: 0.08
0.82 1: 0.07
1.13 1: 0.08
1.38 a: 0.10
0.45 0.055
0.39 0.04
0.42 4: 0.09
Navier Stokes
MSE (x1001)
4.24 0.95
3.91 1: 0.00
3.45 1. 0.24
3.66 1 0.33
2.58 0.28
2,21 1 0.27
2.02 1 0.30
Runtime
0.91 d: 0.08
0.88 4 0.07
1.01 0.09
1.21 4: 0.10
0.51 0.05
0.47 0.04
0.49 4 0.09
Ocean Currents
(x16") 1>
0.02 1: 1.12
5.20 a: 0.98
1.55 a: 0.52
4.75 a: 0.61
3.82 1: 0..09
3.01 1: 0.38
3.34 0.4.
3.56 0.15
3.38 + 0:13
7:11 + 0.31
7.15 0.25
1.61 0.12
1.44 0.09
0.11
Black Sea
MSE (xx001) 1
1.33 4 0.28
1.12 4: 0.35
6.22 0.31
6.41 +: 0.42
4.37 0.33
4:23 d: 0.17
3.91 0.27
Runtime [s]
8.27 1 0.33
8.11 1 0.28
12.81 1 0.19
13.35 0.27
5.87 1. 0.13
5.13 0:21
5.79 1 0.15
greater attention hard-to-learn parts in space but also can let the OSC method implicitly de- velop a better representation of the domain. Empirical quantitative results on the five datasets are shown in Table 5.1. In the heat equation dataset, our approach reduces long-range prediction er- rors by 64% with only 20% of the running time compared with the best baseline model. In the damped wave dataset, our approach reduces errors by 42% with 48% reduction in inference time.
In the Navier-Stokes dataset, our method reduces 31% long-range prediction errors while requiring 37% less time to infer solutions compared to the strongest baseline.
MGN-05C
Neural networks inference time
(GRU)=Acaptiow
OSC inference time Adapatie location time
MGN--05C
(GPU)
MGN-05C
Inference times We show the effectiveness of the COLROW algorithm in accelerating the OSC speed by comparing the OSC method with one of the most commonly used algorithms for efficiently solving linear systems and the OSC with the COLROW solver in Fig. 3.2. Our differentiable OSC and adap- tive collocation only result in a minor overhead over the graph neural network. The overall scheme is still a fraction of the total in terms of inference time as shown in Fig. 5.4, while making GRAPHSPLINENE more accurate.
(Pret, GPU)
MGN-0S5 (Poz, CPU)
MGN
Running Time
Figure 5.1: Inference time benchmark: the OSC and adaptive collocation points result minor overhead; the overall scheme result a fraction of the baseline running time.
5.4 Sensitivity Analysis
Interpolation and collocation methods We demon- strate the efficiency of the OSC method by comparing the combination MGN with different interpolation and collocation methods, including linear interpolation, 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 di- mensions. Results are shown in Fig. 5.2 where we mea- sured the mean square error and running time of 3 second rollouts predictions by varying the number of collocation points from (2 x 2) to (16 x 16). The model MGNC shows the best accuracy while having shorter running time. Even though the linear interpolation can be slightly faster than the OSC, it shows considerable error in the prediction and does not satisfy basic assumptions such as Lipschitz continuity in space.
10
10-1
MSE
10-2
10-1
10-1
10-1
If
Prediction Running Time [s]
MGN Linear
MGN Cubic
MGN OSC MGN Spline
MGN OSC Adaptive
GRAPHSPLINENETS (MGN-OSC+Adaptive) is Pareto-optimal combinations with other interpolation methods.
Figure 5.2:
compared
"We employ for our experiments torch. linalg solve, which uses LU decomposition with partial piv- oting and row interchanges. While this method is numerically stable, it still has a O(r²) complexity with a Orre log(n)) lower bound when employing generic sparse solver algorithms.
8

---

Number of collocation points We study the effect of the number of collocation points on the 3 second rollout prediction error by testing the MGN, MGN+OSC, and MGNaoSS+AAaapiive models.
The MGN is directly trained with the whole domain data.
We use a different number of collocation points (from (2 2) to (28 x 28)) in 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 the number of collocation points, Fig. 5.3 shows that the MGN+OSC and MGN+OSC+Adaptive achieve significant improvements in prediction accuracy over the MGN. The MGN+OSC+Adaptive has a stable better per- formance than the MGN+OSC, and the improvement is larger when there are fewer collocation points. The reason is that with fewer collocation points, the MGN+OSC has insufficient ability to learn the whole domain. With the adaptive collo- cation point, the MGN+OSC+Adaptin can focus hard-to- learning regions during training to obtain overall better predic- tions.
10
MGN- + OSC
MGN OSC Adaptive
MSE
[B 88) 116 x4)
Number of Collecation Points
Figure 5.3: MSE vs the number of collocation points. MGN does not use collocation points, hence constant. The more the collocation points, the better the performance.
Number of rollout steps We show the effectiveness of the OSC method in improving long-range prediction accuracy by comparing the MGN and MGN+OSC model. Fig. 5,4 shows the MGN+OSC can keep stable in long-range rollouts compare with the MGN. The reason is that with the OSC, we can use fewer neural net- work rollout steps to obtain longer-range predictions, which avoids error accumulation during the multi-step rollouts and implicitly learns for compensating integration residual errors. In addition, end-to-end learning lets the neural networks in MGN+OSC learn the states between rollout steps, making the prediction stable and accurate.
MSE
Number of Message Passing Layers We further show- case the sensitivity to the number of message-passing lay- ers while adopting the same baseline processor for As Fig. 5.5 shows, a higher num- ber of message-passing layers generally helps in decreas- ing the MSE at the cost of higher runtime, reflecting find- ings from prior work such as GNS (Alet et al., 2019) and MGN (Pfaff et al., 2021). Notably, GRAPHSPLINENETS exhibits a Pareto-optimal trade-off curve, outperforming MGN in terms of prediction accuracy and computational efficiency different numbers of message-passinng layers. Downsampling mesh nodes can help in robustness and generalization given the larger spatial extent of communication akin to Fortunato et al. (2022); the OSC plays crucial part in handling the problem's 's continuity (e.g., in time) while adaptive collocation effectively helps the model to focus on regions with higher errors thus increasing the model performance.
40
69
Rollout Steps
MEN Adgetie
fair comparison.
Figure 5.4: MSE vs rollout steps. Our full model yields stable long rollouts.
Heat Equation
Mavie-Stoien Equation
Scalability to Higher Dimensions We show- case GRAPHSPLINENTE on a larger 3D dataset, namely the 3D Navier-Stokes equations dataset from PDEBench (Takamoto et al., 2022) (more details in Appendix B.7). The number of nodes in Graph- SplineNets is 4096, which is more than 4x compared to the Black Sea Dataset. Moreover, the 3D Navier- Stokes includes total of 5 variables: density, pres- sure, and velocities on 3 axes. For a fair comparison, we keep the model structure the same as in other ex- periments and compare GRAPHSPLINENT
25
Lopecc
Layes,
Ligers,
1ine
Run
10
MSELoss
M5ELos
Figure 5.5: Pareto plot for MSE and runtime at different # of layers.
against MGN. MGN obtain an MSE of (4.52 0.23) x 10-1 with runtime of more than 80 seconds, while CRRAHNSPLNNNN remarkably outperforms MGN with an MSE of $3.9++0.277 x 10-1 and less than 9 seconds to generate a trajectory an example of which is shown in Fig. 5.6.

---

8
GooundTan Truh
Groand Trsh
Ground Truth
Groand Tush
tet Contee
Sanu
Sindarion
Figure 5.6: Rollout visualization of density for 3D Navier-Stokes equations.
5.5 Qualitative analysis
We visualize the prediction of our model with baselines on the Black Sea dataset Fig. 5.7. Our model can accurately predict even with such complex dynamics from real-world data that result in turbulence phenomena while considerably cutting down the computational costs shown in Table 5.1. We also provide series of qualitative visualizations with different datasets and baselines in Appendix C. Our model has smoother error distribution and more stable long-range prediction. Thanks to the continuous predictions from GRAPHSPLINENETS we can simulate high resolutions without needing additional expensive model inference routines, while the other models only achieve rewer-resolution predictions. Baselines visibly accumulate errors for long-range forecasts; our model lower the error with smoother and accurate predictions space and time.
Grand Trath
GUN
GON-MLP
GRN
MGN
OIN MGN44060
MIGN+08CY bddgen
SARE
Figure 5.7: Visualization rollout forecasts for baselines and TRAAPHSPLINNNET on the Black Sea dataset. First row: predictions of the wind speed norm. Second row: mean square error (MSE). Our full model, namely MGN+OSC(Adaptive) generates more accurate predictions than other baseline surrogate models.
Conclusion and Limitations
In this work, we proposed GRAPHSPLINENETS, novel approach that employs graph neural net- works alongside the orthogonal spline collocation method to enable fast and accurate forecasting of continuous physical processes. We used GNNs obtain predictions based on coarse spatial and temporal grids, and the OSC method to produce predictions at any location in space and time. Our approach demonstrated the ability to generate nigh-resolution predictions faster than previous methods, even without explicit prior knowledge of the underlying differential equations. Addition- ally, we proposed an adaptive collocation sampling strategy that improves the prediction accuracy as the system evolves. We demonstrate how GRAPHSPL. INENETS are robust in predicting complex physics in both simulated and real-world data. We believe this work represents a step forward in new direction in the research area at the intersection of deep learning and dynamical systems that aims at finding fast and accurate learned surrogate models.
A limitation of our approach is that, by utilizing the orthogonal spline collocation (OSC) method for interpolation, we assume some degree of smoothness in the underlying physical processes between collocation points; as such, our approach may struggle when modeling discontinuity points, which would be smoothed-out. Finally, we note that despite benefits, deep learning models including PRRAPHSPLINENN
often lack interpretability. Understanding how the model arrives at its predic- tions or explaining the underlying physics in a transparent manner can challenging, hindering the model's adoption in domains where interpretability is crucial.
10

---

Acknowledgements
We want to express our gratitude towards the anonymous reviewers who greatly helped us improve our paper. This work was supported by the Institute of Information & communications Technology Planning & Evaluation (ITTP) grant funded by the Korean govemneen(MSSTTyO22-000000 velopment of Collective Collaboration Intelligence Framework for Internet of Autonomous Things),
De-
References
F. Alet, A. K. Jeewajee, M. B. Villalonga, A. Rodriguez, T. Lozano-Perez, and L. Kaelbling, Graph element networks: adaptive, structured computation and memory. In International Conference on Machine Learning, pages 212-222. PMLR, 2019.
P. Amodio, J. Cash, G. Roussos, R. Wright, G. Fairweather, I. Gladwell, G. Kraut, and M. Pa- przycki. Almost block diagonal linear systems: sequential and parallel solution techniques, and applications. Numerical linear algebra with applications, 7(5):275-317, 2000.
P. Bauer, A. Thorpe, and G. Brunet. The quiet revolution of numerical weather prediction. Nature, 525 2015.
Berto, S. Massaroli, M. Poli, and J. Park. Neural solvers for fast and accurate numerical optimal control. International Conferences on Learning Representations, 2022.
Berto, C. Hua, J. Park, M. Kim, H. Kim, J. Son, H. Kim, J. Kim, and J. Park. R14co: an extensive reinforcement learning for combinatorial optimization benchmark. arXiv preprint SEIi:2306611100 2023.
W. Bialas and M. Karwan, Two-level linear programming. Management Science, 30(8), 1984.
B. Bialecki. Convergence analysis of orthogonal spline collocation for elliptic boundary value prob- lems. SIAM journal on numerical analysis, 35(2):617-633 1998.
B. Bialecki and G. Fairweather. Orthogonal spline collocation methods for partial differen- tial equations. Journal of Computational and Applied Mathematics, 128(1):55-82, 2001. http://o....01000000116777700010
ISSN 0377-0427. doi: www. www.scieete scienceiieect..... wtt..cco/scceee st//sciencee artt:clle/p31//8770000000000 Numerical Anal- ysis 2000. Vol. VII: Partial Differential Equations.
URL https: //
Bourke. Interpolation methods. Miscellaneous: projection, modelling, rendering, 1(10), 1999.
O. Boussif, D. Assouline, L. Benabbou, and Y. Bengio. Magnet: Mesh agnostic neural pde solver. Advances in Neural Information Processing Systems, 2022.
O. Boussif, G. Boukachab, D. Assouline, S. Massaroli, T. Yuan, L.. Benabbou, and Y. Bengio. What if we enrich day-ahead solar irradiance time series forecasting with spatio-temporal context? Ad- vances in Neural Information Processing Systems, 2023.
A. R. Brink, D. A. Najera-Flores, and C. Martinez. The neural network collocation method for solving partial differential equations. Neural Computing and Applications, 33(11):5591-5600 2021.
M. M. Bronstein, J. Bruna, T. Cohen, and P. Veličković. Geometric deep learning: Grids, groups, graphs, geodesics, and gauges. arXiv preprint arXiv:2104. 13478, 2021.
S. A. Ciliberti, M. Grégoire, J. Staneva, A. Palazov, G. Coppini, R. Lecci, E. Peneva, M. Matreata, V. Marinova, S. Masina, et al. Monitoring and forecasting the ocean state and biogeochemical processes in the black sea: recent developments in the copernieus marine service. Journal of Marine Science and Engineering, 9(10): 146, 2021.
C. De Boor and C. De Boor. practical guide to splines, volume 27. springer- New York, 1978.
J. Diaz, G. Fairweather, and P. Keast. Fortran packages for solving certain almost block diagonal linear systems by modified alternate row and column climination. ACM Transactions on Mathe- matical Software (TOMS). 93))358--373 1983.
G. Fairweather and D. Meade. A survey of spline collocation methods for the numerical solution of differential equations. In Mathematics for large scale computing, pages 297-341. CRC Press, 2020.
11

---

W. Falcon et al. Pytorch lightning, GitHub. Note: https://github. comPHyyDorilLighnnnne/tooooooo lightingg 3(6), 2019.
M. Fortunato, T. Pfaff, P. Wirnsberger, A. Pritzel, and P. Battaglia. Multiscale meshgraphnets. arXiv preprint urXiv:2210.0000 2022.
G. H. Golub and C. Van Loan. Mutrix computations. JHU press, 2013.
H. Guo, X. Zhuang, and T. Rabezuk. A deep collocation method for the bending analysis of kirchhoff plate. Computers, Materials & Continua, 59(2)4333456, 2019. ISSN 1546-2226. doi: 10.32604/
amc.2019.06660 URL http://dx. doi. org/18 10.
32604/ cmc. 2019. 06660.
B. Houska, F. Logist, M. Diehl, and J. Van Impe. A tutorial numerical methods for state and parameter estimation in nonlinear dynamic systems. Identification for Automotive Systems, pages
67-88, 2012.
R. Keisler. Forecasting global weather with graph neural networks arXiv preprint arXiv:2202.07575, 2022.
D. P. Kingma and J. Ba. Adam: A method for stochastic optimization.
arXiv preprint
arXiv:14/2.6991 14/2.690, 2014.
T. N. Kipf and M. Welling. Semi-supervised classification with graph convolutional networks. arXiv preprint ar arXiv:1609.000000
2016.
M. Kirchmeyer, Y. Yin, J. Donå, N. Baskiotis, A. Rakotomamonjy, and P. Gallinari, Generalizing to new physical systems via context-informed dynamics model. arXiv preprint arXiv:2202.01889, 2022.
D. Kremer and B. Hancock. Process simulation in the pharmaceutical industry: a review of some basic physical models. Journal of phurmaceutical sciences, 95(3):517-529, 2006.
Y. J. Lee and C. A. Micchelli. On collocation matrices for interpolation and approximation. Journal of Approximation Theory, 174: 148-181, 2013.
J. Li, C. Hua, Park, H. Ma, V. Dax, and M. J. Kochenderfer. Evolvehypergraph: Group-aware dynamic relational reasoning for trajectory prediction. arXiv preprint arXiv:2208, 05470, 2022a.
Z. Li, N. Kovachki, K. Azizzadenesheli, B. Liu, K. Bhattacharya, A. Stuart, and A. Anandkumar. Fourier neural operator for parametric partial differential equations. International Conference Learning Representations, 2021.
Z. Li, D. Z. Huang, Liu, and A. Anandkumar. Fourier neural operator with learned deformations for pdes general geometries. arXiv preprint 207.v=2207.05209, 2022b.
M. Lienen and S. Günnemann. Learning the dynamics of physical systems from sparse observations with finite element networks. International Conference on Learning Representation, 2022.
S. Lin, J. Kim, C. Hua, S. Kang, and M.-H. Park. Comparing artificial and deep neural network models for prediction of coagulant amount and settled water turbidity: Lessons learned from big data in water treatment operations. Journal of Water Process Engineering, 54:103949, 2023a.
S. Lin, J. Kim, C. Hua, M.-H. Park, and S. Kang. Coagulant dosage determination using deep learning-basse based graph attention multivariate time series forecasting model. Water Research, 232: 119665, 2023b.
A. Logg, K.-A. Mardal, and G. Wells. Automated solution of differential equations by the firite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012.
Z. Long, Lu, X. Ma, and B. Dong. Pde-net: Learning pdes from data. In International Conference on Machine Learning, pages 3208-3216. PMLR, 2018
Marriott S. M. Tan, and N. Marriott. Experiential learning-a case study of the use of computerised stock market trading simulation in finance education. Accounting Education, 24(6):480-497, 2015.
Experiientialleammngg--a
S. Marullo, R. Santoleri, D. Ciani, P. Le Borgne, S. Péré, N. Pinardi, M. Tonani, and G. Nardone. Combining model and geostationary satellite data to reconstruct hourly sst field over the mediter- ranean sea. Remote sensing of environment, 146:1 1-23, 2014.
A. Mücke, R. Engel, J. P. Rachen, R. J. Protheroe, and T. Stanev. Monte carlo simulations of photohadronie processes in astrophysics. Computer Physics Communications, 124(2-3)290-314, 2000.
12

---

B. B. Nardelli, C. Tronconi, A. Pisano, and R. Santoleri. High and ultra-high resolution processing of satellite sea surface temperature data over southern european seas in the framework of myocean project. Remole Sensing of Enviromment, 129:1-16, 2013.
W. L. Oberkampf. Simulation accuracy, uncertainty, and predictive capability: A physical sciences perspective. Computer Simulation Validation: Fundamental Concepts, Methodological Frame- works, and Philosophical Perspectives, pages 19-97, 2019.
J. Park, Berto. A. Jamgochian, M. J. Kochenderfer, and Park. Meta-sysid: A meta-lcarning approach for simultaneous identification and prediction. arXiv preprint arXiv:2206. vrXi:22206.00694, 2022.
A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, 7. Lin, N. Gimelshein, I. Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems, 32, 2019.
J. Pathak, S. Subramanian, P. Harrington, S. Raja, A. Chattopadhyay, M. Mardani, T. Kurth, D. Hall, Z. Li. K. Azizzadenesheli, et al. Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators. arXiv preprint arXiv:2202.11211 2022.
T. Pfaff, M. Fortunato, A. Sanchez-Gionzalez, and P. W. Battaglia. Learning mesh-based simulation with graph networks. International Conference on Learning Representations, 2021.
M. Poli, S. Massaroli, J. Park, A. Yamashita, H. Asama, and J. Park. Graph neural ordinary differ- 2019.
ential equations. arXiv preprint 1rIXiv:911..7555
M. Poli, Massaroli, Berto, J. Park, T. Dao, C. Re, and S. Ermon. Transform once: Efficient operator learning in frequency domain. 2022.
A. Quarteroni and A. Veneziani. Analysis of geometrical multiscale model based on the coupling of ode and pde for blood flow simulations. Multiscale Modeling & Simulation, 1(2):173-195, 2003.
A. Gonzalz, J. Godwin, T. Pfaff, R. Ying. J. Leskovee, and P. Battaglia. Learning to sim- ulate complex physies with graph networks. In International Conference on Machine Learning, pages 8459-8468. PMLR, 2020.
V. Strassen et al. Gaussian elimination is not optimal. Numerische mathematik, 13(4):354-366, 1969.
M. Takamoto, T. Praditia, R. Leiteritz, D. MacKinlay, F. Alesiani, D. Pflüger, and M. Niepert. Pdebench: An extensive benchmark for scientific machine learning. Advances in Neural Infor- mation Processing Systems, 35:1596-1611, 2022.
H. Tuy, A. Migdalas, and P. Värbrand. A global optimization approach for the linear two-level program. Journal of Global Optimization, 1(23). 1993.
M. Wang, D. Zheng, Z. Ye, Q. Gan, M. Li, X. Song, J. Zhou, C. Ma, L. Yu, Y. Gai, T. Xiao, T. He, package
G. Karypis, J. Li, and Z. Zhang. Deep graph library: A graph-centric, for graph neural networks, 2020.
Y.TT. Wu and F. Porté-Agel. Large-eddy simulation of wind-turbine wakes: evaluation of turbine parametrisations. Boundary-layer meteorology. 138(3)345-3366, 2011.
L.-P. Xhonneux, M. Qu, and J. Tang. Continuous graph neural networks. In H. D. III and A. Singh, editors, Proceedings of the 37th International Conference on Machine Learning. volume 119 of Proceedings of Machine Learning Research, pages 10432-10441. PMLR, 13-18 Jul 2020. URL https: / htt:///oooeeddnn..cc mlr.. press/ /v19// xhonneux20a. html.
O. Yadan. Hydra- framework for elegantly configuring complex applications. Github, 2019. URL https: thub. womm///aebboorreercc////aae/eeee
Y. Yin. V. Le Guen, J. Dona, E. de Bézenac, I. Ayed, N. Thome, and P. Gallinari. Augmenting physical models with deep networks for complex dynamics forecasting. Journal of Statistical Mechanics: Theory and Experiment, 1121
13

---

Learning Efficient Surrogate Dynamic Models with Graph Spline Networks Supplementary Material
Table of Contents
A Additional OSC Material
14
A.1 I-D OSC Example
14
A.2 2-D OSC Example
15
A.3 Simple Numerical Example
16
A.4 Numerical Analysis for interpolation and collocation methods
16
B
Additional Experimental Details
18
B.I
Models and implementation details
18
B.2
Heat Equation
19
B.3
Wave Equation
19
B.4
2D Incompressible Navier-Stokes
20
B.5 Ocean Currents
20
B.6
Black Sea
20
B.7 3D Compressible Navier-Stokes Equations
20
B.8 Hardware and Software
22
C Additional Visualizations
23
Additional OSC Material
We further illustrate the OSC method by providing numerical examples this section.
A.1 1-D OSC Example
For simplicity and without loss of generality, consider the function domain unit domain [0, 1] and we set the 1-D ODE problem. We firstly choose the partition points F;, i 0, 3, 2o 0, 1. The number of partition points is + 4. Distance between partition points can be fixed or not fixed. Then, based the Gauss-Legendrr quadrature rule, choose the collocation points. The number of collocation point within one partition is r - = 1, so we have in total N x (r (r-1) = collocation points f,, 0.....3. 3.
3, r = which means we will use a thrce-order three-piece function to simulate
After getting partition points and collocation points, we will construct the simulator. Here we have three partitions; in each partition, assign 2 order polynomial 00.0 + 00,1 r + 20 x
[Fo,
(4a)
01.0 a1, a 222
[21
(4b)
02. + 22 0222 x [ 3'3
(4c)
Note that these three polynomials should be C continuous at the connecting points, i.e., partition points within the domain. For example, Eq. (4a) and Eq. (4b) should be continuous at ø1, then
14

---

can get two equations
00.0 + d0,121 00,201
01,0 + 01,1 01,24
(5)
0 + 40,1 + 200
+ 01,1 201,23)
For boundary condition
[b, I = 20
ú(x) =
(6)
b22 I F3
we can also get two equations
00.0 + +0+0
b1
01.0 + 01,1 + 01,2
= b2
Let us summarize the equations we got so far. Firstly, our undefined polynomials have N (r+1) 9 parameters. The C² continuous condition will create (N - 1) x 2 = 4 equations, and the boundary condition will create 2 equations. Then we have x (r - 1) collocation points. For each collocation point, substitute it to polynomials to get equation. For example, if the ODE is
ú(x) + (z)) = f(z), E [0,
(8)
By substituting collocation point to into the equation, we can get ú(Eo) + i' (Eo) = f(o)
40.0 + 10 + 00 + 40,1 + 200 = f(Eo)
(9)
40.0 + a5o + +) + 25, + 25o) = f(5o)
Now we can know that the number of equations can meet with the degree of freedom of polynomials
Parameters
Boundary
Contimuous
Collocation
(10)
+ 1) x N
+ (N - x 2 N
In this example, generated equations will be constructed to an algebra problem Aa = f where the weight matrix is ABD matrix.
Not 0
Figure A.1: Visualization of an ABD matrix.
0
0
0
0
Eo
+
1
22
+
250
0
A
31
-1
-I1
(1la)
201
0
-1
22
E1
++1
5
25
0
0
1
1
G0,0
b1
Qg.1
f(En)
a
60.2
f
0
(11b)
01,0
01.1
f(L))
[0122
b2
By solving this problem, obtain the simulation results. Note that can employ the COL- ROW sparse solver to solve the system efficiently (see $ 3.4).
A.2 2-D OSC Example
For simplicity and without loss of generality, consider the function domain unit domain |0, 1]x [0, 1], and set N_ N, = 2 3. Partition points and collocation points selection is similar to the I-D OSC method; we have N2 x (r 1)2 collocation points in total. For simplicity, we note the partition points at two dimensions to be the same, i.e., Z,, i = 0, 1, 2. Unlike the 1-D OSC
15

---

method, we choose Hermite bases to describe as the simulator, which keeps C' continuous. As case, the base function point 21 would be
Hi (1( + (1)
[ (zu, z1l
f((a)
(z1-20)²
(0-z2)1z-z1)
(#2-21)²
E (I1, ز25
(12)
(2) -2) 20)- 2 - m0>² 1 [10,
(1
(z1-Bo) (22-2)) -
(x1,
(#2-#1)
We separately assign parameters to basis functions, i.e. Hz = f1 (2) (I) for variable in (mo, _1] [ Vil partition. Then the polynomial in a partition is the multiple combinations of base functions of two dimensions. For example, the polynomial in the partition [zo, [Vo. vi is
fo(x) + 18 (x) + ff) + 6f (1)1
(13)
x I fo(v) + bb,190 (9) + ai.1 (v) + (v))
Now consider the freedom degree of these polynomials. From definition, 1) 24 parameter. Considering boundary conditions, we have 24 - 4 x N = 16 parameters. The number is equal with collocation points N²2 (r- which means we can get an algebra equation by substituting collocation points. Solving this equation, we can get the simulator parameters.
have 2n(r- (n
We similarly have multiple basis functions and parameters the simulation result for the higher dimension OSC method. And then, select partition points and collocation points using the same strategy the 2-D OSC method. The rest algebra equation generating and solving equations parts will not be different.
A.3 Simple Numerical Example
We set N = 3,r = 3 to simulate the problem
u + u' = sin() + (rros((2az) u(0) = u(1)
(14)
can get a simulation solution as follows, which is visualized in Fig. A.2.
6..2x 0.4x² 31.4x x E [0,1/3)
ū(x)= = 1.5 + 1.6z 13.8x + 1x", 2 E [1/3, 2/3)
(15)
28.5 - 100x 108.5x² - 37. E 12/3,1]
1.0
0.5
0.0
-1.
Simulation
Ground Truth
Partition Paints
Collocation Points
Figure A.2: Visualization of an 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 of problems: 1-D linear, 1-D non-linear, 2-D linear, and 2-D non-linear problems. In these experiments, we 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 can directly find the real solution. For non- non-lincar problems, increasing the order of the simulator would be ideal way to get lower
16

---

loss. For example. we set the order of the simulator to 4 for 1-D non non-linear problem and for 2-D non-linear problem. Thanks the efficient calculation of OSC, even though we use higher-order polynomials to simulate, we use less running time to get results.
Table A.1: Error OSC and interpolation methods on different PDEs problems: u(z) x² 2x²² 1.16x² 0.16x (1-D linear), u(x) sin(3xx) (I-D non-linear) u(x, y) r²y² xt ry + xy (2-D linear), u(1 y) = sin (3 (3 (3y (2 2-D non-linear).
MODEL
LINEAR
NON-LINEAR
LINEAR
NON-LINEAR
NEAREST INTERPOLATION
2.3670x
1.7558x 10
1,9882x 10
3.8695 10
LINEAR INTERPOLATION
1,8928 10
8,7731 10
3,4317x10
1.1934X 10
CUBIC INTERPOLATION
3,5232x 10 12
2.2654
2.9117x10
4.5441
1.7239x10
3.4462x10
17

---

B Additional Experimental Details
B.1
Models and implementation details
In this section we introduce additional experimental details.
Model Specific details of model components are introduced below:
Component
Description
Configuration
MGN encoder
three-layer MLP
hidden size= 64
MGN processor
in total 3 processors, each three-layer MLP three-layer MLP
hidden size= 64
MGN decoder
hidden size= 64
Table B.I: Model Components and Configurations.
All MLPs have ReLU I + max(O, x) nonlincarities between layers.
Specific details of applying the OSC method are introduced as follows:
Configuration
Time-oriented OSC
Space-oriented OSC
Polynomial Order
3
Collocation Points
2
2
Table B.2: OSC Configurations.
We note that the number of collocation points in one partition is the same for each spatial dimension; features also share collocation points (i.e., wind speeds and temperature measurements share the same partition).
Regarding the choice of collocation points, for regular boundary problems, initialize colloca- tion points using Gauss-Legendre quadrature, following the recommendation by Bialecki (1998)². This choice is particularly apt for smooth and well-behaved functions over the integration interval, making it suitable for achieving accurate results without oscillatory behavior. In irregular bound- ary scenarios, employ mesh points collocation points. Importantly, the space-oriented OSC technique enables us to obtain values at any spatial position, rendering the entire space available for collocation points, which makes the adaptive collocation points possible.
Training A batch size of 32 was used for all experiments, and the models were trained for up to 5000 epochs with early stopping. We used the Adam optimizer (Kingma and Ba, 2014) with an initial learning rate of 0.001 and step scheduler with 0.85 decay rate every 500 epoch. For all datasets, we used the split of 5 1 for training, validating, and testing for fair comparison. We set the parameter B for adaptive OSC as constant. In particular, we set B = 2 x min partitions width. This ensures that collocations will not shift by more than half of the partition range in one adaptation step, lending stability and robustness to the adaptive mesh.
Space and Time Resolutions We additionally add a more detailed overview of the experiments regarding data setup, especially for space and time resolution, such as total rollout time, space, and time resolution for both OSC space and time in Table B.3.
Note that the rollout time and the other times At refer to the original dataset prediction horizons. For instance, in the Black Sea dataset, the ground truth data provided every day d), while the prediction horizon consists of 10 steps (i.e., total of 10 days). For datasets including more than single scalar value (Ocean Currents and Black Sea have features, while 3D Navier-Stokes have 5), we the OSC is applied feature-wise.
Note that other choices may be made. For instance, Chebyshev grids may offer advantages such as spec- tral convergence properties and the ability to handle function singularities more effectively. This can be an interesting avenue for future research.
18

---

Total Rollout Time
Step (At)
Resolution
Experiment
Ground Truth
Collocation
Ground Truth
Collocation Points
Heat Equation
5
0.1 s
0.2
64 x
12 x 12
Damped Wave
10
s
2
64 x 64
32 x 32
2D Navier-Stowd
10 s
Is
2
64 x
32 x 32
Ocean Currents
10 h
1 h
2 h
73 X 73
36 x 36
Black Sea
10
1 d
2
5000 nodes
1000 nodes
3D Navier-Stokes
0.5
0.05
0.1
262, nodes
32,768 nodes
Table B.3: Experiment details and parameters each dataset.
Example Rollout To further elucidate the rollout procedure with GRAPHŠPLINENETS show- case simple example. For the 2D Navier-Stokes experiment, the ground truth trajectories consist of 10-time steps (10 seconds), with step size of At 1 second. In our model, the neural network au- loregressively performs the rollout 5 times with At = 2 seconds, and the OSC (Orthogonal Spline Collocation) method used interpolate through the skipped steps. In other words, the ground truth consists of
Y, (1, 4, 7,8, 10). while the neural networks in our model provide the interpolated outputs as
Y Imsopian:: iiterpolatee {2, 4, 6, 10),
which are then further refined using the Time-Oriented OSC to produce
Y, TimeOrientedOOSC( mmpppte)
B.2 Heat Equation
The heat equation describes the diffusive process of heat conveyance and can be defined by
Au
(16)
òt
where u denotes the solution to the equation, and A is the Laplacian operator over the domain. In n-dimensional space, be written as:
Pu u
Au
(17)
Ox
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 graph neural network then uses the mesh for training.
B.3 Wave Equation
The damped wave equation can be defined by
8w
Ow
Aw =0
0t
0t
where e is the wave speed and & is the damping coefficient. The state is X (w,
Data generation We consider a spatial domain represented as a x grid and discretize the Laplacian operator. A is implemented using 5 x 5 discrete Laplace operator in simulation; null Neumann boundary condition are imposed for generation. We set c 330 and & 50 similarly to the original implementation in Yin et al. (2021).
19

---

B.4 2D Incompressible Navier-Stokes
The Navier-Stokes equations describe the dynamics of incompressible flows with 2-dimensional PDE. They can be described in vorticity form as:
aw
Øt = vVw + vnw +
(18)
Ve = 0
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 of At = 1 and a time horizon of = 10. We use similar settings as in Yin et al. (2021) and Kirchmeyer et (2022): the space discretized on a 64 x 64 grid and we set (x, v) 0. 0.1(sin(2m(r + v)) + cos(2 ( + v))), where x, are coordinates on the discretized domain. We use viscosity value v = 10-3
B.5
Ocean Currents
The ocean currents dataset (Nardelli et al., 2013; Marullo et al., 2014) is composed of hourly data of ocean currents with an Earth grid horizontal resolution of 1/12 degrees with regular lon- gitude/latitude equirectangular projection". We employ data from 2019-01-01 to 2022-01-01 at depth of 0.49 meters. We consider an area in the Pacific ocean characterized by turbulent currents corresponding to latitude and longitude (-126 -20,, -36 -30). We use currents castward and northward sea water velocities as variables, i.c., uo and vo, respectively as well as water tempera- ture to.
B.6 Black Sea
The Black Sea dataset is composed of daily real-world measurements of ocean currents and tem- peratures (Ciliberti al., 2021). The resolution of the raw data is 1/27 1/36'. We employ data starting from 01/01/2012 until 01/01/2020; we split training, validation and testing with ra- tios of 5:1:1 in the other datasets in the paper sequentially model has not seen data from 2018 or 2019 during training. Of the data points in the grids, less than 50% actually cover the Black Sea due the irregular shape of the sea. To obtain a mesh, we subsample the grid using Delaunay triangulation to 5000 nodes for baselines and 1000 nodes for GRAPHSPLINENETT, which results in a non- uniform mesh with an irregularly-shan shaped boundary akin to Lienen and Günnemann (2022). We can observe detail of the ground truth, hon-adaptive and adaptive meshing alongside GRAPHSPLINENETS errors in Figs. B. B.3 respectively. We use currents eastward and northward sea water velocities as variables, i.e., uo and vo, respectively, as well as water temperature to at single depth of 12.54m. We normalize features by the mean and standard deviation of the training samples. The dataset is available to download from the Copernicus Institute".
the temporal axis; i.e., that the
B.7 3D Compressible Navier-Stokes Equations
These equations describe the dynamics of compressible fluid flows (Takamoto et al., 2022);
0,p+ + V (pv) = 0,
(19a)
p(0 + Vv) Vp + + (/+ v).
+
(19b)
0 + 2
[
po? o v-v 0
+
c+p+
(19c)
where p is the mass density, v is the velocity, p is the gas pressure, = (T 1) the internal energy, I' = 5/3, o' is the viscous stress tensor, 77 the shear and Ç the bulk viscosities. We employ a downsampled version of the dataset from PDEBench of Takamoto et al. (2022), originally 128 x 128 128 and subsampled to 64 64 x 64 up to N₂ 10 for the ground truth data, while the
Available at https://data
http://daaa.mmiin..... marine www.:/aa.......ccuu.....c http:icuucccuuooouu//////GG .ww...BAA..A.. AAAL_ANALLSISFORRCCC. TRAICCASSTPH______1111_ Available at https://daaa http:://dttammariee..m www.mma..araaiiiiiiiio..... http://pww.o.i..e//ooo/////LLLS. MULTIYEAR PHY 200_000000000 004/esscrriiionnnn
20

---

Wind Speed Ground Truth
0.00
0.2
0.3
0.5
0.6
Figure B.1: Wind speed ground truth for the Black Sea datasct.
Without Adaptive Collocation Points
0.00
0.05
0.10
0.15
0.20
0.25
0.30
Mean Squared Error (max: 0.43)
Figure B.2: MSE for GRAPHISPLINENTTS predictions on the Black Sea dataset without adaptive meshing.
21

---

With Adaptive Collocation Points
0.00
0.05
0.10
0.15
0.20
0.25
0.30
Mean Squared Error (max: 0.39)
Figure B.3: MSE for GRAPHSPLINENETS predictions with adaptive meshing. Adaptive meshing allows for capturing more dynamic regions better, evidenced the lower worst emror.
collocation points lay on dynamic grid of 32 x 32 x 32 = 32768 nodes. We use total of 600 trajectories which we divide into training and testing with 90% split done in PDEBench.
B.8 Hardware and Software
Hardware Experiments were carried out on a machine equipped with an INTEL CORE 7900X CPU with 20 threads and NVIDIA RTX A5000 graphic card with 24 GB of VRAM.
Software Software-wise, 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. We employ solid starting point for reproducible deep learning even outside of supervised learning (Berto et al., 2023); the Lightning- Hee----ydrrllliiee combines the PyTorch Lightning library (Falcon et al., 2019) for efficiency with modular Hydra configurations (Yadan, 2019).
main template for our codebase the Lightning-Hyeee
Template which we believe
22

---

C
Additional Visualizations
Ground Truth
GCV
GUN+MLPP
MGN
PGN-
MGN+08C MGN-OSCCAAsptiee
2012
Fianre C GLN-MLP
Vienalization of the Damned Wave datacet GEN
Ground Truth
GCN
MGN
MGN (SC(est
MIG--O OSC MGN OSCCAdattee))
85
20
cB
Figure C.2:
Visualization of the Navier Stokes dataset. [-4
t-3
5-5 t t-7
t-8
Ground Truth
GCN
0.6
GEN
0.4
0.2
Our
0.0
Figure C.3: Wave dataset prediction results. GRAPHSPLINENETS generates more stable and smoother predic- tions compared to baselines.
23

---

Ground Truth
088
GCN
006
-04
GEN
Our
Gounn Trush
0.8
GCN
GEN
-02
Our
Figure C.4: Ocean dataset prediction results. GRAPHSPLINENETT generates more stable and smoother predic- tions.
24

---

1
1
GCN
GEN
Our
-06
CCN
GEN
-02
GCN
Figure C.5: Ocean dataset prediction error. GRAPHSPLINENETS contains the error better than baselines.
25


Reviews:
[{"review_id": "0KR31z08qQ", "reviewer": "Reviewer_B9pj", "summary": "The presented manuscript proposes a novel architecture to predict the evolution of the dynamical system. The proposed method has two key ingredients: orthogonal spline collocation for both time and spatial dimensions and adaptive selection of the collocation points. Extensive experiments demonstrate the performance of the proposed method over alternatives in both prediction quality and runtime.", "questions": "1) how to select the parameter $\\beta$ in eq (3) ? Is it constant or adaptive too?\n2) do you observe any relationship between the number of processor blocks and the resulting error?\n3) why do you not consider the Chebyshev grid for collocation points and the corresponding barycentric Lagrange interpolation approach? Cheb grid should be more accurate than your static approach (Fig. 4.3 left) but less computationally intensive than the adaptive approach.\n4) what item in the loss function $L_s$ or $L_i$ dominates and what are the possible explanations?", "limitations": "The authors provide the limitations of the proposed approach.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 1, "contribution": 3, "strengths": "The presented method is based on the natural encoder-decoder architecture and intermediate processor via message-passing layers. The main strength of the paper is the extensive experimental evaluation and sensitivity analysis of the results to the adjusting hyperparameters. The authors consider 3 datasets generated from different types of PDEs and 2 datasets with empirical data from monitoring Ocean and black sea. Their method shows the best results in all considered tests. Moreover, the proposed approach is faster than the competitors based on graph neural networks.", "weaknesses": "The main weakness of the paper is its structure and presentation of the proposed method. In particular,\n1) Too much attention is paid to the standard method for solving structured linear systems for computing interpolation coefficients. I suggest compressing such detailed descriptions since they are mostly textbook facts. At least a comparison with a general-purpose LU solver from PyTorch is strange since the result of such a comparison is obvious. Thus, section 3.4. can be eliminated or significantly compressed.\n2) section 3.2 introduce the notation for the graph, but later there are no details on how this graph is constructed from snapshots\n3) in the paragraph in lines 133-139 no information about how to select collocation points only the number of collocation points is estimated. So, it is unclear how to extract collocation points from the available samples\n4) Figure 4.1 is entirely confusing. I do not find where is input, output, or ingredients of eq (1) and how they should be processed together to get the prediction. Please, specify what operations are done and in what order starting from the beginning and revise section 4.1 respectively.\n5) also, please move to the main text architecture details of the encoder, decoder, and processor blocks and training setup from the supplementary materials  \n\nI will increase the score if the proper revision of the text will be made in the rebuttal phase.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the constructive review and feedback. We provide responses to your concerns and questions below.\n\n### Too much attention to standard methods\n\nThank you for your suggestion. We agree with you that certain sections detailing the standard method for solving structured linear systems can be condensed. However, the emphasis on our comparison with the general-purpose LU solver from PyTorch was intentional since our aim was to highlight the novelty and efficiency of our implementation: a fully differentiable COLROW method with GPU acceleration which, to the best of our knowledge, is the first of its kind and furthers our contribution. Given its significance, especially in the context of the SciML community, we deemed it beneficial to retain this section. Nonetheless, we will refine it to make it more concise while still preserving its essence based on your suggestion.\n\n### How is the graph constructed from snapshots?\n\nThanks for spotting the missing explanation; we will thoroughly revise this part.  We construct the graph based on the mesh on the snapshot. In regular boundary problems, we create structured mesh based on collocation points. In irregular boundary problems, for example, in the Black Sea dataset, we use the Delaunay Triangulation to create the mesh. \n\n### How do we select collocation points?\n\nFor regular boundary problems, we initialize collocation points using Gauss–Legendre quadrature, following the recommendation by B. Bialecki [3]. This choice is particularly apt for smooth and well-behaved functions over the integration interval, making it suitable for achieving accurate results without oscillatory behavior. As outlined in Section 3.3, the requirement primarily revolves around the number of collocation points. In irregular boundary scenarios, we employ mesh points as collocation points. Importantly, the space-oriented OSC technique enables us to obtain values at any spatial position, rendering the entire space available for collocation points, which makes the adaptive collocation points to be possible. We incorporate these clarifications in Section 3.3.\n\n### Confusing Figure 4.1 and missing architecture details\n\nThank you for the valuable comment! We have revised the figure to depict the flow from input to output more accurately and to integrate the components of eq (1). This updated figure, available on the PDF, now offers a clearer overview of the method and elucidates the inner workings of the neural network. Based on your suggestion, we also included a new supporting figure showing the encoder-processor-decoder structure ”Neural Network (NN)” in the new figure. We will revise Section 4.1 based on your suggestions, including cleaning up the notation and better clarification on the overall process, such as inputs and outputs! \n\n### About $\\beta$\n\n$\\beta$ is constant. In our experiments we set $\\beta$ to be $\\beta = \\frac{1}{2}\\times \\mathrm{min\\ partitions\\ width}$. This ensures that in one adaptation step, collocations will not shift by more than half of the partition range, lending stability and robustness to the adaptive mesh. We will include this clarification in the appendix.\n\n### Error and number of processor layers\n\nWe conducted more experiments to showcase this trend in Table R1 and Figure R2 in the attached PDF. Generally, with an increase in processor blocks, there is a trend of decreased error at the cost of higher runtime, reflecting findings from prior work such as GNS [1] and Meshgraphnets [2]. For the sake of a fair comparison, the baseline processor remained consistent across all models. Notably, GraphSplineNets offers a Pareto optimal trade-off, outperforming Meshgraphnets (MGN) in terms of both prediction accuracy and computational efficiency.\n\n### Using the Chebyshev grid for collocation points\n\nYour question is on point. Figure 4.3 shows a regular grid for illustrative purposes, but this is not the only possible choice. As we mentioned in the previous reply about the approach to selecting collocation points, the Gauss–Legendre quadrature is particularly for functions that exhibit smooth behavior over the interval of integration. It excels in scenarios involving finite intervals, offering accurate results without introducing oscillatory behavior. However, Chebyshev grids may offer advantages such as spectral convergence properties and the ability to handle function singularities more effectively. This can be an interesting avenue for future research.\n\n### Effect on sample points and interpolation points losses\n\nTable 5.1 shows an ablation study on these losses. The MGN+OSC(Post) works in a post-processing way, so it doesn’t utilize the $L_i$; while the MGN+OSC uses both losses and enhances the performance. We can see that while $L_s$ works as a primary contributor to the observed improvements, $L_i$ also plays a meaningful role, enabling the model to encompass a broader spectrum within the feature space and enhance its ability to discern intricate patterns. We will clarify this in the manuscript!\n\n### Final note\n\nThank you for considering the adjustment of the score based on the revisions! However, as per this year's NeurIPS guidelines, modifications to the paper text during the rebuttal phase are not allowed. As such, we kindly request you evaluate our submission based on the clarifications provided in our responses above. We're grateful for your understanding and constructive feedback.\n\n### References\n\n[1] Sanchez-Gonzalez, Alvaro, et al. \"Learning to simulate complex physics with graph networks.\" *International conference on machine learning*. PMLR, 2020.\n\n[2] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia. Learning mesh-based simulation with graph networks. *International Conference on Learning Representations*, 2021.\n\n[3] Bialecki, Bernard, and G. Fairweather. \"Orthogonal spline collocation methods for partial differential equations.\" *Journal of Computational and Applied Mathematics* 128.1-2 (2001): 55-82."}, {"review_id": "X7hwsTOrH4", "reviewer": "Reviewer_3ehZ", "summary": "The authors propose a method based on deep learning that will speed up the prediction of physical systems by reducing the size of the grid and the number of tests to be carried out. The method consists of combining a graph neural network to predict some values with two orthogonal spline collocation, one spatial and one temporal, to predict values continuously.", "questions": "1) Lines 176-177: I think there is a typo, 2-D domain is written twice\n2) Line 203 and equation (2): you write $L$ in the equation and then you write a rounded \"L\", I think you mix the notation\n3) To what extent is the availability of data for physical systems not too restrictive? \n4) Could we imagine replacing GNNs with physical neural networks?", "limitations": "The impact of GNN on the overall method needs to be quantified a little more, the level of performance required by the neural network needs to be defined. This would make it easier to adapt the method, for example by changing the GNN model for another one.\nWe can make the same point with OSC, trying to combine GNN with an alternative to show what the use of OSC brings specifically, in terms of runtime and performance.\nThis would help justify the choice of these two methods over others.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is very clear and well written. The graphics are a great complement to the text and help to make it easier to understand. \nThe choice of graph neural networks seems judicious in the context of a physical system that admits links and dynamic interactions between variables. As for orthogonal spline collocation, the theoretical guarantees and computational complexity justify their use. The results show improved predictions and increased speed on various physical problems. The authors are careful to point out the limitation of their approach.", "weaknesses": "See limitations and questions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and for acknowledging the strengths of our work in soundness, presentation, and contribution! We address your concerns point by point in the following sections.\n\n### Availability of data for physical systems\n\nYour question delves into an important consideration. The availability of data for physical systems can be task-dependent, a common trait shared with various deep-learning models. Our approach yields superior results with fewer sample points than conventional surrogate models. The adaptability of our approach contributes to the observed robustness and accuracy improvements. The extent of data availability hinges on factors such as task scale and accuracy requirements. In a practical context, our collocation points can be likened to sensors. Furthermore, the adaptive collocation points mimic movable sensors, which is a practical arrangement in which a finite number of sample points has to be adjusted to maximize the amount of information. The scalability of the system influences sensor demand—the larger the scale, the greater the need. This flexible sensor deployment addresses varying levels of data availability, offering enhanced accuracy while optimizing resource usage.\n\n### About replacing GNNs with physical neural networks\n\nIf you mean physical neural networks in the sense of adjustable hardware-based networks: we believe it can be an interesting area for future research for our method that could enable way faster runtimes! However, let us also answer in case you meant physics-informed neural networks (PINNs). In this case, the objective would be different for PINNs and GraphSplineNets: in PINNs the scope is, given an equation, we fit parameters to satisfy its conditions; in our case - deep surrogate models - we optimize parameters to fit given data, which may or may not be easily described by an equation.\n\n### Impact of GNN discussion\n\nWe appreciate your valuable feedback. The GNN plays a pivotal role as the foundation of our methodology, enabling us to effectively consider interactions between various regions while accommodating irregular meshes. It serves as a fundamental tool to capture intricate relationships and patterns within the data. For scenarios involving regular domains and boundaries, alternative models such as the Fourier Neural Operators (FNO) could be explored for space-based collocation points. Meanwhile, we indeed have flexibility in model selection for the time aspect. For instance, non-graph models could be employed, broadening the applicability of our approach.\n\n### About typos\n\nThanks for spotting the typos! We will fix them in the manuscript."}, {"review_id": "DVq8nQNsCG", "reviewer": "Reviewer_jrMD", "summary": "The authors introduce GraphSplineNets, providing a framework for combining MPNNs like MGNs with interpolation/collocation methods. Their contributions also include two methods to surpass the performance of a naive combination: end-to-end training and collocation points adaptation. Temporal and spatial downsampling allow the MPNN to run faster, while the collocation method makes possible *increased* inference-time resolution despite the NN's downsampled domain.  ", "questions": "1. Why use only 3 MP layers in baseline MGNs? This gives the GraphSplineNets method an unfair advantage in that it can communicate over a larger spatial extent (because it is operating on a coarser mesh). Relatedly, the baseline MGN results (e.g., those in table 5.1) are rather poor compared to results shown in the MGN paper, which recommends using 15 MP layers.\n\n2. If you keep on increasing the collocation points, will GSNs eventually provide the same solution as MGNs? Assuming this is the case, Figure 5.3 suggests that more collocation points is better, but it might be nice to see when the relationship flips and the methods converge.\n\n3. How do GSNs perform on larger scale problems? In particular, can they compete with approaches designed to be more scalable like MultiScale MeshGraphNets (Fortunato et al., 2022)?\n\n4. Line 27: It's true that traditional solvers can struggle, but it might make sense to support this point by citing a paper that analyzes traditional solvers rather than a paper that analyzes ML solvers.\n\n5. In the Appendix, there seem to be some minor typos/errors in the equations. E.g., equations 7 and 11 seem to assume 2 partitions when the example has 3 partitions. This did not affect my score, but cleaning it up would be great!\n", "limitations": "The authors addressed some key limitations. I might add that no large-scale or 3D problems were considered, which hinders understanding of the benefits of the proposed method at the scales where its contribution (efficiency) might matter most.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Significance:\n\n1. This is an important area, and the paper's core idea is interesting and well presented (as a proof of concept, at least). The paper offers a lot of infrastructure (code and ideas) important to future steps at the intersection of GNN surrogates and collocation.\n2. The inference speed benefit of GSNs is intriguing because applying GNNs at scale can be difficult.\n\nClarity:\n\n3. The motivation for this work is communicated well. This work's connection to other approaches that deal with GNN scalability is clearly established. \n4. The explanation of the GNN-based OSC method is very well done and nicely illustrated. \n\nQuality:\n\n5. OSC seems like a potentially promising approach for dealing with the efficiency challenges of GNN simulations.\n6. The end-to-end optimization approach (differentiating through the OSC) is well justified and intuitive. \n\nOriginality:\n\n7. Plugging GNNs into an OSC method, allowing temporal coarsening of MPNN (e.g. MGN) activity to increase efficiency, and optimizations of this joint GNN+OSC method (end-to-end training, adaptive points).\n8. Giving automatic differentiation support to an ABD matrix equation solving routine. \n", "weaknesses": "My main concern is the quality of the evaluation of the method. I think this paper would gain significance if it showed that GSN's performance benefits hold in more challenging settings (by demonstrating why the benefits of GSN will generally be present, by evaluating on harder problems and against better tuned baselines, etc.). In other words, addressing some of the following weaknesses might help other researchers better understand the potential relevance of this method.\n\n1. The paper leaves unclear why its methods beat baseline methods' errors (i.e., is it because downsampled mesh nodes have access to information from farther away?), if the errors it achieves are practically relevant (i.e., sufficiently low for the model to be used), and if it beats stronger baselines (i.e., MGN with more message passing layers or MultiScale MGN). \n\n2. Moreover, no large-scale or 3D problems were considered, and this is at odds with the description of the proposed method as one that addresses efficiency challenges of existing approaches.\n\nMinor concerns: \n\n3. Figure 4.1 is nice, but I think at $t=0$, you should have a $y_3$ instead of a $y_0$. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the constructive and detailed review and thoughtful feedback! We provide responses to your concerns and questions below.\n\n### Why do GraphSplineNets work well?\n\nThank you for your question! This is an important point for our paper and as such we would like to clarify it. \n\nAs you correctly noticed, we empirically found that downsampling mesh nodes can help in robustness and generalization; the OSC plays a crucial part in handling the problem’s continuity (e.g. in time) while adaptive collocation effectively helps the model to focus on some regions with higher errors.  This can be in part attributed to the fact that GraphSplineNets can communicate to a larger spatial extent; but not only, as demonstrated in our additional experiments.\n\nWe include in the PDF a new figure Figure R2 and table Table R1 that show a comparison of GraphSplineNets against MGN with 3 message passing layers and 15. As we can notice from the comparison, a higher number of message-passing layers generally helps in decreasing the MSE at the cost of higher runtime, reflecting findings from prior work such as GNS[1] and Meshgraphnets[2]. For the sake of a fair comparison, the baseline processor remained consistent across all models. Notably, GraphSplineNets shows a Pareto optimal trade-off, outperforming Meshgraphnets (MGN) in terms of both prediction accuracy and computational efficiency at different numbers of message-passing layers. Given that the neural network structure is the same for both our GraphSplineNets and MGN, we attribute the increase in accuracy (and importantly, at a reduced computational cost) to our proposed method.\n\n### Large-scale and 3D experiments\n\nWe agree with you! While we focused mostly on simulated and real data from 2D domains in the main paper, given our background in fluid dynamics and weather forecasting, we do agree that large-scale 3D experiments can better demonstrate the scalability of our approach. As such, we have prepared a new experiment - namely, the 3D Navier-Stokes equations dataset from PDEBench [5]. We use the 3D Navier-Stokes dataset with a total of 600 trajectories which we divide into training and testing with a 90% split as done in PDEBench. The number of nodes in GraphSplineNets is 4096, which is more than 4 times higher than in our previous hardest experiment (the Black Sea Dataset). Moreover, the 3D Navier Stokes include a total of 5 variables: density, pressure and velocities on 3 axes. For fair comparison, we keep the model structure the same as in other experiments.\n\nWe show in the attached PDF in Figure R3 and Table R2 that our GraphSplineNets can successfully scale to complex, large-scale, 3D data in both accuracy and importantly speed - compared to the MGN baseline, we were able to obtain around 10x increase in inference speed! \n\n### Effect of increasing the number of collocation points\n\nThanks for your question. Increasing the collocation points brings GSNs closer to the behavior of  MGNs. When the number of collocation points matches the sample points, only the time-oriented OSC works. Conversely, if the neural network in GSNs predicts every steps, only the space-oriented OSC in GSNs. In Figure 5.3, we show that more collocation points typically yield lower errors, but this will cause an increase in computational time. It's also imperative to recognize that the extended rollout accuracy in Figure 5.3 for MGN+OSC is attributed to the incorporation of the time-oriented OSC.\n\n### Citation for traditional solvers\n\nThanks for your suggestion. We will cite a paper analyzing traditional solvers instead. We are thinking about adding [3, 4]; however, if you have a better citation(s) in mind, we would be glad to hear more about it.\n\n### Typos\n\n- Figure 4.1: we are impressed by your attention to detail! Thanks, we have updated the figure, which can be retrieved on the rebuttal PDF.\n- Equations 7 and 11: thanks, we will clean up the equations based on your and other reviewers feedback!\n\n### References\n\n[1] Sanchez-Gonzalez, Alvaro, et al. \"Learning to simulate complex physics with graph networks.\" *International conference on machine learning*. PMLR, 2020.\n\n[2] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia. Learning mesh-based simulation with graph networks. *International Conference on Learning Representations*, 2021.\n\n[3] Kremer, D. M., and B. C. Hancock. \"Process simulation in the pharmaceutical industry: a review of some basic physical models.\" *Journal of pharmaceutical sciences* 95.3 (2006): 517-529.\n\n[4] Oberkampf, William L. \"Simulation accuracy, uncertainty, and predictive capability: A physical sciences perspective.\" *Computer Simulation Validation: Fundamental Concepts, Methodological Frameworks, and Philosophical Perspectives*(2019): 69-97.\n\n[5] Takamoto, Makoto, et al. \"PDEBench: An extensive benchmark for scientific machine learning.\" *Advances in Neural Information Processing Systems* 35 (2022): 1596-1611."}, {"review_id": "p84OoH00W0", "reviewer": "Reviewer_Rs2G", "summary": "The paper builds upon a recent advances in deep-learning based collocation methods for physical simulation to improve prediction accuracy in continuous space and time. The authors proposed a hybrid model based on graph neural network and differentiable orthogonal spline collocation method. Interpolation in time as well as in space is adopted to enable employing coarser temporal grids to speed up prediction. The paper also introduces a new adaptive collocation sampling strategy and techniques for training such as not needing to learn interpolators or collocation weights to decrease the computational overhead. Experiments in various PDE settings demonstrate that the proposed model, GRAPHSPLINENETS outperforms state-of-the-art baselines in both accuracy and runtime.", "questions": "* Is OSC inference time in Figure 5.1 associated with both spacial and time oriented collocation? \n* Does adaptive collocation sampling also perform for important time interval?", "limitations": "Yes, the limitations are discussed in the main text. ", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper introduces some techniques, a novel loss function comprising OSC function, COLROW algorithm to efficiently solve linear equations associated with OSC methods, and adaptive collocation sampling to allow for prioritized sampling of important location regions. Wide variety of experiments are conducted, and the proposed method outperforms the baselines.", "weaknesses": "Although the authors claim that GraphSplineNets forecast continuous responses in time and space, comparison of GraphSplineNets are limited to space continuous baselines. Also, as in the main text, only a fewer points are sampled as collocation points in time — It is still difficult to assess the advantage of incorporating time-oriented OSC as well as space-oriented into the architecture. It is also unclear how is dataset for interpolation points error in the equation (2) generated. \n\nThe explanation on Figure 5.1 is too brief for me to understand. Does the size of input of MGN is larger than that of MGN + OSC (12x12)? In appearance it looks strange to me that the running time of MGN w/o OSC higher than MGN w/ OSC. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review and for acknowledging the strengths of our work in soundness, presentation, and contribution! We address your concerns point by point in the following sections.\n\n### Why is collocation in time important?\n\nYour observation regarding the limited sampling of collocation points in time is on point. The incorporation of time-oriented OSC aims to facilitate predictions over extended future intervals in a single step - the reason why fewer points are sampled in time is that too large steps could eventually hurt the solutions in data-driven settings, much like it would in a normal time-stepping of an ODE via explicit methods.  By interpolating through few time steps, the advantage of our approach is twofold: it enhances computational efficiency by reducing the number of necessary prediction steps and potentially minimizes error accumulation that can arise from more frequent, smaller time-step predictions since the network can implicitly learn to compensate for time stepping error.\n\n### How is the dataset for interpolation points error generated?\n\nWe use interpolations on available data to match equation (2) errors. Therefore, the data is not generated but “given” as a dataset (more information about datasets used in Appendix B), as is common in supervised surrogate modeling. \n\n### Runtime comparison in Figure 5.1\n\nThanks for your question! We will clarify this in the manuscript as well. The OSC inference time in this figure is associated with both spatial and time-oriented collocation. MGN without OSC has higher runtime because MGN+OSC effectively *downsamples* the domain while still capturing the dynamics thanks to the interaction with OSC. The downsampling allows for noticeably faster runtimes and more robustness thanks to the OSC, whose overhead is minor compared to the total runtime savings. \n\n### Adaptive Collocation in Time\n\nThis is a clever idea. While we did not test it in the paper, we believe it is an exciting avenue for future work. We also have some ideas we will put down in the paper regarding this point - particularly how to handle non-regular time stepping due to the adaptive collocation in time. In this case, an additional parameter should be kept track of, namely time $\\Delta t _i$ for collocation point $i$ in time; otherwise, the network would try to predict the change at a constant time step $\\Delta t$ given by the base dataset. Each collocation point would have its own $\\Delta t_i$ and the output of GraphSplineNets can be modulated by this variable to perform time stepping effectively. Thanks for your suggestion!"}]

Author rebuttal:
We thank all reviewers for providing valuable and constructive feedback on our work; we sincerely appreciate the time and effort dedicated to reviewing our manuscript. We are glad to receive such an overall positive evaluation on all aspects of the submission, including soundness, presentation, and contribution. In response to each reviewer's comments, we have provided a detailed explanation and made necessary revisions to improve the quality of our work. 

We additionally include three common responses to summarize three major points that came out during the rebuttal:

- **Scaling GraphSplineNets to large-scale datasets**
- **Why do GraphSplineNets work well?**
- **Improvements to our manuscript**

## Scaling GraphSplineNets to large-scale datasets

Based on reviewer jrMD comments, we tested GraphSplineNets on a large-scale 3D dataset, namely, the 3D Navier-Stokes equations dataset from PDEBench [1]. We show in the attached PDF in Figure R3 and Table R2 that our GraphSplineNets can successfully scale to complex, large-scale, 3D data outperforming the baseline in both accuracy and, importantly, speed. Compared to the MeshGraphNet [2] baseline, we obtained an almost 10x increase in inference speed, while still improving accuracy. We believe GraphSplineNets would be even more impactful for efficiently learning surrogate models in even larger-scale complex datasets.

## Why do GraphSplineNets work well?

This is an important point raised by the reviewers we would like to clarify.  As correctly noted by jrMD, a reason can partly be attributed to the fact that GraphSplineNets can effectively have a longer spatial context with few message-passing layers since we downsample the domain. As also suggested by B9pj, we conducted an additional study on the effect of the number of processors (i.e., message passing layers) for both GraphSplineNets and the baseline MeshGraphNets [2] shown in Table R1 and Figure R2 of the attached PDF. Notably, GraphSplineNets shows a Pareto optimal trade-off, outperforming MGN in prediction accuracy and computational efficiency at the same number of message-passing layers. We conclude that the performance is not only due to the number of message-passing layers but also to our proposed training strategy with both time and adaptive space OSC, as demonstrated in Table 5.1 with the three variants of our method.

## Improvements to our manuscript

Based on all reviewers’ feedback, especially B9pj, we made several improvements to the manuscript and figures. We revised the main Figure 4.1 in the PDF file. Per this year's NeurIPS guidelines, modifications to the paper text during the rebuttal phase are not allowed; thus, we will upload the updated manuscript incorporating your constructive feedback after the rebuttal.

### References

[1] Takamoto, Makoto, et al. "PDEBench: An extensive benchmark for scientific machine learning." *Advances in Neural Information Processing Systems* 35 (2022): 1596-1611.

[2] T. Pfaff, M. Fortunato, A. Sanchez-Gonzalez, and P. W. Battaglia. Learning mesh-based simulation with graph networks. *International Conference on Learning Representations*, 2021.

Meta-review: The paper presents a novel architecture to predict the evolution of the dynamical system. The reviewers all agree that this is a valuable and interesting result.
Decision: Accept (poster)